faxion.F90

Go to the documentation of this file.

! Oct-2012 P. Marguinaud 64b LFI
! Jan-2011 P. Marguinaud Thread-safe FA
SUBROUTINE faxion_fort                                           &
&                     (fa,  pchame, kpuiss, kdimnc, klcham, pmin,  &
&                      pmax, knbits, ldarpe, pecart, ldmlam,       &
&                      knozpa, kstrof, ktronc, kxlopa )
USE fa_mod, ONLY : fa_com
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 l'erreur totale par le codage GRIB
!     sur un champ en coefficients spectraux, en mode "complex packing".
!**
!    Arguments :
!    ( Tableau ) PCHAME (Entree) ==> Champ en coef. spectraux en entree;
!                KPUISS (Entree) ==> Puissance de laplacien a utiliser;
!                KDIMNC (Entree) ==> Taille de la zone non codee;
!                KLCHAM (Entree) ==> Longueur totale du champ en c.spx.;
!                PMIN   (Entree) ==> Minimum precalcule du champ module;
!                PMAX   (Entree) ==> Maximum      "      "    "     "  ;
!                KNBITS (Entree) ==> Nombre de bits par valeur a coder;
!                LDARPE (Entree) ==> Vrai si GRIB non standard;
!                PECART (Sortie) ==> Erreur relative maximum commise,
!                                    en valeur absolue;
!                LDMLAM (Entree) ==> Vrai si fichier aladin
!    (Tableau)   KNOZPA (Entree) ==> Nombre d'onde zonal maxi par parallele;
!                                    (du pole nord vers l'equateur seulement)
!                KSTROF (Entree) ==> Sous-troncature non compactee EFFECTIVE (coef. spectraux)
!                KTRONC (Entree) ==> Troncature (NS pour Aladin)
!                KXLOPA (Entree) ==> Nombre maximum de longitudes par cercle de latitudes
!
!     Modifications
!     -------------
!
!  Juillet 1998, J. Clochard, SCEM/TTI/DAO:
!
!    -Le critere a minimiser est mis sous la forme de fonctions
!    -Pour le cas LDARPE=.TRUE., utilisation des memes gardes-fous
!     que dans CODEGA.
!
!  Septem. 2000, D. Paradis, SCEM/TTI/DEV:
!
!    -Ajustement de la puissance Dolby pour des nb d'onde (JN)
!     uniquement inferieurs a MIN(FA%MTRONC(),(FA%NXLOPA()-1)/3): on
!     s'affranchit de l'influence des grands nb d'onde sur des
!     grilles lineaires pour le choix du Dolby.
!
!  Mars 2003, R. El Khatib CNRM/GMAP:
!
!    -Optimisation : ISAMAX remplace par MAXLOC
!
!  March 2010: J. Masek - fix of precomputed optimal Laplacian power
!
!
TYPE(fa_com) :: FA
INTEGER (KIND=JPLIKB) KPUISS, KDIMNC, KLCHAM, KNBITS, KSTROF
INTEGER (KIND=JPLIKB) KTRONC, KXLOPA
!
INTEGER (KIND=JPLIKB) KNOZPA (fa%jpxind)
!
REAL (KIND=JPDBLR) PMIN, PMAX, PECART
REAL (KIND=JPDBLR) PCHAME(klcham)
REAL (KIND=JPDBLR) ZERR  (klcham)
REAL (KIND=JPDBLR) ZECART_LOC
!
LOGICAL LDARPE, LDMLAM
!
INTEGER (KIND=JPLIKB) JN, JIND, ISCALE, J
INTEGER (KIND=JPLIKB) INDICE, IPUISX, IOFF, IM
INTEGER (KIND=JPLIKB) INDLAP, IRAPOR, IPUISR
INTEGER (KIND=JPLIKB) IMLIM, IEXP, IMANT, IPUIS2
INTEGER (KIND=JPLIKB) IDEB, IFIN, ITRDOL, ILCHADO
INTEGER (KIND=JPLIKB) IEXP32, IMANT32
!
REAL (KIND=JPDBLR) ZREFER, ZDIFFR, ZCOMPA, ZAUXI1, ZAUXI2, ZS

REAL(KIND=JPRB) :: ZHOOK_HANDLE

IF (lhook) CALL dr_hook('FAXION_MT',0,zhook_handle)
!**
!     0.  -  INITIALISATIONS FAITES AU PREMIER APPEL DU SOUS-PROGRAMME.
!-----------------------------------------------------------------------
!
IF (fa%FAXION_LLPREA) THEN
  fa%FAXION_ZEPSIL=tiny(fa%FAXION_ZEPSIL)
  fa%FAXION_ISCALX=99
  fa%FAXION_LLPREA=.false.
ENDIF
!**
!     1.  -  CALCUL DIRECT, EN ESSAYANT DE MINIMISER LES OPERATIONS.
!            ( a l'aide des puissances precalculees dans FA%XLAP2D )
!-----------------------------------------------------------------------
!
IF (ldarpe) THEN
!
  zdiffr=pmax-pmin
!
  IF ( zdiffr .LE. fa%FAXION_ZEPSIL ) THEN
    zauxi1=min( abs(pmin), abs(pmax) )
    IF ( zauxi1 .LE. fa%FAXION_ZEPSIL ) zauxi1=0._jpdblr
    pmax=sign(zauxi1,pmax)
    pmin=pmax
    zauxi1=0._jpdblr
    zauxi2=0._jpdblr
  ELSE
    zcompa=real(2**knbits-1, jpdblr)
    zauxi1=zcompa/zdiffr
    zauxi2=zdiffr/zcompa
  ENDIF
!
  zrefer=pmin
!
ELSE
!
  CALL confi (pmin,iexp32,imant32,zrefer)
  iexp=iexp32
  imant=imant32

  zs = (pmax-zrefer)/(2**(knbits+1)-1)
  zauxi1=1._jpdblr
  zauxi2=2._jpdblr
  IF (zs.NE.0._jpdblr) zs = log(zs)/log(zauxi2) + zauxi2
  iscale = min(int(zs,jplikb),                 &
&               int(zs+sign(zauxi1,zs),jplikb))
  iscale = max(-fa%FAXION_ISCALX,min(fa%FAXION_ISCALX,iscale))
  zauxi1 = zauxi2**(-iscale)
  zauxi2 = zauxi2**iscale
!
ENDIF
!
!  Calcul de la borne superieure de la plage des nb d'onde
!  servant pour l'ajustement du Dolby (borne inferieure=KSTROF)
!  valable uniquement dans le cas Arpege pour ne pas prendre
!  en compte les coeff spectraux relatifs aux nb d'onde non
!  quadratiques (grille lineaire).
!
IF (ldarpe.AND..NOT.ldmlam) THEN
  itrdol = min( ktronc , (kxlopa-1)/3 )
  ilchado = (itrdol+1)**2
ENDIF
!
IF (kpuiss.EQ.0) THEN
!
  IF (ldmlam) THEN 
!$OMP PARALLEL DO PRIVATE(JN,JIND) IF(FA%LOPENMP)
    DO jn=0,ktronc
    DO jind=knozpa(2*jn+3),knozpa(2*jn+4)
    zerr(jind)=zcrit0(jind)
    ENDDO
    ENDDO
!$OMP END PARALLEL DO
  ELSE
    DO j=kdimnc+1,ilchado
    zerr(j)=zcrit0(j)
    ENDDO
  ENDIF
!
ELSE
  ipuisx=abs(kpuiss)
!
  IF (kpuiss.GT.0) THEN
    indice=1
  ELSE
    indice=0
  ENDIF
!
  IF (ipuisx.LE.fa%JPUILA) THEN
!
!        Dans ce cas precis, on pourrait aussi remplacer la division
!     par ZCOEFF par une multiplication par FA%XLAP2D(J,IPUISX,1-INDICE).
!     Il y aurait gain en calcul, mais reference memoire supplementaire.
!
    IF (ldmlam) THEN
#ifndef RS6K
!$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP) &
!$OMP&  IF(FA%LOPENMP)
#endif
      DO jn=1,ktronc
      DO jind=knozpa(2*jn+3)+4,knozpa(2*jn+4)
      ioff=jind-knozpa(2*jn+3)
      im=ioff/4 
      indlap=((jn-1)*fa%JPXTRO)+im
      zerr(jind)=zcritr(jind,fa%XLAP2DA(indlap,ipuisx,indice))
      ENDDO
      ENDDO
#ifndef RS6K
!$OMP END PARALLEL DO
#endif
    ELSE
      DO j=kdimnc+1,ilchado
      zerr(j)=zcritr(j,fa%XLAP2D(j,ipuisx,indice))
      ENDDO
    ENDIF
!
  ELSEIF (ipuisx.LE.2*fa%JPUILA) THEN
    ipuis2=ipuisx/2
!
    IF (ipuisx.EQ.2*ipuis2) THEN
!
      IF (ldmlam) THEN
#ifndef RS6K
!$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP) &
!$OMP&  IF(FA%LOPENMP)
#endif
        DO jn=1,ktronc
        DO jind=knozpa(2*jn+3)+4,knozpa(2*jn+4)
        ioff=jind-knozpa(2*jn+3)
        im=ioff/4 
        indlap=((jn-1)*fa%JPXTRO)+im
        zerr(jind)=zcritr(jind,                         &
&                   fa%XLAP2DA(indlap,ipuis2,indice)**2)
        ENDDO
        ENDDO
