fapula.F90

Go to the documentation of this file.

! Oct-2012 P. Marguinaud 64b LFI
! Jan-2011 P. Marguinaud Thread-safe FA
SUBROUTINE fapula_fort                                 &
&                     (fa,  krep, krang, pspec, kpulap)
USE fa_mod, ONLY : fa_com, jpniil
USE parkind1, ONLY : jprb
USE yomhook , ONLY : lhook, dr_hook
USE lfi_precision
IMPLICIT NONE
!****
!        Sous-programme INTERNE du logiciel de Fichiers ARPEGE:
!     calcul de la PUissance de LAplacien "optimale" pour aplatir
!     le spectre hors de la sous-troncature.
!
!**
!    Arguments : KREP   (Sortie) ==> Code-reponse du sous-programme;
!                KRANG  (Entree) ==> Rang de l'unite logique;
!    ( Tableau ) PSPEC  (Entree) ==> Champ en coef. spectraux en entree;
!                KPULAP (Sortie) ==> Puissance de laplacien calculee.
!*
!
!     Creation
!     --------
!  Avril 1999, D. Paradis, SCEM/TTI/DEV:
!     Ce ss-programme est derive de la routine CALCOP du logiciel GRIBEX
!     du CEPMMT. Il retourne la puissance de laplacien (exprimee en un
!     nombre entier de milliemes) qui aplatit le mieux le spectre des
!     coefficients spectraux (hors sous-troncature).
!     Cette puissance de laplacien, P, sera l'exposant du facteur
!     .(n*(n+1))**P   pour les champs ARPEGE ou n est le nombre d'onde total
!     .(n**2+m**2)**P pour les champs ALADIN ou n est le nombre d'onde meridien
!                                            et m est le nombre d'onde zonal
!     Ce facteur sera applique aux coefficients spectraux hors sous-troncature
!     pour en reduire l'amplitude avant le codage GRIB.
!
!
!     Method.
!     -------
!
!     For ARPEGE case, consider F(n,m) = (n*(n+1))**(-P) * G(n,m),
!     where F(n,m) is the original spectral field and n the total wavenumber.
!     For ALADIN case (LAM), consider F(n,m) = (n**2+m**2)**(-P) * G(n,m),
!     where F(n,m) is the original spectral field, n the meridian
!     wavenumber and m the zonal wavenumber.
!     The aim is to minimize G in a 1 norm with respect to P.  This can only
!     partially be achieved.
!
!     1) For ARPEGE case, what we do is 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 function before fitting
!             W(n) = 1.0 / (n-ISTRON+1)
!
!     2) For ALADIN case, what we do is to compute H(in), where
!         H(in) = max(F(n,m))      where n and m verify: k = n**2 + m**2
!           and "in" is the rank of k among all the reached values by k:
!           we are only interested in k values which verify
!           k = n**2 + m**2, with (n,m) belonging to spectrum without
!           the non-compacted sub-truncation area. We call inmax the number
!           of these k values. "in" belongs to interval 1,inmax.
!
!        We then perform a least square fit for the equation
!            log(H(in)) = beta0+beta1*log(k).
!
!        To ensure a better fit for the lower end of the spectrum, we
!        apply an (arbitrary) weighting function before fitting
!           W(in) = 1.0 / in
!
!     Reference.
!     ----------
!
!     Seber, G.A.F. (1979). Linear Regression Analyses.
!     John Wiley and Sons
!
!     ECMWF Research Department documentation of the IFS
!
!
!     Modifications
!     -------------
!
!
!
!     Subroutine arguments
!
TYPE(fa_com) :: FA
INTEGER (KIND=JPLIKB) KREP, KRANG, KPULAP
!
REAL (KIND=JPDBLR) PSPEC(*)
!
!     Local arguments
!
INTEGER (KIND=JPLIKB) IRANGC, ITRONC, IMTRONC, INMAX, INMIN
INTEGER (KIND=JPLIKB) IKMAX, IDEB, IFIN, JN, JM, JK, IK, IN, IM
INTEGER (KIND=JPLIKB) IMLIM, JIND, IOFF, INUMER
INTEGER (KIND=JPLIKB) INIMES, ISTRON
INTEGER (KIND=JPLIKB), DIMENSION(:), ALLOCATABLE :: ITAB1
INTEGER (KIND=JPLIKB), DIMENSION(:), ALLOCATABLE :: ITAB2
!
REAL (KIND=JPDBLR) ZEPS, ZZ, ZXMW, ZYMW, ZWSUM, ZX, ZY
REAL (KIND=JPDBLR) ZP, ZBETA1, ZSUM1, ZSUM2
REAL (KIND=JPDBLR), DIMENSION(:), ALLOCATABLE    :: ZW, ZNORM
!
LOGICAL LLMLAM
!
CHARACTER(LEN=FA%JPLMES) CLMESS 
CHARACTER(LEN=FA%JPLSPX) CLNSPR
LOGICAL                  LLFATA

