fitspectrum_mod.F90

Go to the documentation of this file.

MODULE fitspectrum_mod

USE parkind1  ,ONLY : jpim     ,jprb
USE yomhook   ,ONLY : lhook,   dr_hook

IMPLICIT NONE

PRIVATE
PUBLIC fitspectrum

CONTAINS

SUBROUTINE fitspectrum(KSMAX,KSPEC2G,PSPEC,KSTRUNK,PBETA0,PBETA1)



!**** *FITSPECTRUM*  - linear fit of spectrum to (n*(n+1))**P

!     Purpose.  GRIB allows for "complex packing" of spectral
!     --------  coeficients. This means that a sub-spectrum
!               is not packed and for the rest a multiplying
!               factor of (n*(n+1))**p, where n is the total
!               wave-number, is applied before packing. This routine
!               produces an optimal P (-PBETA1) based on the field
!               itself.

!**   Interface.
!     ----------
!        *CALL* *FITSPECTRUM(...)

!        Explicit arguments :   KSMAX   - spectral field truncation
!        --------------------   KSPEC2G - size of spectral field array
!                               PSPEC   - spectral field
!                               KSTRUNK - truncation of "non-packed" 
!                                         part of spectrum
!                               PBETA0  - linear fit par.  (output)
!                               PBETA1  - linear fit par.  (output)

!        Implicit arguments :  
!        --------------------

!     Method.
!     -------
!        If F(n,m) = (n*(n+1))**(-P) * G(n,m) where F(n,m) is the
!     original spectral field and n the total wavenumber the aim
!     is to minimize G in a 1 norm with respect to P.This can only
!     partially be achieved. What we do is to to compute H(n) where
!     H(n)=max(F(n,m)) with respect to m. We then perform a least 
!     square fit for the equation log(H(n))=beta0+beta1*log(n*(n+1)).
!     To ensure a better fit for the lower end of the spectrum we
!     apply an (arbitrary) weighting funtion W(n)=1./(n-kstrunk+1)
!     before fitting.

!     Externals.
!     ----------

!     Reference.
!     ----------
!        Seber, G.A.F. (1979). Linear Regression Analyses.
!        John Wiley and Sons

!        ECMWF Research Department documentation of the IFS

!     Author.
!     -------
!        Mats Hamrud *ECMWF*

!     Modifications.
!     --------------
!       Original : 94-05-05
!       Modified 95-06-06 by L.Isaksen : Reordering of spectral arrays
!       Modified 02-04-03 Y. Tremolet  : Resolution different from NSMAX
!        M.Hamrud      01-Oct-2003 CY28 Cleaning
!     ------------------------------------------------------------------


INTEGER(KIND=JPIM),INTENT(IN)    :: KSPEC2G 
INTEGER(KIND=JPIM),INTENT(IN)    :: KSMAX 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSPEC(kspec2g) 
INTEGER(KIND=JPIM),INTENT(IN)    :: KSTRUNK 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBETA0 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBETA1 
REAL(KIND=JPRB) :: ZW(0:ksmax),ZNORM(0:ksmax)

INTEGER(KIND=JPIM) :: INM, ISMIN, JIR, JM, JN
INTEGER(KIND=JPIM) :: IASM0G(0:ksmax)

REAL(KIND=JPRB) :: Z, ZEPS, ZSUM1, ZSUM2, ZWSUM, ZX, ZXMW, ZY, ZYMW, ZZ
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

IF (lhook) CALL dr_hook('FITSPECTRUM',0,zhook_handle)
inm=1
DO jm=0,ksmax
  iasm0g(jm)=inm
  inm=inm+(ksmax+1-jm)*2
ENDDO

!*       1.    BEST FIT.
!              ---------

!*       1.1   INITIAL COMPUTATIONS

ismin = kstrunk+1
zeps  = epsilon(z)

!*       1.2   SET UP WEIGHTS

zz=REAL(ksmax-ismin+1,jprb)
DO jn=ismin,ksmax
  zw(jn) = zz/REAL(jn-ismin+1,jprb)
ENDDO

!*       1.3   COMPUTE NORMS

DO jn=ismin,ksmax
  znorm(jn) = 0.0_jprb
ENDDO
DO jm=0,ksmax
  DO jir=0,min(1,jm)
    DO jn=max(ismin,jm),ksmax
      inm=iasm0g(jm)+(jn-jm)*2+jir
      znorm(jn) = max(znorm(jn),abs(pspec(inm)))
    ENDDO
  ENDDO
ENDDO
DO jn=ismin,ksmax
  znorm(jn)=max(znorm(jn),zeps)
  IF(znorm(jn) == zeps) THEN
    zw(jn) = 100._jprb*zeps
  ENDIF
ENDDO

!*       1.4   LINEAR FIT

zxmw  = 0.0_jprb
zymw  = 0.0_jprb
zwsum = 0.0_jprb
DO jn=ismin,ksmax
  zx    = log(REAL(JN*(JN+1),JPRB))
  zy    = log(znorm(jn))
  zxmw  = zxmw+zx*zw(jn)
  zymw  = zymw+zy*zw(jn)
  zwsum = zwsum+zw(jn)
ENDDO
zxmw  = zxmw/zwsum
zymw  = zymw/zwsum
zsum1 = 0.0_jprb
zsum2 = 0.0_jprb
DO jn=ismin,ksmax
  zx    = log(REAL(JN*(JN+1),JPRB))
  zy    = log(znorm(jn))
  zsum1 = zsum1+zw(jn)*(zy-zymw)*(zx-zxmw)
  zsum2 = zsum2+zw(jn)*(zx-zxmw)**2
ENDDO

pbeta1 = zsum1/zsum2
pbeta0 = zymw-pbeta1*zxmw

!     ------------------------------------------------------------------

IF (lhook) CALL dr_hook('FITSPECTRUM',1,zhook_handle)
END SUBROUTINE fitspectrum

END MODULE fitspectrum_mod