#ifndef RS6K
!$OMP END PARALLEL DO
#endif
      ELSE
        DO j=kdimnc+1,ilchado
        zerr(j)=zcritr(j,fa%XLAP2D(j,ipuis2,indice)**2)
        ENDDO
      ENDIF
!
    ELSE
!
    IF (ldmlam) THEN
#ifndef RS6K
!$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP) &
!$OMP&  IF(FA%LOPENMP)
#endif
       DO jn=1,ktronc
       DO jind=knozpa(2*jn+3)+4,knozpa(2*jn+4)
        ioff=jind-knozpa(2*jn+3)
        im=ioff/4 
        indlap=((jn-1)*fa%JPXTRO)+im
        zerr(jind)=zcritr(jind,fa%XLAP2DA(indlap,fa%JPUILA,indice) &
&             *fa%XLAP2DA(indlap,ipuisx-fa%JPUILA,indice))
       ENDDO
       ENDDO
#ifndef RS6K
!$OMP END PARALLEL DO
#endif
    ELSE
        DO j=kdimnc+1,ilchado
        zerr(j)=zcritr(j,fa%XLAP2D(j,fa%JPUILA,indice) &
&             *fa%XLAP2D(j,ipuisx-fa%JPUILA,indice))
        ENDDO
    ENDIF