!**
!     1.  -  CONTROLES DES PARAMETRES D'APPEL, INITIALISATIONS.
!-----------------------------------------------------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('FAPULA_MT',0,zhook_handle)
IF (krang.LE.0.OR.krang.GT.fa%JPNXFA) THEN
  krep=-66
  GOTO 1001
ENDIF
!
istron=fa%FICHIER(krang)%NSTROF
irangc=fa%FICHIER(krang)%NUCADR
itronc=fa%CADRE(irangc)%MTRONC
llmlam=fa%CADRE(irangc)%LIMLAM
!
IF (llmlam) imtronc=fa%CADRE(irangc)%NOZPAR(2)
IF (itronc.LE.istron) THEN
  krep=-88
  GOTO 1001
ELSEIF (llmlam.AND.imtronc.LE.istron) THEN
  krep=-88
  GOTO 1001
ELSEIF (llmlam.AND.(imtronc.GT.3*itronc.OR. &
&    itronc.GT.3*imtronc)) THEN
! Il s'agit d'un garde-fou, modifiable (ne pas oublier FARCIS et FACSIM)
  krep=-114
  GOTO 1001
ELSE       
  krep=0
ENDIF
!
zeps  = 1.0e-15_jpdblr
!
IF (llmlam) THEN
  ikmax = itronc*itronc + imtronc*imtronc
!
!  Pour le cas ALADIN, il faut creer les tableaux ITAB1 et ITAB2 :
!  ITAB1(JK) = 0 sauf pour JK tel qu'il existe (n,m) appartenant au spectre
!    hors de la sous-troncature non-compactee verifiant JK = n**2 + m**2:
!    on a alors ITAB1(JK) = 1
!  ITAB2 regroupe toutes les valeurs de JK precedentes.
!
!  ITAB1 est necessaire a la creation de ITAB2 et ITAB2 permet de boucler
!  directement sur les valeurs de JK utiles, et de dimensionner au plus
!  juste les tableaux ZW et ZNORM.
!
  ALLOCATE ( itab1(ikmax) )
  itab1 = 0
!DP#ifdef SCALAR
  ALLOCATE ( itab2((itronc-1)*(imtronc-1)) )
!DP#else
!DP        ALLOCATE ( ITAB2(0:(ITRONC-1)*(IMTRONC-1)) )
!DP        ITAB2 = 0
!DP#endif
  DO jm = 1,imtronc
    imlim=istron-jm
    ideb=max(fa%CADRE(irangc)%NOMPAR(2*jm+3)+4*(1+imlim), &
&            fa%CADRE(irangc)%NOMPAR(2*jm+3)+4)
    ifin=fa%CADRE(irangc)%NOMPAR(2*jm+4)
    DO jind=ideb,ifin
      ioff=jind-fa%CADRE(irangc)%NOMPAR(2*jm+3)
      jn=ioff/4
      jk=(jn**2 + jm**2)
!  On conserve les valeurs de JK atteintes
      itab1(jk) = 1
    ENDDO
  ENDDO
!  On construit le tableau ITAB2 des valeurs de JK atteintes
  ik = 0
  DO jk = 1,ikmax
!DP#ifdef SCALAR
    IF (itab1(jk).GT.0) THEN
      ik = ik+1
      itab2(ik) = jk
      itab1(jk) = ik
    ENDIF
!DP#else
!DP       IK = IK + ITAB1(JK)
!DP       ITAB2(IK) = ITAB2(IK) + JK*ITAB1(JK)
!DP       ITAB1(JK) = IK
!DP#endif
  ENDDO
  inmin = 1
  inmax = ik
!
!  CAS ARPEGE
!
ELSE
  inmin = istron + 1
  inmax = itronc
ENDIF
!
!**
!     2.  -  CALCUL DES NORMES ET DES POIDS
!-----------------------------------------------------------------------
!
ALLOCATE ( zw(0:inmax), znorm(0:inmax) )
znorm = 0.0_jpdblr
!
IF (llmlam) THEN
!
!  CAS ALADIN
!
  DO jm = 1,imtronc
    imlim=istron-jm
    ideb=max(fa%CADRE(irangc)%NOMPAR(2*jm+3)+4*(1+imlim), &
&            fa%CADRE(irangc)%NOMPAR(2*jm+3)+4)
    ifin=fa%CADRE(irangc)%NOMPAR(2*jm+4)
    DO jind=ideb,ifin
      ioff=jind-fa%CADRE(irangc)%NOMPAR(2*jm+3)
      jn=ioff/4
      jk=(jn**2 + jm**2)
      znorm(itab1(jk)) = max(znorm(itab1(jk)), abs(pspec(jind)))
    ENDDO
  ENDDO
  DEALLOCATE ( itab1 )
