Ignore:
Timestamp:
Mar 2, 2015, 5:11:07 PM (9 years ago)
Author:
lguez
Message:

Bug fix in fxhyp. There was some bad logic for the shifting of
longitude grids rlonuv, rlonv, rlonm025 and rlonp025 toward [- pi,
pi]. In some cases, one of the four grids was not shifted and the
others were. For example, you could see the bug with: iim = 48, clon =
20, grossismx = 3, dzoomx = 0.15. The bad logic involved the variable
is2 in the loop on ik = 1, 4. The loop included some tests depending
on ik. The whole thing was quite convoluted. Rewrote fxhyp. In
particular, replaced the loop on ik by a call to a new procedure,
invert_zoom_x. fxhyp.F was in dyn3d, it is replaced by fxhyp_m.F90
which is in dyn3d_common. Removed several arguments of fxhyp: zoom
parameters are accessed through serre.h; rlonm025 and rlonp025 were
not used in inigeom; min and max of longitude steps are written inside
fxhyp.

Some simplifications and modernizations in fyhyp. In particular,
several arguments are removed: zoom parameters are accessed through
serre.h; yprimv was not used in inigeom; min and max of latitude steps
are written inside fyhyp.

Removed now useless intermediary procedure fxyhyper.

Changed default value of dzoomx and dzoomy from 0 to 0.2. dzoom[xy] is
only needed when grossism[xy] > 1 and then it should be > 0.

dzoom[xy] can no longer be the extent of the zoom in degrees: it must
be expressed as the fraction of the total domain.