!
    ENDIF
!
  ELSE
    irapor=1+(ipuisx-1)/fa%JPUILA
    ipuisr=ipuisx/irapor
!
    IF (ipuisx.EQ.irapor*ipuisr) THEN
!
      IF (ldmlam) THEN
#ifndef RS6K
!$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP) &
!$OMP&  IF(FA%LOPENMP)
#endif
        DO jn=1,ktronc
        DO jind=knozpa(2*jn+3)+4,knozpa(2*jn+4)
        ioff=jind-knozpa(2*jn+3)
        im=ioff/4 
        indlap=((jn-1)*fa%JPXTRO)+im
        zerr(jind)=zcritr(jind,                          &
&               fa%XLAP2DA(indlap,ipuisr,indice)**irapor)
        ENDDO
        ENDDO
#ifndef RS6K
!$OMP END PARALLEL DO
#endif
      ELSE
        DO j=kdimnc+1,ilchado
        zerr(j)=zcritr(j,fa%XLAP2D(j,ipuisr,indice)**irapor)
        ENDDO
      ENDIF
!
    ELSE
!
     IF (ldmlam) THEN
#ifndef RS6K
!$OMP PARALLEL DO PRIVATE(JN,JIND,IOFF,IM,INDLAP) &
!$OMP&  IF(FA%LOPENMP)
#endif
        DO jn=1,ktronc
        DO jind=knozpa(2*jn+3)+4,knozpa(2*jn+4)
        ioff=jind-knozpa(2*jn+3)
        im=ioff/4 
        indlap=((jn-1)*fa%JPXTRO)+im
        zerr(jind)=zcritr(jind,                                   &
&           fa%XLAP2DA(indlap,fa%JPUILA,indice)**(irapor-1)        &
&          *fa%XLAP2DA(indlap,ipuisx-fa%JPUILA*(irapor-1),indice))
        ENDDO
        ENDDO
#ifndef RS6K
!$OMP END PARALLEL DO
#endif
     ELSE
        DO j=kdimnc+1,ilchado
        zerr(j)=zcritr(j,                                      &
&             fa%XLAP2D(j,fa%JPUILA,indice)**(irapor-1)        &
&            *fa%XLAP2D(j,ipuisx-fa%JPUILA*(irapor-1),indice))
        ENDDO
      ENDIF