!
!  CAS ARPEGE
!
ELSE
!
  DO jn=0,inmax
    DO jm=-jn,jn
      im=abs(jm)
      IF (jm < 0) THEN
        jind=fa%CADRE(irangc)%NDIM0GG(im)+(jn-im)*2 +1
      ELSE
        jind=fa%CADRE(irangc)%NDIM0GG(im)+(jn-im)*2
      ENDIF
      znorm(jn) = max(znorm(jn), abs(pspec(jind)))
    ENDDO
  ENDDO
!
ENDIF
!
zz = REAL(inmax-inmin+1,jpdblr)
DO in = inmin, inmax
  zw(in) = zz / REAL(in-inmin+1,jpdblr)
!
!     Ensure the norms have a value which is not too small in case of
!     problems with math functions (e.g. LOG).
!
  IF(znorm(in).LT.zeps) THEN
    znorm(in) = zeps
    zw(in) = 100._jpdblr*zeps
  ENDIF
ENDDO
!*
!     3.  -  AJUSTEMENT PAR LES MOINDRES CARRES POUR DETERMINER LA PENTE
!            (ZBETA1).
!-----------------------------------------------------------------------
!
!
!
zxmw  = 0.0_jpdblr
zymw  = 0.0_jpdblr
zwsum = 0.0_jpdblr
!
!     Sum the weighted X and Ys.
DO jn = inmin, inmax
  IF (llmlam) THEN
    zx =  log(REAL(ITAB2(JN),JPDBLR))
  ELSE
    zx =  log(REAL(JN*(JN+1),JPDBLR))
  ENDIF
  zy    = log(znorm(jn))
  zxmw  = zxmw+zx*zw(jn)
  zymw  = zymw+zy*zw(jn)
  zwsum = zwsum+zw(jn)
ENDDO
!
!     Form mean weighted X and Y.
zxmw  = zxmw / zwsum
zymw  = zymw / zwsum
zsum1 = 0.0_jpdblr
zsum2 = 0.0_jpdblr
!
!     Perform a least square fit for the equation
DO jn = inmin, inmax
  IF (llmlam) THEN
    zx =  log(REAL(ITAB2(JN),JPDBLR))
  ELSE
    zx =  log(REAL(JN*(JN+1),JPDBLR))
  ENDIF
  zy    = log(znorm(jn))
  zsum1 = zsum1+zw(jn)*(zy-zymw)*(zx-zxmw)
  zsum2 = zsum2+zw(jn)*(zx-zxmw)**2
ENDDO
!
!
DEALLOCATE ( zw, znorm )
IF (llmlam) DEALLOCATE ( itab2 )
!
!     On peut en tirer la pente recherchee:
zbeta1 = zsum1 / zsum2
!
!     Et finalement, la puissance de laplacien, arrondie en un
!     nombre entier de milliemes.
zp = -zbeta1
zp = max(-9.999_jpdblr, min(9.999_jpdblr, zp))
kpulap = nint( zp * 1000.0_jpdblr, kind=jplikb )
!**
!    10.  -  PHASE TERMINALE : MESSAGERIE EVENTUELLE,
!            VIA LE SOUS-PROGRAMME "FAIPAR" .
!-----------------------------------------------------------------------
!
1001 CONTINUE
llfata=llmoer(krep,krang)
!
IF (fa%LFAMOP.OR.llfata) THEN
  inimes=2
  clnspr='FAPULA'
  inumer=jpniil
  WRITE (unit=clmess,fmt='(''KREP='',I4,'', KRANG='',I4,  &
&         '', ISTRON='',I4,'', KPULAP='',I6)')            &
&         krep,krang,istron,kpulap
  CALL faipar_fort                                        &
&                 (fa, inumer,inimes,krep,.false.,clmess, &
&                  clnspr,clnspr,.false.)
ENDIF
!
IF (lhook) CALL dr_hook('FAPULA_MT',1,zhook_handle)

CONTAINS

#include "facom2.llmoer.h"

END SUBROUTINE

!INTF KREP            OUT                               
!INTF KRANG         IN                                  
!INTF PSPEC         IN    DIMS=*                        
!INTF KPULAP          OUT