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

Go to the documentation of this file.
00001 C Jan-2011 P. Marguinaud Thread-safe FA
00002       SUBROUTINE FAXION_MT (FA,  PCHAME, KPUISS, KDIMNC, KLCHAM, PMIN, 
00003      S                      PMAX, KNBITS, LDARPE, PECART, LDMLAM, 
00004      S                      KNOZPA, KSTROF, KTRONC, KXLOPA )
00005       USE FA_MOD, ONLY : FA_COM
00006       USE PARKIND1, ONLY : JPRB
00007       USE YOMHOOK , ONLY : LHOOK, DR_HOOK
00008 C****
00009 C        Sous-programme INTERNE du logiciel de Fichiers ARPEGE:
00010 C     calcul de l'erreur totale par le codage GRIB
00011 C     sur un champ en coefficients spectraux, en mode "complex packing".
00012 C**
00013 C    Arguments :
00014 C    ( Tableau ) PCHAME (Entree) ==> Champ en coef. spectraux en entree;
00015 C                KPUISS (Entree) ==> Puissance de laplacien a utiliser;
00016 C                KDIMNC (Entree) ==> Taille de la zone non codee;
00017 C                KLCHAM (Entree) ==> Longueur totale du champ en c.spx.;
00018 C                PMIN   (Entree) ==> Minimum precalcule du champ module;
00019 C                PMAX   (Entree) ==> Maximum      "      "    "     "  ;
00020 C                KNBITS (Entree) ==> Nombre de bits par valeur a coder;
00021 C                LDARPE (Entree) ==> Vrai si GRIB non standard;
00022 C                PECART (Sortie) ==> Erreur relative maximum commise,
00023 C                                    en valeur absolue;
00024 C                LDMLAM (Entree) ==> Vrai si fichier aladin
00025 C    (Tableau)   KNOZPA (Entree) ==> Nombre d'onde zonal maxi par parallele;
00026 C                                    (du pole nord vers l'equateur seulement)
00027 C                KSTROF (Entree) ==> Sous-troncature non compactee EFFECTIVE (coef. spectraux)
00028 C                KTRONC (Entree) ==> Troncature (NS pour Aladin)
00029 C                KXLOPA (Entree) ==> Nombre maximum de longitudes par cercle de latitudes
00030 C
00031 C     Modifications
00032 C     -------------
00033 C
00034 C  Juillet 1998, J. Clochard, SCEM/TTI/DAO:
00035 C
00036 C    -Le critere a minimiser est mis sous la forme de fonctions
00037 C    -Pour le cas LDARPE=.TRUE., utilisation des memes gardes-fous
00038 C     que dans CODEGA.
00039 C
00040 C  Septem. 2000, D. Paradis, SCEM/TTI/DEV:
00041 C
00042 C    -Ajustement de la puissance Dolby pour des nb d'onde (JN)
00043 C     uniquement inferieurs a MIN(FA%MTRONC(),(FA%NXLOPA()-1)/3): on
00044 C     s'affranchit de l'influence des grands nb d'onde sur des
00045 C     grilles lineaires pour le choix du Dolby.
00046 C
00047 C  Mars 2003, R. El Khatib CNRM/GMAP:
00048 C
00049 C    -Optimisation : ISAMAX remplace par MAXLOC
00050 C
00051 C  March 2010: J. Masek - fix of precomputed optimal Laplacian power
00052 #include "precision.h"
00053 C
00054 C
00055       TYPE(FA_COM) :: FA
00056       INTEGER KPUISS, KDIMNC, KLCHAM, KNBITS, KSTROF
00057       INTEGER KTRONC, KXLOPA
00058 C
00059       INTEGER KNOZPA (FA%JPXIND)
00060       INTEGER ILOC(1)
00061 C
00062       REAL (KIND=JPDBLR) PMIN, PMAX, PECART
00063       REAL (KIND=JPDBLR) PCHAME(KLCHAM)
00064       REAL (KIND=JPDBLR) ZERR  (KLCHAM)
00065       REAL (KIND=JPDBLR) ZECART_LOC
00066 C
00067       LOGICAL LDARPE, LDMLAM
00068 C
00069       INTEGER JN, JIND, ISCALE, J, INDICE, IPUISX, IOFF, IM
00070       INTEGER INDLAP, IRAPOR, IPUISR, IMLIM, IEXP, IMANT, IPUIS2
00071       INTEGER IZZZZZ, IDEB, IFIN, ITRDOL, ILCHADO
00072 C
00073       REAL (KIND=JPDBLR) ZREFER, ZDIFFR, ZCOMPA, ZAUXI1, ZAUXI2, ZS
00074       REAL (KIND=JPDBLR) ZYPR, ZYP0, ZCRITR, ZCRIT0, ZCOEFF
00075 #include "facom_mt.h"
00076 C
00077 C
00078 C
00079 C       Fonctions intermediaires utilisees pour alleger l'ecriture.
00080 C
00081       REAL(KIND=JPRB) :: ZHOOK_HANDLE
00082       ZYPR ( IZZZZZ, ZCOEFF ) = ZREFER + ZAUXI2 * ANINT ( 
00083      S                  ZAUXI1 * ( PCHAME(IZZZZZ)/ZCOEFF - ZREFER ) )
00084 C
00085       ZYP0 ( IZZZZZ ) = ZREFER + ZAUXI2 * ANINT (
00086      S                  ZAUXI1 * ( PCHAME(IZZZZZ) - ZREFER ) )
00087 C
00088 C       Fonction "critere" a minimiser lorsque la modulation de la
00089 C       puissance de laplacien est possible.
00090 C
00091       ZCRITR ( IZZZZZ, ZCOEFF ) = ( ZYPR(IZZZZZ,ZCOEFF)*ZCOEFF - 
00092      S                      PCHAME(IZZZZZ) )**2
00093 C
00094 C       Meme chose, mais pour la puissance 0.
00095 C
00096       ZCRIT0 ( IZZZZZ ) = (ZYP0(IZZZZZ) - PCHAME(IZZZZZ) )**2
00097       IF (LHOOK) CALL DR_HOOK('FAXION_MT',0,ZHOOK_HANDLE)
00098 C**
00099 C     0.  -  INITIALISATIONS FAITES AU PREMIER APPEL DU SOUS-PROGRAMME.
00100 C-----------------------------------------------------------------------
00101 C
00102       IF (FA%FAXION_LLPREA) THEN
00103         FA%FAXION_ZEPSIL=TINY(FA%FAXION_ZEPSIL)
00104         FA%FAXION_ISCALX=99
00105         FA%FAXION_LLPREA=.FALSE.
00106       ENDIF
00107 C**
00108 C     1.  -  CALCUL DIRECT, EN ESSAYANT DE MINIMISER LES OPERATIONS.
00109 C            ( a l'aide des puissances precalculees dans FA%XLAP2D )
00110 C-----------------------------------------------------------------------
00111 C
00112       IF (LDARPE) THEN
00113 C
00114         ZDIFFR=PMAX-PMIN
00115 C
00116         IF ( ZDIFFR .LE. FA%FAXION_ZEPSIL ) THEN
00117           ZAUXI1=MIN ( ABS (PMIN), ABS (PMAX) )
00118           IF ( ZAUXI1 .LE. FA%FAXION_ZEPSIL ) ZAUXI1=0.
00119           PMAX=SIGN (ZAUXI1,PMAX)
00120           PMIN=PMAX
00121           ZAUXI1=0.
00122           ZAUXI2=0.
00123         ELSE
00124           ZCOMPA=2**KNBITS-1
00125           ZAUXI1=ZCOMPA/ZDIFFR
00126           ZAUXI2=ZDIFFR/ZCOMPA
00127         ENDIF
00128 C
00129         ZREFER=PMIN
00130 C
00131       ELSE
00132 C
00133         CALL CONFI (PMIN,IEXP,IMANT,ZREFER)
00134         ZS = (PMAX-ZREFER)/(2**(KNBITS+1)-1)
00135         ZAUXI1=1.
00136         ZAUXI2=2.
00137         IF (ZS.NE.0.) ZS = LOG(ZS)/LOG(ZAUXI2) + ZAUXI2
00138         ISCALE = MIN(INT(ZS),INT(ZS+SIGN(ZAUXI1,ZS)))
00139         ISCALE = MAX(-FA%FAXION_ISCALX,MIN(FA%FAXION_ISCALX,ISCALE))
00140         ZAUXI1 = ZAUXI2**(-ISCALE)
00141         ZAUXI2 = ZAUXI2**ISCALE
00142 C
00143       ENDIF
00144 C
00145 C  Calcul de la borne superieure de la plage des nb d'onde
00146 C  servant pour l'ajustement du Dolby (borne inferieure=KSTROF)
00147 C  valable uniquement dans le cas Arpege pour ne pas prendre
00148 C  en compte les coeff spectraux relatifs aux nb d'onde non
00149 C  quadratiques (grille lineaire).
00150 C
00151       IF (LDARPE.AND..NOT.LDMLAM) THEN
00152         ITRDOL = MIN ( KTRONC , (KXLOPA-1)/3 )
00153         ILCHADO = (ITRDOL+1)**2
00154       ENDIF
00155 C
00156       IF (KPUISS.EQ.0) THEN
00157 C
00158         IF (LDMLAM) THEN 
00159 !$OMP PARALLEL DO PRIVATE(JN,JIND) IF(FA%LOPENMP)
00160           DO 1010 JN=0,KTRONC
00161           DO 1010 JIND=KNOZPA(2*JN+3),KNOZPA(2*JN+4)
00162           ZERR(JIND)=ZCRIT0(JIND)
00163  1010     CONTINUE
00164 !$OMP END PARALLEL DO
00165         ELSE
00166           DO 101 J=KDIMNC+1,ILCHADO
00167           ZERR(J)=ZCRIT0(J)
00168   101     CONTINUE
00169         ENDIF
00170 C
00171       ELSE
00172         IPUISX=IABS (KPUISS)
00173 C
00174         IF (KPUISS.GT.0) THEN
00175           INDICE=1
00176         ELSE
00177           INDICE=0
00178         ENDIF
00179 C
00180         IF (IPUISX.LE.FA%JPUILA) THEN
00181 C
00182 C        Dans ce cas precis, on pourrait aussi remplacer la division
00183 C     par ZCOEFF par une multiplication par FA%XLAP2D(J,IPUISX,1-INDICE).
00184 C     Il y aurait gain en calcul, mais reference memoire supplementaire.
00185 C
00186           IF (LDMLAM) THEN
00187 #if !defined(RS6K)
00188 !$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP,ZCOEFF)
00189 !$OMP&  IF(FA%LOPENMP)
00190 #endif
00191             DO 1020 JN=1,KTRONC
00192             DO 1020 JIND=KNOZPA(2*JN+3)+4,KNOZPA(2*JN+4)
00193             IOFF=JIND-KNOZPA(2*JN+3)
00194             IM=IOFF/4 
00195             INDLAP=((JN-1)*FA%JPXTRO)+IM
00196             ZERR(JIND)=ZCRITR(JIND,FA%XLAP2DA(INDLAP,IPUISX,INDICE))
00197  1020       CONTINUE
00198 #if !defined(RS6K)
00199 !$OMP END PARALLEL DO
00200 #endif
00201           ELSE
00202             DO 102 J=KDIMNC+1,ILCHADO
00203             ZERR(J)=ZCRITR(J,FA%XLAP2D(J,IPUISX,INDICE))
00204   102       CONTINUE
00205           ENDIF
00206 C
00207         ELSEIF (IPUISX.LE.2*FA%JPUILA) THEN
00208           IPUIS2=IPUISX/2
00209 C
00210           IF (IPUISX.EQ.2*IPUIS2) THEN
00211 C
00212             IF (LDMLAM) THEN
00213 #if !defined(RS6K)
00214 !$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP,ZCOEFF)
00215 !$OMP&  IF(FA%LOPENMP)
00216 #endif
00217               DO 1030 JN=1,KTRONC
00218               DO 1030 JIND=KNOZPA(2*JN+3)+4,KNOZPA(2*JN+4)
00219               IOFF=JIND-KNOZPA(2*JN+3)
00220               IM=IOFF/4 
00221               INDLAP=((JN-1)*FA%JPXTRO)+IM
00222               ZERR(JIND)=ZCRITR(JIND,
00223      S                   FA%XLAP2DA(INDLAP,IPUIS2,INDICE)**2)
00224  1030         CONTINUE
00225 #if !defined(RS6K)
00226 !$OMP END PARALLEL DO
00227 #endif
00228             ELSE
00229               DO 103 J=KDIMNC+1,ILCHADO
00230               ZERR(J)=ZCRITR(J,FA%XLAP2D(J,IPUIS2,INDICE)**2)
00231   103         CONTINUE
00232             ENDIF
00233 C
00234           ELSE
00235 C
00236           IF (LDMLAM) THEN
00237 #if !defined(RS6K)
00238 !$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP,ZCOEFF)
00239 !$OMP&  IF(FA%LOPENMP)
00240 #endif
00241              DO 1040 JN=1,KTRONC
00242              DO 1040 JIND=KNOZPA(2*JN+3)+4,KNOZPA(2*JN+4)
00243               IOFF=JIND-KNOZPA(2*JN+3)
00244               IM=IOFF/4 
00245               INDLAP=((JN-1)*FA%JPXTRO)+IM
00246               ZERR(JIND)=ZCRITR(JIND,FA%XLAP2DA(INDLAP,FA%JPUILA,INDICE)
00247      S             *FA%XLAP2DA(INDLAP,IPUISX-FA%JPUILA,INDICE))
00248  1040         CONTINUE
00249 #if !defined(RS6K)
00250 !$OMP END PARALLEL DO
00251 #endif
00252           ELSE
00253               DO 104 J=KDIMNC+1,ILCHADO
00254               ZERR(J)=ZCRITR(J,FA%XLAP2D(J,FA%JPUILA,INDICE)
00255      S             *FA%XLAP2D(J,IPUISX-FA%JPUILA,INDICE))
00256   104       CONTINUE
00257           ENDIF
00258 C
00259           ENDIF
00260 C
00261         ELSE
00262           IRAPOR=1+(IPUISX-1)/FA%JPUILA
00263           IPUISR=IPUISX/IRAPOR
00264 C
00265           IF (IPUISX.EQ.IRAPOR*IPUISR) THEN
00266 C
00267             IF (LDMLAM) THEN
00268 #if !defined(RS6K)
00269 !$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP,ZCOEFF)
00270 !$OMP&  IF(FA%LOPENMP)
00271 #endif
00272               DO 1050 JN=1,KTRONC
00273               DO 1050 JIND=KNOZPA(2*JN+3)+4,KNOZPA(2*JN+4)
00274               IOFF=JIND-KNOZPA(2*JN+3)
00275               IM=IOFF/4 
00276               INDLAP=((JN-1)*FA%JPXTRO)+IM
00277               ZERR(JIND)=ZCRITR(JIND,
00278      S               FA%XLAP2DA(INDLAP,IPUISR,INDICE)**IRAPOR)
00279  1050         CONTINUE
00280 #if !defined(RS6K)
00281 !$OMP END PARALLEL DO
00282 #endif
00283             ELSE
00284               DO 105 J=KDIMNC+1,ILCHADO
00285               ZERR(J)=ZCRITR(J,FA%XLAP2D(J,IPUISR,INDICE)**IRAPOR)
00286   105         CONTINUE
00287             ENDIF
00288 C
00289           ELSE
00290 C
00291            IF (LDMLAM) THEN
00292 #if !defined(RS6K)
00293 !$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP,ZCOEFF)
00294 !$OMP&  IF(FA%LOPENMP)
00295 #endif
00296               DO 1060 JN=1,KTRONC
00297               DO 1060 JIND=KNOZPA(2*JN+3)+4,KNOZPA(2*JN+4)
00298               IOFF=JIND-KNOZPA(2*JN+3)
00299               IM=IOFF/4 
00300               INDLAP=((JN-1)*FA%JPXTRO)+IM
00301               ZERR(JIND)=ZCRITR(JIND,
00302      S           FA%XLAP2DA(INDLAP,FA%JPUILA,INDICE)**(IRAPOR-1)
00303      S          *FA%XLAP2DA(INDLAP,IPUISX-FA%JPUILA*(IRAPOR-1),INDICE))
00304  1060         CONTINUE
00305 #if !defined(RS6K)
00306 !$OMP END PARALLEL DO
00307 #endif
00308            ELSE
00309               DO 106 J=KDIMNC+1,ILCHADO
00310               ZERR(J)=ZCRITR(J,
00311      S             FA%XLAP2D(J,FA%JPUILA,INDICE)**(IRAPOR-1)
00312      S            *FA%XLAP2D(J,IPUISX-FA%JPUILA*(IRAPOR-1),INDICE))
00313   106         CONTINUE
00314             ENDIF
00315 C
00316           ENDIF
00317 C
00318         ENDIF
00319 C
00320       ENDIF
00321 C
00322       PECART=0.
00323       IF (LDMLAM) THEN
00324 !$OMP PARALLEL IF(FA%LOPENMP)
00325 !$OMP& PRIVATE(JN,IMLIM,IDEB,IFIN,JIND,ZECART_LOC)
00326 !$OMP& REDUCTION(+:PECART)
00327 !$OMP DO
00328         DO 1070 JN=1,KTRONC
00329         IMLIM=KSTROF-JN
00330         IDEB=MAX(KNOZPA(2*JN+3)+4,KNOZPA(2*JN+3)+4*(1+IMLIM))
00331         IFIN=KNOZPA(2*JN+4)
00332         ZECART_LOC=0.
00333         DO 1071 JIND=IDEB,IFIN
00334           ZECART_LOC=ZECART_LOC+ZERR(JIND)
00335  1071   CONTINUE
00336         PECART=PECART+ZECART_LOC
00337 C
00338  1070   CONTINUE
00339 !$OMP END DO 
00340 !$OMP END PARALLEL
00341       ELSE
00342         DO J=KDIMNC+1,ILCHADO
00343           PECART=PECART+ZERR(J)
00344         ENDDO
00345       ENDIF
00346 C
00347       IF (FA%LFAMOP) THEN
00348         IF (LDARPE.AND..NOT.LDMLAM) THEN
00349           WRITE (UNIT=FA%NULOUT,FMT=*)
00350      S           'FAXION: KPUISS=', KPUISS,', KDIMNC=',KDIMNC,
00351      S           ', KLCHAM=', KLCHAM, ', ITRDOL=', ITRDOL, ', ILCHADO=',
00352      S           ILCHADO, ', PMIN=', PMIN, ', PMAX=', PMAX
00353         ELSE
00354           WRITE (UNIT=FA%NULOUT,FMT=*)
00355      S           'FAXION: KPUISS=', KPUISS,', KDIMNC=',KDIMNC,
00356      S           ', KLCHAM=', KLCHAM,
00357      S           ', PMIN=', PMIN, ', PMAX=', PMAX
00358         ENDIF
00359         WRITE (UNIT=FA%NULOUT,FMT=*)'FAXION: PECART=', PECART
00360       ENDIF
00361 C
00362       IF (LHOOK) CALL DR_HOOK('FAXION_MT',1,ZHOOK_HANDLE)
00363       END
00364 
00365