!
    ENDIF
!
  ENDIF
!
ENDIF
!
pecart=0._jpdblr
IF (ldmlam) THEN
!$OMP PARALLEL IF(FA%LOPENMP)                         &
!$OMP& PRIVATE(JN,IMLIM,IDEB,IFIN,JIND,ZECART_LOC)    &
!$OMP& REDUCTION(+:PECART)
!$OMP DO
  DO jn=1,ktronc
  imlim=kstrof-jn
  ideb=max(knozpa(2*jn+3)+4,knozpa(2*jn+3)+4*(1+imlim))
  ifin=knozpa(2*jn+4)
  zecart_loc=0._jpdblr
  DO jind=ideb,ifin
    zecart_loc=zecart_loc+zerr(jind)
  ENDDO
  pecart=pecart+zecart_loc
!
  ENDDO
!$OMP END DO 
!$OMP END PARALLEL
ELSE
  DO j=kdimnc+1,ilchado
    pecart=pecart+zerr(j)
  ENDDO
ENDIF
!
IF (fa%LFAMOP) THEN
  IF (ldarpe.AND..NOT.ldmlam) THEN
    WRITE (unit=fa%NULOUT,fmt=*)                                    &
&           'FAXION: KPUISS=', kpuiss,', KDIMNC=',kdimnc,           &
&           ', KLCHAM=', klcham, ', ITRDOL=', itrdol, ', ILCHADO=', &
&           ilchado, ', PMIN=', pmin, ', PMAX=', pmax
  ELSE
    WRITE (unit=fa%NULOUT,fmt=*)                          &
&           'FAXION: KPUISS=', kpuiss,', KDIMNC=',kdimnc, &
&           ', KLCHAM=', klcham,                          &
&           ', PMIN=', pmin, ', PMAX=', pmax
  ENDIF
  WRITE (unit=fa%NULOUT,fmt=*)'FAXION: PECART=', pecart
ENDIF
!
IF (lhook) CALL dr_hook('FAXION_MT',1,zhook_handle)

CONTAINS

REAL (KIND=JPDBLR) FUNCTION zypr (IZZZZZ, ZCOEFF)
INTEGER (KIND=JPLIKB) :: IZZZZZ
REAL (KIND=JPDBLR)    :: ZCOEFF
zypr = zrefer + zauxi2 * anint(zauxi1 * (pchame(izzzzz)/zcoeff - zrefer))
END FUNCTION

REAL (KIND=JPDBLR) FUNCTION zyp0 (IZZZZZ)
INTEGER (KIND=JPLIKB) :: IZZZZZ
zyp0 = zrefer + zauxi2 * anint(zauxi1 * (pchame(izzzzz) - zrefer))
END FUNCTION

!
!       Fonction "critere" a minimiser lorsque la modulation de la
!       puissance de laplacien est possible.
!
REAL (KIND=JPDBLR) FUNCTION zcritr (IZZZZZ, ZCOEFF)
INTEGER (KIND=JPLIKB) :: IZZZZZ
REAL (KIND=JPDBLR)    :: ZCOEFF
zcritr = (zypr(izzzzz,zcoeff)*zcoeff - pchame(izzzzz))**2
END FUNCTION

!
!       Meme chose, mais pour la puissance 0.
!
REAL (KIND=JPDBLR) FUNCTION zcrit0 (IZZZZZ)
INTEGER (KIND=JPLIKB) :: IZZZZZ
zcrit0 = (zyp0(izzzzz) - pchame(izzzzz))**2
END FUNCTION

END SUBROUTINE

!INTF PCHAME        IN    DIMS=KLCHAM                   
!INTF KPUISS        IN                                  
!INTF KDIMNC        IN                                  
!INTF KLCHAM        IN                                  
!INTF PMIN          IN                                  
!INTF PMAX          IN                                  
!INTF KNBITS        IN                                  
!INTF LDARPE        IN                                  
!INTF PECART          OUT                               
!INTF LDMLAM        IN                                  
!INTF KNOZPA        IN    DIMS=FA%JPXIND                
!INTF KSTROF        IN                                  
!INTF KTRONC        IN                                  
!INTF KXLOPA        IN