SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/LIB/XRD38/FA/mt/fapula

Go to the documentation of this file.
00001 C Jan-2011 P. Marguinaud Thread-safe FA
00002       SUBROUTINE FAPULA_MT (FA,  KREP, KRANG, PSPEC, KPULAP )
00003       USE FA_MOD, ONLY : FA_COM
00004       USE PARKIND1, ONLY : JPRB
00005       USE YOMHOOK , ONLY : LHOOK, DR_HOOK
00006 C****
00007 C        Sous-programme INTERNE du logiciel de Fichiers ARPEGE:
00008 C     calcul de la PUissance de LAplacien "optimale" pour aplatir
00009 C     le spectre hors de la sous-troncature.
00010 C
00011 C**
00012 C    Arguments : KREP   (Sortie) ==> Code-reponse du sous-programme;
00013 C                KRANG  (Entree) ==> Rang de l'unite logique;
00014 C    ( Tableau ) PSPEC  (Entree) ==> Champ en coef. spectraux en entree;
00015 C                KPULAP (Sortie) ==> Puissance de laplacien calculee.
00016 C*
00017 C
00018 C     Creation
00019 C     --------
00020 C  Avril 1999, D. Paradis, SCEM/TTI/DEV:
00021 C     Ce ss-programme est derive de la routine CALCOP du logiciel GRIBEX
00022 C     du CEPMMT. Il retourne la puissance de laplacien (exprimee en un
00023 C     nombre entier de milliemes) qui aplatit le mieux le spectre des
00024 C     coefficients spectraux (hors sous-troncature).
00025 C     Cette puissance de laplacien, P, sera l'exposant du facteur
00026 C     .(n*(n+1))**P   pour les champs ARPEGE ou n est le nombre d'onde total
00027 C     .(n**2+m**2)**P pour les champs ALADIN ou n est le nombre d'onde meridien
00028 C                                            et m est le nombre d'onde zonal
00029 C     Ce facteur sera applique aux coefficients spectraux hors sous-troncature
00030 C     pour en reduire l'amplitude avant le codage GRIB.
00031 C
00032 C
00033 C     Method.
00034 C     -------
00035 C
00036 C     For ARPEGE case, consider F(n,m) = (n*(n+1))**(-P) * G(n,m),
00037 C     where F(n,m) is the original spectral field and n the total wavenumber.
00038 C     For ALADIN case (LAM), consider F(n,m) = (n**2+m**2)**(-P) * G(n,m),
00039 C     where F(n,m) is the original spectral field, n the meridian
00040 C     wavenumber and m the zonal wavenumber.
00041 C     The aim is to minimize G in a 1 norm with respect to P.  This can only
00042 C     partially be achieved.
00043 C
00044 C     1) For ARPEGE case, what we do is to compute H(n), where
00045 C         H(n) = max(F(n,m))    with respect to m.
00046 C
00047 C        We then perform a least square fit for the equation
00048 C            log(H(n)) = beta0+beta1*log(n*(n+1)).
00049 C
00050 C        To ensure a better fit for the lower end of the spectrum, we
00051 C        apply an (arbitrary) weighting function before fitting
00052 C             W(n) = 1.0 / (n-ISTRON+1)
00053 C
00054 C     2) For ALADIN case, what we do is to compute H(in), where
00055 C         H(in) = max(F(n,m))      where n and m verify: k = n**2 + m**2
00056 C           and "in" is the rank of k among all the reached values by k:
00057 C           we are only interested in k values which verify
00058 C           k = n**2 + m**2, with (n,m) belonging to spectrum without
00059 C           the non-compacted sub-truncation area. We call inmax the number
00060 C           of these k values. "in" belongs to interval 1,inmax.
00061 C
00062 C        We then perform a least square fit for the equation
00063 C            log(H(in)) = beta0+beta1*log(k).
00064 C
00065 C        To ensure a better fit for the lower end of the spectrum, we
00066 C        apply an (arbitrary) weighting function before fitting
00067 C           W(in) = 1.0 / in
00068 C
00069 C     Reference.
00070 C     ----------
00071 C
00072 C     Seber, G.A.F. (1979). Linear Regression Analyses.
00073 C     John Wiley and Sons
00074 C
00075 C     ECMWF Research Department documentation of the IFS
00076 C
00077 C
00078 C     Modifications
00079 C     -------------
00080 C
00081 #include "precision.h"
00082 C
00083 C
00084 C     Subroutine arguments
00085 C
00086       TYPE(FA_COM) :: FA
00087       INTEGER KREP, KRANG, KPULAP
00088 C
00089       REAL (KIND=JPDBLR) PSPEC(*)
00090 C
00091 C     Local arguments
00092 C
00093       INTEGER IRANGC, ITRONC, IMTRONC, INMAX, INMIN
00094       INTEGER IKMAX, IDEB, IFIN, JN, JM, JK, IK, IN
00095       INTEGER IMLIM, JIND, IOFF, INUMER, INIMES, ISTRON
00096       INTEGER, DIMENSION(:), ALLOCATABLE    :: ITAB1, ITAB2
00097 C
00098       REAL (KIND=JPDBLR) ZEPS, ZZ, ZXMW, ZYMW, ZWSUM, ZX, ZY
00099       REAL (KIND=JPDBLR) ZK, ZP, ZBETA1, ZSUM1, ZSUM2
00100       REAL (KIND=JPDBLR), PARAMETER :: NINE = 9.999
00101       REAL (KIND=JPDBLR), DIMENSION(:), ALLOCATABLE    :: ZW, ZNORM
00102 C
00103       LOGICAL LLMLAM
00104 C
00105 #include "facom2.h"
00106 #include "facom_mt.h"
00107 C**
00108 C     1.  -  CONTROLES DES PARAMETRES D'APPEL, INITIALISATIONS.
00109 C-----------------------------------------------------------------------
00110 C
00111       REAL(KIND=JPRB) :: ZHOOK_HANDLE
00112       IF (LHOOK) CALL DR_HOOK('FAPULA_MT',0,ZHOOK_HANDLE)
00113       IF (KRANG.LE.0.OR.KRANG.GT.FA%JPNXFA) THEN
00114         KREP=-66
00115         GOTO 1001
00116       ENDIF
00117 C
00118       ISTRON=FA%NSTROF(KRANG)
00119       IRANGC=FA%NUCADR(KRANG)
00120       ITRONC=FA%MTRONC(IRANGC)
00121       LLMLAM=FA%LIMLAM(IRANGC)
00122 C
00123       IF (LLMLAM) IMTRONC=FA%NOZPAR(2,IRANGC)
00124       IF (ITRONC.LE.ISTRON) THEN
00125         KREP=-88
00126         GOTO 1001
00127       ELSEIF (LLMLAM.AND.IMTRONC.LE.ISTRON) THEN
00128         KREP=-88
00129         GOTO 1001
00130       ELSEIF (LLMLAM.AND.(IMTRONC.GT.3*ITRONC.OR.
00131      S    ITRONC.GT.3*IMTRONC)) THEN
00132 C Il s'agit d'un garde-fou, modifiable (ne pas oublier FARCIS et FACSIM)
00133         KREP=-114
00134         GOTO 1001
00135       ELSE       
00136         KREP=0
00137       ENDIF
00138 C
00139       ZEPS  = 1.0D-15
00140 C
00141       IF (LLMLAM) THEN
00142         IKMAX = ITRONC*ITRONC + IMTRONC*IMTRONC
00143 C
00144 C  Pour le cas ALADIN, il faut creer les tableaux ITAB1 et ITAB2 :
00145 C  ITAB1(JK) = 0 sauf pour JK tel qu'il existe (n,m) appartenant au spectre
00146 C    hors de la sous-troncature non-compactee verifiant JK = n**2 + m**2:
00147 C    on a alors ITAB1(JK) = 1
00148 C  ITAB2 regroupe toutes les valeurs de JK precedentes.
00149 C
00150 C  ITAB1 est necessaire a la creation de ITAB2 et ITAB2 permet de boucler
00151 C  directement sur les valeurs de JK utiles, et de dimensionner au plus
00152 C  juste les tableaux ZW et ZNORM.
00153 C
00154         ALLOCATE ( ITAB1(IKMAX) )
00155         ITAB1 = 0
00156 CDP#ifdef SCALAR
00157         ALLOCATE ( ITAB2((ITRONC-1)*(IMTRONC-1)) )
00158 CDP#else
00159 CDP        ALLOCATE ( ITAB2(0:(ITRONC-1)*(IMTRONC-1)) )
00160 CDP        ITAB2 = 0
00161 CDP#endif
00162         DO JM = 1,IMTRONC
00163           IMLIM=ISTRON-JM
00164           IDEB=MAX(FA%NOMPAR(2*JM+3,IRANGC)+4*(1+IMLIM),
00165      S            FA%NOMPAR(2*JM+3,IRANGC)+4)
00166           IFIN=FA%NOMPAR(2*JM+4,IRANGC)
00167           DO JIND=IDEB,IFIN
00168             IOFF=JIND-FA%NOMPAR(2*JM+3,IRANGC)
00169             JN=IOFF/4
00170             JK=(JN**2 + JM**2)
00171 C  On conserve les valeurs de JK atteintes
00172             ITAB1(JK) = 1
00173           ENDDO
00174         ENDDO
00175 C  On construit le tableau ITAB2 des valeurs de JK atteintes
00176         IK = 0
00177         DO JK = 1,IKMAX
00178 CDP#ifdef SCALAR
00179           IF (ITAB1(JK).GT.0) THEN
00180             IK = IK+1
00181             ITAB2(IK) = JK
00182             ITAB1(JK) = IK
00183           ENDIF
00184 CDP#else
00185 CDP       IK = IK + ITAB1(JK)
00186 CDP       ITAB2(IK) = ITAB2(IK) + JK*ITAB1(JK)
00187 CDP       ITAB1(JK) = IK
00188 CDP#endif
00189         ENDDO
00190         INMIN = 1
00191         INMAX = IK
00192 C
00193 C  CAS ARPEGE
00194 C
00195       ELSE
00196         INMIN = ISTRON + 1
00197         INMAX = ITRONC
00198       ENDIF
00199 C
00200 C**
00201 C     2.  -  CALCUL DES NORMES ET DES POIDS
00202 C-----------------------------------------------------------------------
00203 C
00204       ALLOCATE ( ZW(INMAX), ZNORM(INMAX) )
00205       ZNORM = 0.0D0
00206 C
00207       IF (LLMLAM) THEN
00208 C
00209 C  CAS ALADIN
00210 C
00211         DO JM = 1,IMTRONC
00212           IMLIM=ISTRON-JM
00213           IDEB=MAX(FA%NOMPAR(2*JM+3,IRANGC)+4*(1+IMLIM),
00214      S            FA%NOMPAR(2*JM+3,IRANGC)+4)
00215           IFIN=FA%NOMPAR(2*JM+4,IRANGC)
00216           DO JIND=IDEB,IFIN
00217             IOFF=JIND-FA%NOMPAR(2*JM+3,IRANGC)
00218             JN=IOFF/4
00219             JK=(JN**2 + JM**2)
00220             ZNORM(ITAB1(JK)) = MAX(ZNORM(ITAB1(JK)), ABS(PSPEC(JIND)))
00221           ENDDO
00222         ENDDO
00223         DEALLOCATE ( ITAB1 )
00224 C
00225 C  CAS ARPEGE
00226 C
00227       ELSE
00228 C
00229         DO JN = INMIN , INMAX
00230           IDEB = JN**2+1
00231           IFIN = (1+JN)**2
00232           DO JIND = IDEB , IFIN
00233             ZNORM(JN) = MAX(ZNORM(JN), ABS(PSPEC(JIND)))
00234           ENDDO
00235         ENDDO
00236 C
00237       ENDIF
00238 C
00239       ZZ = REAL(INMAX-INMIN+1,JPDBLR)
00240       DO IN = INMIN, INMAX
00241         ZW(IN) = ZZ / REAL(IN-INMIN+1,JPDBLR)
00242 C
00243 C     Ensure the norms have a value which is not too small in case of
00244 C     problems with math functions (e.g. LOG).
00245 C
00246         IF(ZNORM(IN).LT.ZEPS) THEN
00247           ZNORM(IN) = ZEPS
00248           ZW(IN) = 100.D0*ZEPS
00249         ENDIF
00250       ENDDO
00251 C*
00252 C     3.  -  AJUSTEMENT PAR LES MOINDRES CARRES POUR DETERMINER LA PENTE
00253 C            (ZBETA1).
00254 C-----------------------------------------------------------------------
00255 C
00256 C
00257  300  CONTINUE
00258 C
00259       ZXMW  = 0.0D0
00260       ZYMW  = 0.0D0
00261       ZWSUM = 0.0D0
00262 C
00263 C     Sum the weighted X and Ys.
00264       DO JN = INMIN, INMAX
00265         IF (LLMLAM) THEN
00266           ZX =  LOG(REAL(ITAB2(JN),JPDBLR))
00267         ELSE
00268           ZX =  LOG(REAL(JN*(JN+1),JPDBLR))
00269         ENDIF
00270         ZY    = LOG(ZNORM(JN))
00271         ZXMW  = ZXMW+ZX*ZW(JN)
00272         ZYMW  = ZYMW+ZY*ZW(JN)
00273         ZWSUM = ZWSUM+ZW(JN)
00274       ENDDO
00275 C
00276 C     Form mean weighted X and Y.
00277       ZXMW  = ZXMW / ZWSUM
00278       ZYMW  = ZYMW / ZWSUM
00279       ZSUM1 = 0.0D0
00280       ZSUM2 = 0.0D0
00281 C
00282 C     Perform a least square fit for the equation
00283       DO JN = INMIN, INMAX
00284         IF (LLMLAM) THEN
00285           ZX =  LOG(REAL(ITAB2(JN),JPDBLR))
00286         ELSE
00287           ZX =  LOG(REAL(JN*(JN+1),JPDBLR))
00288         ENDIF
00289         ZY    = LOG(ZNORM(JN))
00290         ZSUM1 = ZSUM1+ZW(JN)*(ZY-ZYMW)*(ZX-ZXMW)
00291         ZSUM2 = ZSUM2+ZW(JN)*(ZX-ZXMW)**2
00292       ENDDO
00293 C
00294 C
00295       DEALLOCATE ( ZW, ZNORM )
00296       IF (LLMLAM) DEALLOCATE ( ITAB2 )
00297 C
00298 C     On peut en tirer la pente recherchee:
00299       ZBETA1 = ZSUM1 / ZSUM2
00300 C
00301 C     Et finalement, la puissance de laplacien, arrondie en un
00302 C     nombre entier de milliemes.
00303       ZP = -ZBETA1
00304       ZP = MAX(-NINE, MIN(NINE, ZP))
00305       KPULAP = NINT( ZP * 1000.0 )
00306 C**
00307 C    10.  -  PHASE TERMINALE : MESSAGERIE EVENTUELLE,
00308 C            VIA LE SOUS-PROGRAMME "FAIPAR" .
00309 C-----------------------------------------------------------------------
00310 C
00311  1001 CONTINUE
00312       LLFATA=LLMOER (KREP,KRANG)
00313 C
00314       IF (FA%LFAMOP.OR.LLFATA) THEN
00315         INIMES=2
00316         CLNSPR='FAPULA'
00317         INUMER=FA%JPNIIL
00318         WRITE (UNIT=CLMESS,FMT='(''KREP='',I4,'', KRANG='
00319 ',I4,     S         '', ISTRON='',I4,'', KPULAP='',I6)')
00320      S     KREP,KRANG,ISTRON,KPULAP
00321         CALL FAIPAR_MT (FA, INUMER,INIMES,KREP,.FALSE.,CLMESS,
00322      S               CLNSPR,CLNSPR,.FALSE.)
00323       ENDIF
00324 C
00325       IF (LHOOK) CALL DR_HOOK('FAPULA_MT',1,ZHOOK_HANDLE)
00326       END
00327