Location:
LMDZ5/trunk/libf/dyn3d_common
Files:
2 added
1 deleted
1 edited
3 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d_common/coefpoly_m.F90

    r2217 r2218  
    1 !
    2 ! $Header$
    3 !
    4       SUBROUTINE coefpoly ( Xf1, Xf2, Xprim1, Xprim2, xtild1,xtild2 ,
    5      ,                                          a0,a1,a2,a3         )
    6       IMPLICIT NONE
    7 c
    8 c   ...  Auteur :   P. Le Van  ...
    9 c
    10 c
    11 c    Calcul des coefficients a0, a1, a2, a3 du polynome de degre 3 qui
    12 c      satisfait aux 4 equations  suivantes :
     1module coefpoly_m
    132
    14 c    a0 + a1*xtild1 + a2*xtild1*xtild1 + a3*xtild1*xtild1*xtild1 = Xf1
    15 c    a0 + a1*xtild2 + a2*xtild2*xtild2 + a3*xtild2*xtild2*xtild2 = Xf2
    16 c               a1  +     2.*a2*xtild1 +     3.*a3*xtild1*xtild1 = Xprim1
    17 c               a1  +     2.*a2*xtild2 +     3.*a3*xtild2*xtild2 = Xprim2
     3  IMPLICIT NONE
    184
    19 c  On en revient a resoudre un systeme de 4 equat.a 4 inconnues a0,a1,a2,a3
     5contains
    206
    21       REAL(KIND=8) Xf1, Xf2,Xprim1,Xprim2, xtild1,xtild2, xi
    22       REAL(KIND=8) Xfout, Xprim
    23       REAL(KIND=8) a1,a2,a3,a0, xtil1car, xtil2car,derr,x1x2car
     7  SUBROUTINE coefpoly(xf1, xf2, xprim1, xprim2, xtild1, xtild2, a0, a1, a2, a3)
    248
    25       xtil1car = xtild1 * xtild1
    26       xtil2car = xtild2 * xtild2
     9    ! From LMDZ4/libf/dyn3d/coefpoly.F, version 1.1.1.1 2004/05/19 12:53:05
    2710
    28       derr= 2. *(Xf2-Xf1)/( xtild1-xtild2)
     11    ! Author: P. Le Van
    2912
    30       x1x2car = ( xtild1-xtild2)*(xtild1-xtild2)
     13    ! Calcul des coefficients a0, a1, a2, a3 du polynôme de degré 3 qui
     14    ! satisfait aux 4 équations suivantes :
    3115
    32       a3 = (derr + Xprim1+Xprim2 )/x1x2car
    33       a2     = ( Xprim1 - Xprim2 + 3.* a3 * ( xtil2car-xtil1car ) )    /
    34      /           (  2.* ( xtild1 - xtild2 )  )
     16    ! a0 + a1 * xtild1 + a2 * xtild1**2 + a3 * xtild1**3 = Xf1
     17    ! a0 + a1 * xtild2 + a2 * xtild2**2 + a3 * xtild2**3 = Xf2
     18    ! a1 + 2. * a2 * xtild1 + 3. * a3 * xtild1**2 = Xprim1
     19    ! a1 + 2. * a2 * xtild2 + 3. * a3 * xtild2**2 = Xprim2
    3520
    36       a1     = Xprim1 -3.* a3 * xtil1car     -2.* a2 * xtild1
    37       a0     =  Xf1 - a3 * xtild1* xtil1car -a2 * xtil1car - a1 *xtild1
     21    ! (passe par les points (Xf(it), xtild(it)) et (Xf(it + 1),
     22    ! xtild(it + 1))
    3823
    39       RETURN
    40       END
     24    ! On en revient à resoudre un système de 4 équations à 4 inconnues
     25    ! a0, a1, a2, a3.
     26
     27    DOUBLE PRECISION, intent(in):: xf1, xf2, xprim1, xprim2, xtild1, xtild2
     28    DOUBLE PRECISION, intent(out):: a0, a1, a2, a3
     29
     30    ! Local:
     31    DOUBLE PRECISION xtil1car, xtil2car, derr, x1x2car
     32
     33    !------------------------------------------------------------
     34
     35    xtil1car = xtild1 * xtild1
     36    xtil2car = xtild2 * xtild2
     37
     38    derr = 2. * (xf2-xf1)/(xtild1-xtild2)
     39
     40    x1x2car = (xtild1-xtild2) * (xtild1-xtild2)
     41
     42    a3 = (derr+xprim1+xprim2)/x1x2car
     43    a2 = (xprim1-xprim2+3. * a3 * (xtil2car-xtil1car))/(2. * (xtild1-xtild2))
     44
     45    a1 = xprim1 - 3. * a3 * xtil1car - 2. * a2 * xtild1
     46    a0 = xf1 - a3 * xtild1 * xtil1car - a2 * xtil1car - a1 * xtild1
     47
     48  END SUBROUTINE coefpoly
     49
     50end module coefpoly_m

  • LMDZ5/trunk/libf/dyn3d_common/fxhyp_m.F90

    r2217 r2218  
    1 !
    2 ! $Id$
    3 !
    4 c
    5 c
    6        SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau ,
    7      , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
    8      , champmin,champmax                                               )
    9 
    10 c      Auteur :  P. Le Van
    11 
    12        IMPLICIT NONE
    13 
    14 c    Calcule les longitudes et derivees dans la grille du GCM pour une
    15 c     fonction f(x) a tangente  hyperbolique  .
    16 c
    17 c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
    18 c     dzoom  etant  la distance totale de la zone du zoom
    19 c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
    20 c
    21 c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
    22 c   ********************************************************************
    23 
    24 
    25        INTEGER nmax, nmax2
    26        PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
    27 c
    28        LOGICAL scal180
    29        PARAMETER ( scal180 = .TRUE. )
    30 
    31 c      scal180 = .TRUE.  si on veut avoir le premier point scalaire pour   
    32 c      une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
    33 c      sinon scal180 = .FALSE.
    34 
    35 #include "dimensions.h"
    36 #include "paramet.h"
    37        
    38 c     ......  arguments  d'entree   .......
    39 c
    40        REAL xzoomdeg,dzooma,tau,grossism
    41 
    42 c    ......   arguments  de  sortie  ......
    43 
    44        REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
    45      ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
    46 
    47 c     .... variables locales  ....
    48 c
    49        REAL   dzoom
    50        REAL(KIND=8) xlon(iip1),xprimm(iip1),xuv
    51        REAL(KIND=8) xtild(0:nmax2)
    52        REAL(KIND=8) fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
    53        REAL(KIND=8) Xf(0:nmax2),xxpr(0:nmax2)
    54        REAL(KIND=8) xvrai(iip1),xxprim(iip1)
    55        REAL(KIND=8) pi,depi,epsilon,xzoom,fa,fb
    56        REAL(KIND=8) Xf1, Xfi , a0,a1,a2,a3,xi2
    57        INTEGER i,it,ik,iter,ii,idif,ii1,ii2
    58        REAL(KIND=8) xi,xo1,xmoy,xlon2,fxm,Xprimin
    59        REAL(KIND=8) champmin,champmax,decalx
    60        INTEGER is2
    61        SAVE is2
    62 
    63        REAL(KIND=8) heavyside
    64 
    65        pi       = 2. * ASIN(1.)
    66        depi     = 2. * pi
    67        epsilon  = 1.e-3
    68        xzoom    = xzoomdeg * pi/180.
    69 c
    70        if (iim==1) then
    71 
    72           rlonm025(1)=-pi/2.
    73           rlonv(1)=0.
    74           rlonu(1)=pi
    75           rlonp025(1)=pi/2.
    76           rlonm025(2)=rlonm025(1)+depi
    77           rlonv(2)=rlonv(1)+depi
    78           rlonu(2)=rlonu(1)+depi
    79           rlonp025(2)=rlonp025(1)+depi
    80 
    81           xprimm025(:)=1.
    82           xprimv(:)=1.
    83           xprimu(:)=1.
    84           xprimp025(:)=1.
    85           champmin=depi
    86           champmax=depi
    87           return
    88 
    89        endif
    90 
    91            decalx   = .75
    92        IF( grossism.EQ.1..AND.scal180 )  THEN
    93            decalx   = 1.
    94        ENDIF
    95 
    96        WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx
    97 c
    98        IF( dzooma.LT.1.)  THEN
    99          dzoom = dzooma * depi
    100        ELSEIF( dzooma.LT. 25. ) THEN
    101          WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
    102      ,menter et relancer ! '
    103          STOP 1
    104        ELSE
    105          dzoom = dzooma * pi/180.
    106        ENDIF
    107 
    108        WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
    109        WRITE(6,24) xzoom,grossism,tau,dzoom
    110 
    111        DO i = 0, nmax2
    112         xtild(i) = - pi + REAL(i) * depi /nmax2
    113        ENDDO
    114 
    115        DO i = nmax, nmax2
    116 
    117        fa  = tau*  ( dzoom/2.  - xtild(i) )
    118        fb  = xtild(i) *  ( pi - xtild(i) )
    119 
    120          IF( 200.* fb .LT. - fa )   THEN
    121            fhyp ( i) = - 1.
    122          ELSEIF( 200. * fb .LT. fa ) THEN
    123            fhyp ( i) =   1.
    124          ELSE
    125             IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN
    126                 IF(   200.*fb + fa.LT.1.e-10 )  THEN
    127                     fhyp ( i ) = - 1.
    128                 ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN
    129                     fhyp ( i )  =   1.
    130                 ENDIF
    131             ELSE
    132                     fhyp ( i )  =  TANH ( fa/fb )
    133             ENDIF
    134          ENDIF
    135         IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
    136         IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
    137 
    138        ENDDO
    139 
    140 cc  ....  Calcul  de  beta  ....
    141 
    142        ffdx = 0.
    143 
    144        DO i = nmax +1,nmax2
    145 
    146        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
    147        fa  = tau*  ( dzoom/2.  - xmoy )
    148        fb  = xmoy *  ( pi - xmoy )
    149 
    150        IF( 200.* fb .LT. - fa )   THEN
    151          fxm = - 1.
    152        ELSEIF( 200. * fb .LT. fa ) THEN
    153          fxm =   1.
    154        ELSE
    155             IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN
    156                 IF(   200.*fb + fa.LT.1.e-10 )  THEN
    157                     fxm   = - 1.
    158                 ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN
    159                     fxm   =   1.
    160                 ENDIF
    161             ELSE
    162                     fxm   =  TANH ( fa/fb )
    163             ENDIF
    164        ENDIF
    165 
    166        IF ( xmoy.EQ. 0. )  fxm  =  1.
    167        IF ( xmoy.EQ. pi )  fxm  = -1.
    168 
    169        ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
    170 
    171        ENDDO
    172 
    173         beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
    174 
    175        IF( 2.*beta - grossism.LE. 0.)  THEN
    176         WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
    177      ,tine fxhyp est mauvaise ! '
    178         WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
    179      , ' et relancer ! ***  '
    180         CALL ABORT_GCM("FXHYP", "", 1)
    181        ENDIF
    182 c
    183 c   .....  calcul  de  Xprimt   .....
    184 c
    185        
    186        DO i = nmax, nmax2
    187         Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
    188        ENDDO
    189 c   
    190        DO i =  nmax+1, nmax2
    191         Xprimt( nmax2 - i ) = Xprimt( i )
    192        ENDDO
    193 c
    194 
    195 c   .....  Calcul  de  Xf     ........
    196 
    197        Xf(0) = - pi
    198 
    199        DO i =  nmax +1, nmax2
    200 
    201        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
    202        fa  = tau*  ( dzoom/2.  - xmoy )
    203        fb  = xmoy *  ( pi - xmoy )
    204 
    205        IF( 200.* fb .LT. - fa )   THEN
    206          fxm = - 1.
    207        ELSEIF( 200. * fb .LT. fa ) THEN
    208          fxm =   1.
    209        ELSE
    210          fxm =  TANH ( fa/fb )
    211        ENDIF
    212 
    213        IF ( xmoy.EQ. 0. )  fxm =  1.
    214        IF ( xmoy.EQ. pi )  fxm = -1.
    215        xxpr(i)    = beta + ( grossism - beta ) * fxm
    216 
    217        ENDDO
    218 
    219        DO i = nmax+1, nmax2
    220         xxpr(nmax2-i+1) = xxpr(i)
    221        ENDDO
    222 
    223         DO i=1,nmax2
    224          Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
    225         ENDDO
    226 
    227 
    228 c    *****************************************************************
    229 c
    230 
    231 c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
    232 c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
    233 c
    234       WRITE(6,18)
    235 c
    236       DO 5000  ik = 1, 4
    237 
    238        IF( ik.EQ.1 )        THEN
    239          xuv =  -0.25
    240        ELSE IF ( ik.EQ.2 )  THEN
    241          xuv =   0.
    242        ELSE IF ( ik.EQ.3 )  THEN
    243          xuv =   0.50
    244        ELSE IF ( ik.EQ.4 )  THEN
    245          xuv =   0.25
    246        ENDIF
    247 
    248       xo1   = 0.
    249 
    250       ii1=1
    251       ii2=iim
    252       IF(ik.EQ.1.and.grossism.EQ.1.) THEN
    253         ii1 = 2
    254         ii2 = iim+1
    255       ENDIF
    256       DO 1500 i = ii1, ii2
    257 
    258       xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim)
    259 
    260       Xfi    = xlon2
    261 c
    262       DO 250 it =  nmax2,0,-1
    263       IF( Xfi.GE.Xf(it))  GO TO 350
    264 250   CONTINUE
    265 
    266       it = 0
    267 
    268 350   CONTINUE
    269 
    270 c    ......  Calcul de   Xf(xi)    ......
    271 c
    272       xi  = xtild(it)
    273 
    274       IF(it.EQ.nmax2)  THEN
    275        it       = nmax2 -1
    276        Xf(it+1) = pi
    277       ENDIF
    278 c  .....................................................................
    279 c
    280 c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
    281 c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
    282 c          et (Xf(it+1),xtild(it+1) )
    283 
    284        CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
    285      ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
    286 
    287        Xf1     = Xf(it)
    288        Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
    289 
    290        DO 500 iter = 1,300
    291         xi = xi - ( Xf1 - Xfi )/ Xprimin
    292 
    293         IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
    294          xo1      = xi
    295          xi2      = xi * xi
    296          Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
    297          Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
    298 500   CONTINUE
    299         WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
    300           STOP 6
    301 550   CONTINUE
    302 
    303        xxprim(i) = depi/ ( REAL(iim) * Xprimin )
    304        xvrai(i)  =  xi + xzoom
    305 
    306 1500   CONTINUE
    307 
    308 
    309        IF(ik.EQ.1.and.grossism.EQ.1.)  THEN
    310          xvrai(1)    = xvrai(iip1)-depi
    311          xxprim(1)   = xxprim(iip1)
    312        ENDIF
    313 
    314        DO i = 1 , iim
    315         xlon(i)     = xvrai(i)
    316         xprimm(i)   = xxprim(i)
    317        ENDDO
    318        DO i = 1, iim -1
    319         IF( xvrai(i+1). LT. xvrai(i) )  THEN
    320          WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
    321      ,  ')'
    322         STOP 7
    323         ENDIF
    324        ENDDO
    325 c
    326 c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
    327 c   ........................................................................
    328 
    329        champmin =  1.e12
    330        champmax = -1.e12
    331        DO i = 1, iim
    332         champmin = MIN( champmin,xvrai(i) )
    333         champmax = MAX( champmax,xvrai(i) )
    334        ENDDO
    335 
    336       IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN
    337                 GO TO 1600
    338       ELSE
    339        WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
    340      ,  ' et pi '
    341 c
    342         IF( xzoom.LE.0.)  THEN
    343           IF( ik.EQ. 1 )  THEN
    344           DO i = 1, iim
    345            IF( xvrai(i).GE. - pi )  GO TO 80
    346           ENDDO
    347             WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
    348             STOP 8
    349  80       CONTINUE
    350           is2 = i
    351           ENDIF
    352 
    353           IF( is2.NE. 1 )  THEN
    354             DO ii = is2 , iim
    355              xlon  (ii-is2+1) = xvrai(ii)
    356              xprimm(ii-is2+1) = xxprim(ii)
    357             ENDDO
    358             DO ii = 1 , is2 -1
    359              xlon  (ii+iim-is2+1) = xvrai(ii) + depi
    360              xprimm(ii+iim-is2+1) = xxprim(ii)
    361             ENDDO
    362           ENDIF
    363         ELSE
    364           IF( ik.EQ.1 )  THEN
    365            DO i = iim,1,-1
    366              IF( xvrai(i).LE. pi ) GO TO 90
    367            ENDDO
    368              WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
    369               STOP 9
    370  90        CONTINUE
    371             is2 = i
    372           ENDIF
    373            idif = iim -is2
    374            DO ii = 1, is2
    375             xlon  (ii+idif) = xvrai(ii)
    376             xprimm(ii+idif) = xxprim(ii)
    377            ENDDO
    378            DO ii = 1, idif
    379             xlon (ii)  = xvrai (ii+is2) - depi
    380             xprimm(ii) = xxprim(ii+is2)
    381            ENDDO
    382          ENDIF
    383       ENDIF
    384 c
    385 c     .........   Fin  de la reorganisation   ............................
    386 
    387  1600    CONTINUE
    388 
    389 
    390          xlon  ( iip1)  = xlon(1) + depi
    391          xprimm( iip1 ) = xprimm (1 )
    392        
    393          DO i = 1, iim+1
    394          xvrai(i) = xlon(i)*180./pi
    395          ENDDO
    396 
    397          IF( ik.EQ.1 )  THEN
    398 c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
    399 c          WRITE(6,18)
    400 c          WRITE(6,68) xvrai
    401 c          WRITE(6,*) ' XPRIM k ',ik
    402 c          WRITE(6,566)  xprimm
    403 
    404            DO i = 1,iim +1
    405              rlonm025(i) = xlon( i )
    406             xprimm025(i) = xprimm(i)
    407            ENDDO
    408          ELSE IF( ik.EQ.2 )  THEN
    409 c          WRITE(6,18)
    410 c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
    411 c          WRITE(6,68) xvrai
    412 c          WRITE(6,*) ' XPRIM k ',ik
    413 c          WRITE(6,566)  xprimm
    414 
    415            DO i = 1,iim + 1
    416              rlonv(i) = xlon( i )
    417             xprimv(i) = xprimm(i)
    418            ENDDO
    419 
    420          ELSE IF( ik.EQ.3)   THEN
    421 c          WRITE(6,18)
    422 c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
    423 c          WRITE(6,68) xvrai
    424 c          WRITE(6,*) ' XPRIM ik ',ik
    425 c          WRITE(6,566)  xprimm
    426 
    427            DO i = 1,iim + 1
    428              rlonu(i) = xlon( i )
    429             xprimu(i) = xprimm(i)
    430            ENDDO
    431 
    432          ELSE IF( ik.EQ.4 )  THEN
    433 c          WRITE(6,18)
    434 c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
    435 c          WRITE(6,68) xvrai
    436 c          WRITE(6,*) ' XPRIM ik ',ik
    437 c          WRITE(6,566)  xprimm
    438 
    439            DO i = 1,iim + 1
    440              rlonp025(i) = xlon( i )
    441             xprimp025(i) = xprimm(i)
    442            ENDDO
    443 
    444          ENDIF
    445 
    446 5000    CONTINUE
    447 c
    448        WRITE(6,18)
    449 c
    450 c    ...........  fin  de la boucle  do 5000      ............
    451 
    452         DO i = 1, iim
    453          xlon(i) = rlonv(i+1) - rlonv(i)
    454         ENDDO
    455         champmin =  1.e12
    456         champmax = -1.e12
    457         DO i = 1, iim
    458          champmin = MIN( champmin, xlon(i) )
    459          champmax = MAX( champmax, xlon(i) )
    460         ENDDO
    461          champmin = champmin * 180./pi
    462          champmax = champmax * 180./pi
    463 
    464 18     FORMAT(/)
    465 24     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
    466 68     FORMAT(1x,7f9.2)
    467 566    FORMAT(1x,7f9.4)
    468 
    469        RETURN
    470        END
     1module fxhyp_m
     2
     3  IMPLICIT NONE
     4
     5contains
     6
     7  SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
     8
     9    ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
     10    ! Author: P. Le Van, from formulas by R. Sadourny
     11
     12    ! Calcule les longitudes et dérivées dans la grille du GCM pour
     13    ! une fonction f(x) à dérivée tangente hyperbolique.
     14
     15    ! Il vaut mieux avoir : grossismx \times dzoom < pi
     16
     17    ! Le premier point scalaire pour une grille regulière (grossismx =
     18    ! 1., taux=0., clon=0.) est à - 180 degrés.
     19
     20    use arth_m, only: arth
     21    use invert_zoom_x_m, only: invert_zoom_x, nmax
     22    use nrtype, only: pi, pi_d, twopi, twopi_d
     23    use principal_cshift_m, only: principal_cshift
     24
     25    include "dimensions.h"
     26    ! for iim
     27
     28    include "serre.h"
     29    ! for clon, grossismx, dzoomx, taux
     30
     31    REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
     32    real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
     33
     34    ! Local:
     35    real rlonm025(iim + 1), rlonp025(iim + 1)
     36    REAL dzoom, step
     37    real d_rlonv(iim)
     38    DOUBLE PRECISION xtild(0:2 * nmax)
     39    DOUBLE PRECISION fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax)
     40    DOUBLE PRECISION Xf(0:2 * nmax), xxpr(2 * nmax)
     41    DOUBLE PRECISION fa, fb
     42    INTEGER i, is2
     43    DOUBLE PRECISION xmoy, fxm
     44
     45    !----------------------------------------------------------------------
     46
     47    print *, "Call sequence information: fxhyp"
     48
     49    test_iim: if (iim==1) then
     50       rlonv(1)=0.
     51       rlonu(1)=pi
     52       rlonv(2)=rlonv(1)+twopi
     53       rlonu(2)=rlonu(1)+twopi
     54
     55       xprimm025(:)=1.
     56       xprimv(:)=1.
     57       xprimu(:)=1.
     58       xprimp025(:)=1.
     59    else test_iim
     60       test_grossismx: if (grossismx == 1.) then
     61          step = twopi / iim
     62
     63          xprimm025(:iim) = step
     64          xprimp025(:iim) = step
     65          xprimv(:iim) = step
     66          xprimu(:iim) = step
     67
     68          rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim)
     69          rlonm025(:iim) = rlonv(:iim) - 0.25 * step
     70          rlonp025(:iim) = rlonv(:iim) + 0.25 * step
     71          rlonu(:iim) = rlonv(:iim) + 0.5 * step
     72       else test_grossismx
     73          dzoom = dzoomx * twopi_d
     74          xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1)
     75
     76          ! Compute fhyp:
     77          DO i = nmax, 2 * nmax
     78             fa = taux * (dzoom / 2. - xtild(i))
     79             fb = xtild(i) * (pi_d - xtild(i))
     80
     81             IF (200. * fb < - fa) THEN
     82                fhyp(i) = - 1.
     83             ELSE IF (200. * fb < fa) THEN
     84                fhyp(i) = 1.
     85             ELSE
     86                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
     87                   IF (200. * fb + fa < 1e-10) THEN
     88                      fhyp(i) = - 1.
     89                   ELSE IF (200. * fb - fa < 1e-10) THEN
     90                      fhyp(i) = 1.
     91                   END IF
     92                ELSE
     93                   fhyp(i) = TANH(fa / fb)
     94                END IF
     95             END IF
     96
     97             IF (xtild(i) == 0.) fhyp(i) = 1.
     98             IF (xtild(i) == pi_d) fhyp(i) = -1.
     99          END DO
     100
     101          ! Calcul de beta
     102
     103          ffdx = 0.
     104
     105          DO i = nmax + 1, 2 * nmax
     106             xmoy = 0.5 * (xtild(i-1) + xtild(i))
     107             fa = taux * (dzoom / 2. - xmoy)
     108             fb = xmoy * (pi_d - xmoy)
     109
     110             IF (200. * fb < - fa) THEN
     111                fxm = - 1.
     112             ELSE IF (200. * fb < fa) THEN
     113                fxm = 1.
     114             ELSE
     115                IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN
     116                   IF (200. * fb + fa < 1e-10) THEN
     117                      fxm = - 1.
     118                   ELSE IF (200. * fb - fa < 1e-10) THEN
     119                      fxm = 1.
     120                   END IF
     121                ELSE
     122                   fxm = TANH(fa / fb)
     123                END IF
     124             END IF
     125
     126             IF (xmoy == 0.) fxm = 1.
     127             IF (xmoy == pi_d) fxm = -1.
     128
     129             ffdx = ffdx + fxm * (xtild(i) - xtild(i-1))
     130          END DO
     131
     132          print *, "ffdx = ", ffdx
     133          beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d)
     134          print *, "beta = ", beta
     135
     136          IF (2. * beta - grossismx <= 0.) THEN
     137             print *, 'Bad choice of grossismx, taux, dzoomx.'
     138             print *, 'Decrease dzoomx or grossismx.'
     139             STOP 1
     140          END IF
     141
     142          ! calcul de Xprimt
     143          Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp
     144          xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1)
     145
     146          ! Calcul de Xf
     147
     148          DO i = nmax + 1, 2 * nmax
     149             xmoy = 0.5 * (xtild(i-1) + xtild(i))
     150             fa = taux * (dzoom / 2. - xmoy)
     151             fb = xmoy * (pi_d - xmoy)
     152
     153             IF (200. * fb < - fa) THEN
     154                fxm = - 1.
     155             ELSE IF (200. * fb < fa) THEN
     156                fxm = 1.
     157             ELSE
     158                fxm = TANH(fa / fb)
     159             END IF
     160
     161             IF (xmoy == 0.) fxm = 1.
     162             IF (xmoy == pi_d) fxm = -1.
     163             xxpr(i) = beta + (grossismx - beta) * fxm
     164          END DO
     165
     166          xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1)
     167
     168          Xf(0) = - pi_d
     169
     170          DO i=1, 2 * nmax - 1
     171             Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1))
     172          END DO
     173
     174          Xf(2 * nmax) = pi_d
     175
     176          call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), &
     177               xprimm025(:iim), xuv = - 0.25d0)
     178          call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), &
     179               xuv = 0d0)
     180          call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), &
     181               xuv = 0.5d0)
     182          call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), &
     183               xprimp025(:iim), xuv = 0.25d0)
     184       end if test_grossismx
     185
     186       is2 = 0
     187
     188       IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
     189            .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
     190          IF (clon <= 0.) THEN
     191             is2 = 1
     192
     193             do while (rlonm025(is2) < - pi .and. is2 < iim)
     194                is2 = is2 + 1
     195             end do
     196
     197             if (rlonm025(is2) < - pi) then
     198                print *, 'Rlonm025 plus petit que - pi !'
     199                STOP 1
     200             end if
     201          ELSE
     202             is2 = iim
     203
     204             do while (rlonm025(is2) > pi .and. is2 > 1)
     205                is2 = is2 - 1
     206             end do
     207
     208             if (rlonm025(is2) > pi) then
     209                print *, 'Rlonm025 plus grand que pi !'
     210                STOP 1
     211             end if
     212          END IF
     213       END IF
     214
     215       call principal_cshift(is2, rlonm025, xprimm025)
     216       call principal_cshift(is2, rlonv, xprimv)
     217       call principal_cshift(is2, rlonu, xprimu)
     218       call principal_cshift(is2, rlonp025, xprimp025)
     219
     220       forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
     221       print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, &
     222            "degrees"
     223       print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, &
     224            "degrees"
     225
     226       ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
     227       DO i = 1, iim + 1
     228          IF (rlonp025(i) < rlonv(i)) THEN
     229             print *, 'rlonp025(', i, ') = ', rlonp025(i)
     230             print *, "< rlonv(", i, ") = ", rlonv(i)
     231             STOP 1
     232          END IF
     233
     234          IF (rlonv(i) < rlonm025(i)) THEN
     235             print *, 'rlonv(', i, ') = ', rlonv(i)
     236             print *, "< rlonm025(", i, ") = ", rlonm025(i)
     237             STOP 1
     238          END IF
     239
     240          IF (rlonp025(i) > rlonu(i)) THEN
     241             print *, 'rlonp025(', i, ') = ', rlonp025(i)
     242             print *, "> rlonu(", i, ") = ", rlonu(i)
     243             STOP 1
     244          END IF
     245       END DO
     246    end if test_iim
     247
     248  END SUBROUTINE fxhyp
     249
     250end module fxhyp_m

  • LMDZ5/trunk/libf/dyn3d_common/fyhyp_m.F90

    r2217 r2218  
    1 !
    2 ! $Id$
    3 !
    4 c
    5 c
    6        SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau  , 
    7      ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
    8      ,  champmin,champmax                                            )
    9 
    10 cc    ...  Version du 01/04/2001 ....
    11 
    12        IMPLICIT NONE
    13 c
    14 c    ...   Auteur :  P. Le Van  ...
    15 c
    16 c    .......    d'apres  formulations  de R. Sadourny  .......
    17 c
    18 c     Calcule les latitudes et derivees dans la grille du GCM pour une
    19 c     fonction f(y) a tangente  hyperbolique  .
    20 c
    21 c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
    22 c     dzoom  etant  la distance totale de la zone du zoom ( en radians )
    23 c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom   
    24 c
    25 c
    26 c N.B : Il vaut mieux avoir : grossism * dzoom  <  pi/2  (radians) ,en lati.
    27 c      ********************************************************************
    28 c
    29 c
    30 #include "dimensions.h"
    31 #include "paramet.h"
    32 
    33        INTEGER      nmax , nmax2
    34        PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
    35 c
    36 c
    37 c     .......  arguments  d'entree    .......
    38 c
    39        REAL yzoomdeg, grossism,dzooma,tau
    40 c         ( rentres  par  run.def )
    41 
    42 c     .......  arguments  de sortie   .......
    43 c
    44        REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm),
    45      , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
    46 
    47 c
    48 c     .....     champs  locaux    .....
    49 c
    50      
    51        REAL   dzoom
    52        REAL(KIND=8) ylat(jjp1), yprim(jjp1)
    53        REAL(KIND=8) yuv
    54        REAL(KIND=8) yt(0:nmax2)
    55        REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
    56        SAVE Ytprim, yt,Yf
    57        REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2)
    58        REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
    59        REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm
    60        REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
    61        REAL(KIND=8) yfi,Yf1,ffdy
    62        REAL(KIND=8) ypn,deply,y00
    63        SAVE y00, deply
    64 
    65        INTEGER i,j,it,ik,iter,jlat
    66        INTEGER jpn,jjpn
    67        SAVE jpn
    68        REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m
    69        REAL(KIND=8) fa(0:nmax2),fb(0:nmax2)
    70        REAL y0min,y0max
    71 
    72        REAL(KIND=8)     heavyside
    73 
    74        pi       = 2. * ASIN(1.)
    75        depi     = 2. * pi
    76        pis2     = pi/2.
    77        pisjm    = pi/ REAL(jjm)
    78        epsilon  = 1.e-3
    79        y0       =  yzoomdeg * pi/180.
    80 
    81        IF( dzooma.LT.1.)  THEN
    82          dzoom = dzooma * pi
    83        ELSEIF( dzooma.LT. 12. ) THEN
    84          WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
    85      ,menter et relancer ! '
    86          STOP 1
     1module fyhyp_m
     2
     3  IMPLICIT NONE
     4
     5contains
     6
     7  SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
     8
     9    ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
     10
     11    ! Author: P. Le Van, from analysis by R. Sadourny
     12
     13    ! Calcule les latitudes et dérivées dans la grille du GCM pour une
     14    ! fonction f(y) à dérivée tangente hyperbolique.
     15
     16    ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
     17
     18    use coefpoly_m, only: coefpoly
     19
     20    include "dimensions.h"
     21    ! for jjm
     22
     23    include "serre.h"
     24    ! for clat, grossismy, dzoomy, tauy
     25
     26    REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1)
     27    REAL, intent(out):: rlatv(jjm)
     28    real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)
     29
     30    ! Local:
     31
     32    DOUBLE PRECISION champmin, champmax
     33    INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
     34    REAL dzoom ! distance totale de la zone du zoom (en radians)
     35    DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
     36    DOUBLE PRECISION yuv
     37    DOUBLE PRECISION, save:: yt(0:nmax2)
     38    DOUBLE PRECISION fhyp(0:nmax2), beta
     39    DOUBLE PRECISION, save:: ytprim(0:nmax2)
     40    DOUBLE PRECISION fxm(0:nmax2)
     41    DOUBLE PRECISION, save:: yf(0:nmax2)
     42    DOUBLE PRECISION yypr(0:nmax2)
     43    DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
     44    DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm
     45    DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
     46    DOUBLE PRECISION yfi, yf1, ffdy
     47    DOUBLE PRECISION ypn, deply, y00
     48    SAVE y00, deply
     49
     50    INTEGER i, j, it, ik, iter, jlat
     51    INTEGER jpn, jjpn
     52    SAVE jpn
     53    DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m
     54    DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
     55    REAL y0min, y0max
     56
     57    DOUBLE PRECISION heavyside
     58
     59    !-------------------------------------------------------------------
     60
     61    print *, "Call sequence information: fyhyp"
     62
     63    pi = 2.*asin(1.)
     64    pis2 = pi/2.
     65    pisjm = pi/real(jjm)
     66    epsilon = 1e-3
     67    y0 = clat*pi/180.
     68    dzoom = dzoomy*pi
     69    print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
     70    print *, y0, grossismy, tauy, dzoom
     71
     72    DO i = 0, nmax2
     73       yt(i) = -pis2 + real(i)*pi/nmax2
     74    END DO
     75
     76    heavyy0m = heavyside(-y0)
     77    heavyy0 = heavyside(y0)
     78    y0min = 2.*y0*heavyy0m - pis2
     79    y0max = 2.*y0*heavyy0 + pis2
     80
     81    fa = 999.999
     82    fb = 999.999
     83
     84    DO i = 0, nmax2
     85       IF (yt(i)<y0) THEN
     86          fa(i) = tauy*(yt(i)-y0 + dzoom/2.)
     87          fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))
     88       ELSE IF (yt(i)>y0) THEN
     89          fa(i) = tauy*(y0-yt(i) + dzoom/2.)
     90          fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0)
     91       END IF
     92
     93       IF (200.*fb(i)<-fa(i)) THEN
     94          fhyp(i) = -1.
     95       ELSE IF (200.*fb(i)<fa(i)) THEN
     96          fhyp(i) = 1.
    8797       ELSE
    88          dzoom = dzooma * pi/180.
    89        ENDIF
    90 
    91        WRITE(6,18)
    92        WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
    93        WRITE(6,24) y0,grossism,tau,dzoom
    94 
    95        DO i = 0, nmax2
    96         yt(i) = - pis2  + REAL(i)* pi /nmax2
    97        ENDDO
    98 
    99        heavyy0m = heavyside( -y0 )
    100        heavyy0  = heavyside(  y0 )
    101        y0min    = 2.*y0*heavyy0m - pis2
    102        y0max    = 2.*y0*heavyy0  + pis2
    103 
    104        fa = 999.999
    105        fb = 999.999
    106        
    107        DO i = 0, nmax2
    108         IF( yt(i).LT.y0 )  THEN
    109          fa (i) = tau*  (yt(i)-y0+dzoom/2. )
    110          fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
    111         ELSEIF ( yt(i).GT.y0 )  THEN
    112          fa(i) =   tau *(y0-yt(i)+dzoom/2. )
    113          fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 )
    114        ENDIF
    115        
    116        IF( 200.* fb(i) .LT. - fa(i) )   THEN
    117          fhyp ( i) = - 1.
    118        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
    119          fhyp ( i) =   1.
    120        ELSE 
    121          fhyp(i) =  TANH ( fa(i)/fb(i) )
    122        ENDIF
    123 
    124        IF( yt(i).EQ.y0 )  fhyp(i) = 1.
    125        IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
    126 
    127        ENDDO
    128 
    129 cc  ....  Calcul  de  beta  ....
    130 c
    131        ffdy   = 0.
    132 
    133        DO i = 1, nmax2
    134         ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
    135         IF( ymoy.LT.y0 )  THEN
    136          fa(i)= tau * ( ymoy-y0+dzoom/2.)
    137          fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
    138         ELSEIF ( ymoy.GT.y0 )  THEN
    139          fa(i)= tau * ( y0-ymoy+dzoom/2. )
    140          fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
    141         ENDIF
    142 
    143         IF( 200.* fb(i) .LT. - fa(i) )    THEN
    144          fxm ( i) = - 1.
    145         ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
    146          fxm ( i) =   1.
    147         ELSE
    148          fxm(i) =  TANH ( fa(i)/fb(i) )
    149         ENDIF
    150          IF( ymoy.EQ.y0 )  fxm(i) = 1.
    151          IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
    152          ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
    153 
    154         ENDDO
    155 
    156         beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
    157 
    158        IF( 2.*beta - grossism.LE. 0.)  THEN
    159 
    160         WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
    161      ,tine fyhyp est mauvaise ! '
    162         WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
    163      , ' et relancer ! ***  '
    164         CALL ABORT_GCM("FYHYP", "", 1)
    165 
    166        ENDIF
    167 c
    168 c   .....  calcul  de  Ytprim   .....
    169 c
    170        
    171        DO i = 0, nmax2
    172         Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
    173        ENDDO
    174 
    175 c   .....  Calcul  de  Yf     ........
    176 
    177        Yf(0) = - pis2
    178        DO i = 1, nmax2
    179         yypr(i)    = beta + ( grossism - beta ) * fxm(i)
    180        ENDDO
    181 
    182        DO i=1,nmax2
    183         Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
    184        ENDDO
    185 
    186 c    ****************************************************************
    187 c
    188 c   .....   yuv  = 0.   si calcul des latitudes  aux pts.  U  .....
    189 c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
    190 c
    191       WRITE(6,18)
    192 c
    193       DO 5000  ik = 1,4
    194 
    195        IF( ik.EQ.1 )  THEN
    196          yuv  = 0.
    197          jlat = jjm + 1
    198        ELSE IF ( ik.EQ.2 )  THEN
    199          yuv  = 0.5
    200          jlat = jjm
    201        ELSE IF ( ik.EQ.3 )  THEN
    202          yuv  = 0.25
    203          jlat = jjm
    204        ELSE IF ( ik.EQ.4 )  THEN
    205          yuv  = 0.75
    206          jlat = jjm
    207        ENDIF
    208 c
    209        yo1   = 0.
    210        DO 1500 j =  1,jlat
    211         yo1   = 0.
    212         ylon2 =  - pis2 + pisjm * ( REAL(j)  + yuv  -1.) 
    213         yfi    = ylon2
    214 c
    215        DO 250 it =  nmax2,0,-1
    216         IF( yfi.GE.Yf(it))  GO TO 350
    217 250    CONTINUE
    218        it = 0
    219 350    CONTINUE
    220 
    221        yi = yt(it)
    222        IF(it.EQ.nmax2)  THEN
    223         it       = nmax2 -1
    224         Yf(it+1) = pis2
    225        ENDIF
    226 c  .................................................................
    227 c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yi) 
    228 c      .....           et   Y'(yi)                             .....
    229 c  .................................................................
    230 
    231        CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
    232      ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )     
    233 
    234        Yf1     = Yf(it)
    235        Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
    236 
    237        DO 500 iter = 1,300
    238          yi = yi - ( Yf1 - yfi )/ Yprimin
    239 
    240         IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
    241          yo1      = yi
    242          yi2      = yi * yi
    243          Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
    244          Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
    245 500   CONTINUE
    246         WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
    247          STOP 2
    248 550   CONTINUE
    249 c
    250        Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
    251        yprim(j)  = pi / ( jjm * Yprimin )
    252        yvrai(j)  = yi
    253 
    254 1500    CONTINUE
    255 
    256        DO j = 1, jlat -1
    257         IF( yvrai(j+1). LT. yvrai(j) )  THEN
    258          WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
    259      ,  ')'
    260          STOP 3
    261         ENDIF
    262        ENDDO
    263 
    264        WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
    265      , ,' et  pi/2 '
    266 c
    267         IF( ik.EQ.1 )   THEN
    268            ypn = pis2
    269           DO j = jlat,1,-1
    270            IF( yvrai(j).LE. ypn ) GO TO 1502
    271           ENDDO
    272 1502     CONTINUE
    273 
    274          jpn   = j
    275          y00   = yvrai(jpn)
    276          deply = pis2 -  y00
    277         ENDIF
    278 
    279          DO  j = 1, jjm +1 - jpn
    280            ylatt (j)  = -pis2 - y00  + yvrai(jpn+j-1)
    281            yprimm(j)  = yprim(jpn+j-1)
    282          ENDDO
    283 
    284          jjpn  = jpn
    285          IF( jlat.EQ. jjm ) jjpn = jpn -1
    286 
    287          DO j = 1,jjpn
    288           ylatt (j + jjm+1 -jpn) = yvrai(j) + deply
    289           yprimm(j + jjm+1 -jpn) = yprim(j)
    290          ENDDO
    291 
    292 c      ***********   Fin de la reorganisation     *************
    293 c
    294  1600   CONTINUE
    295 
     98          fhyp(i) = tanh(fa(i)/fb(i))
     99       END IF
     100
     101       IF (yt(i)==y0) fhyp(i) = 1.
     102       IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
     103    END DO
     104
     105    ! Calcul de beta
     106
     107    ffdy = 0.
     108
     109    DO i = 1, nmax2
     110       ymoy = 0.5*(yt(i-1) + yt(i))
     111       IF (ymoy<y0) THEN
     112          fa(i) = tauy*(ymoy-y0 + dzoom/2.)
     113          fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)
     114       ELSE IF (ymoy>y0) THEN
     115          fa(i) = tauy*(y0-ymoy + dzoom/2.)
     116          fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0)
     117       END IF
     118
     119       IF (200.*fb(i)<-fa(i)) THEN
     120          fxm(i) = -1.
     121       ELSE IF (200.*fb(i)<fa(i)) THEN
     122          fxm(i) = 1.
     123       ELSE
     124          fxm(i) = tanh(fa(i)/fb(i))
     125       END IF
     126       IF (ymoy==y0) fxm(i) = 1.
     127       IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
     128       ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
     129    END DO
     130
     131    beta = (grossismy*ffdy-pi)/(ffdy-pi)
     132
     133    IF (2. * beta - grossismy <= 0.) THEN
     134       print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
     135            // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
     136            // 'dzoomy et relancer.'
     137       STOP 1
     138    END IF
     139
     140    ! calcul de Ytprim
     141
     142    DO i = 0, nmax2
     143       ytprim(i) = beta + (grossismy-beta)*fhyp(i)
     144    END DO
     145
     146    ! Calcul de Yf
     147
     148    yf(0) = -pis2
     149    DO i = 1, nmax2
     150       yypr(i) = beta + (grossismy-beta)*fxm(i)
     151    END DO
     152
     153    DO i = 1, nmax2
     154       yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1))
     155    END DO
     156
     157    ! yuv = 0. si calcul des latitudes aux pts. U
     158    ! yuv = 0.5 si calcul des latitudes aux pts. V
     159
     160    loop_ik: DO ik = 1, 4
     161       IF (ik==1) THEN
     162          yuv = 0.
     163          jlat = jjm + 1
     164       ELSE IF (ik==2) THEN
     165          yuv = 0.5
     166          jlat = jjm
     167       ELSE IF (ik==3) THEN
     168          yuv = 0.25
     169          jlat = jjm
     170       ELSE IF (ik==4) THEN
     171          yuv = 0.75
     172          jlat = jjm
     173       END IF
     174
     175       yo1 = 0.
    296176       DO j = 1, jlat
    297           ylat(j) =  ylatt( jlat +1 -j )
    298          yprim(j) = yprimm( jlat +1 -j )
    299        ENDDO
    300  
    301         DO j = 1, jlat
    302          yvrai(j) = ylat(j)*180./pi
    303         ENDDO
    304 
    305         IF( ik.EQ.1 )  THEN
    306 c         WRITE(6,18)
    307 c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
    308 c         WRITE(6,68) (yvrai(j),j=1,jlat)
    309 cc         WRITE(6,*) ' YPRIM '
    310 cc         WRITE(6,445) ( yprim(j),j=1,jlat)
    311 
    312           DO j = 1, jlat
    313             rrlatu(j) =  ylat( j )
    314            yyprimu(j) = yprim( j )
    315           ENDDO
    316 
    317         ELSE IF ( ik.EQ. 2 )  THEN
    318 c         WRITE(6,18)
    319 c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
    320 c         WRITE(6,68) (yvrai(j),j=1,jlat)
    321 cc         WRITE(6,*)' YPRIM '
    322 cc         WRITE(6,445) ( yprim(j),j=1,jlat)
    323 
    324           DO j = 1, jlat
    325             rrlatv(j) =  ylat( j )
    326            yyprimv(j) = yprim( j )
    327           ENDDO
    328 
    329         ELSE IF ( ik.EQ. 3 )  THEN
    330 c         WRITE(6,18)
    331 c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
    332 c         WRITE(6,68) (yvrai(j),j=1,jlat)
    333 cc         WRITE(6,*) ' YPRIM '
    334 cc         WRITE(6,445) ( yprim(j),j=1,jlat)
    335 
    336           DO j = 1, jlat
    337             rlatu2(j) =  ylat( j )
    338            yprimu2(j) = yprim( j )
    339           ENDDO
    340 
    341         ELSE IF ( ik.EQ. 4 )  THEN
    342 c         WRITE(6,18)
    343 c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
    344 c         WRITE(6,68)(yvrai(j),j=1,jlat)
    345 cc         WRITE(6,*) ' YPRIM '
    346 cc         WRITE(6,68) ( yprim(j),j=1,jlat)
    347 
    348           DO j = 1, jlat
    349             rlatu1(j) =  ylat( j )
    350            yprimu1(j) = yprim( j )
    351           ENDDO
    352 
    353         ENDIF
    354 
    355 5000   CONTINUE
    356 c
    357         WRITE(6,18)
    358 c
    359 c  .....     fin de la boucle  do 5000 .....
    360 
    361         DO j = 1, jjm
    362          ylat(j) = rrlatu(j) - rrlatu(j+1)
    363         ENDDO
    364         champmin =  1.e12
    365         champmax = -1.e12
    366         DO j = 1, jjm
    367          champmin = MIN( champmin, ylat(j) )
    368          champmax = MAX( champmax, ylat(j) )
    369         ENDDO
    370          champmin = champmin * 180./pi
    371          champmax = champmax * 180./pi
    372 
    373 24     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
    374 18      FORMAT(/)
    375 68      FORMAT(1x,7f9.2)
    376 
    377         RETURN
    378         END
     177          yo1 = 0.
     178          ylon2 = -pis2 + pisjm*(real(j) + yuv-1.)
     179          yfi = ylon2
     180
     181          it = nmax2
     182          DO while (it >= 1 .and. yfi < yf(it))
     183             it = it - 1
     184          END DO
     185
     186          yi = yt(it)
     187          IF (it==nmax2) THEN
     188             it = nmax2 - 1
     189             yf(it + 1) = pis2
     190          END IF
     191
     192          ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi)
     193          ! et Y'(yi)
     194
     195          CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
     196               yt(it), yt(it + 1), a0, a1, a2, a3)
     197
     198          yf1 = yf(it)
     199          yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
     200
     201          iter = 1
     202          DO
     203             yi = yi - (yf1-yfi)/yprimin
     204             IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit
     205             yo1 = yi
     206             yi2 = yi*yi
     207             yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi
     208             yprimin = a1 + 2.*a2*yi + 3.*a3*yi2
     209          END DO
     210          if (abs(yi-yo1) > epsilon) then
     211             print *, 'Pas de solution.', j, ylon2
     212             STOP 1
     213          end if
     214
     215          yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
     216          yprim(j) = pi/(jjm*yprimin)
     217          yvrai(j) = yi
     218       END DO
     219
     220       DO j = 1, jlat - 1
     221          IF (yvrai(j + 1)<yvrai(j)) THEN
     222             print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', &
     223                  j, ')'
     224             STOP 1
     225          END IF
     226       END DO
     227
     228       print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2'
     229
     230       IF (ik==1) THEN
     231          ypn = pis2
     232          DO j = jjm + 1, 1, -1
     233             IF (yvrai(j)<=ypn) exit
     234          END DO
     235
     236          jpn = j
     237          y00 = yvrai(jpn)
     238          deply = pis2 - y00
     239       END IF
     240
     241       DO j = 1, jjm + 1 - jpn
     242          ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1)
     243          yprimm(j) = yprim(jpn + j-1)
     244       END DO
     245
     246       jjpn = jpn
     247       IF (jlat==jjm) jjpn = jpn - 1
     248
     249       DO j = 1, jjpn
     250          ylatt(j + jjm + 1-jpn) = yvrai(j) + deply
     251          yprimm(j + jjm + 1-jpn) = yprim(j)
     252       END DO
     253
     254       ! Fin de la reorganisation
     255
     256       DO j = 1, jlat
     257          ylat(j) = ylatt(jlat + 1-j)
     258          yprim(j) = yprimm(jlat + 1-j)
     259       END DO
     260
     261       DO j = 1, jlat
     262          yvrai(j) = ylat(j)*180./pi
     263       END DO
     264
     265       IF (ik==1) THEN
     266          DO j = 1, jjm + 1
     267             rlatu(j) = ylat(j)
     268             yyprimu(j) = yprim(j)
     269          END DO
     270       ELSE IF (ik==2) THEN
     271          DO j = 1, jjm
     272             rlatv(j) = ylat(j)
     273          END DO
     274       ELSE IF (ik==3) THEN
     275          DO j = 1, jjm
     276             rlatu2(j) = ylat(j)
     277             yprimu2(j) = yprim(j)
     278          END DO
     279       ELSE IF (ik==4) THEN
     280          DO j = 1, jjm
     281             rlatu1(j) = ylat(j)
     282             yprimu1(j) = yprim(j)
     283          END DO
     284       END IF
     285    END DO loop_ik
     286
     287    DO j = 1, jjm
     288       ylat(j) = rlatu(j) - rlatu(j + 1)
     289    END DO
     290    champmin = 1e12
     291    champmax = -1e12
     292    DO j = 1, jjm
     293       champmin = min(champmin, ylat(j))
     294       champmax = max(champmax, ylat(j))
     295    END DO
     296    champmin = champmin*180./pi
     297    champmax = champmax*180./pi
     298
     299    DO j = 1, jjm
     300       IF (rlatu1(j) <= rlatu2(j)) THEN
     301          print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
     302          STOP 13
     303       ENDIF
     304
     305       IF (rlatu2(j) <= rlatu(j+1)) THEN
     306          print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
     307          STOP 14
     308       ENDIF
     309
     310       IF (rlatu(j) <= rlatu1(j)) THEN
     311          print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
     312          STOP 15
     313       ENDIF
     314
     315       IF (rlatv(j) <= rlatu2(j)) THEN
     316          print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
     317          STOP 16
     318       ENDIF
     319
     320       IF (rlatv(j) >= rlatu1(j)) THEN
     321          print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
     322          STOP 17
     323       ENDIF
     324
     325       IF (rlatv(j) >= rlatu(j)) THEN
     326          print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
     327          STOP 18
     328       ENDIF
     329    ENDDO
     330
     331    print *, 'Latitudes'
     332    print 3, champmin, champmax
     333
     3343   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
     335         ' d environ ', f0.2, ' degres ', /, &
     336         ' alors que la maille en dehors de la zone du zoom est ', &
     337         "d'environ ", f0.2, ' degres ')
     338
     339  END SUBROUTINE fyhyp
     340
     341end module fyhyp_m

  • LMDZ5/trunk/libf/dyn3d_common/inigeom.F

    r1952 r2218  
    1616c
    1717c
     18      use fxhyp_m, only: fxhyp
     19      use fyhyp_m, only: fyhyp
    1820      IMPLICIT NONE
    1921c
    … …  
    264266      WRITE(6,*)'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
    265267 
    266        CALL fxyhyper( clat, grossismy, dzoomy, tauy    ,
    267      ,                clon, grossismx, dzoomx, taux    ,
    268      , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2  ,
    269      , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
    270 
     268      CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
     269      CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
    271270 
    272271      ENDIF
Note: See TracChangeset for help on using the changeset viewer.