Changeset 214 in lmdz_wrf for trunk/tools

Timestamp:

Jan 12, 2015, 6:56:17 PM (10 years ago)

Author:

lfita

Message:

Trying to avoid hyperbolic-tangent for the zoom...

Location:

trunk/tools

Files:

• iniaqua.py (modified) (25 diffs)
• lmdz_const.py (modified) (3 diffs)

trunk/tools/iniaqua.py

@@ -19 +19 @@
 filekinds = ['CF', 'WRF']
 
-## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z param 
+## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo
+def fxy(dx, dy): 
+    """! 
+       ! $Id: fxy.F 1403 2010-07-01 09:02:53Z fairhead $ 
+       ! 
+       c     Auteur  :  P. Le Van 
+       c 
+       c     Calcul  des longitudes et des latitudes  pour une fonction f(x,y) 
+       c           a tangente sinusoidale et eventuellement avec zoom  . 
+       c 
+       c 
+    """ 
+    fname ='fxy' 
+ 
+#c    ......  calcul  des  latitudes  et de y'   ..... 
+#c 
+    vrlatu = np.zeros((dy+1), dtype=np.float) 
+    vyprimu = np.zeros((dy+1), dtype=np.float) 
+         
+    vrlatu = ffy(np.arange(dy+1)*1. + 1.) 
+    vyprimu = ffy(np.arange(dy+1)*1. + 1.) 
+         
+    vrlatv = np.zeros((dy), dtype=np.float) 
+    vrlatu1 = np.zeros((dy), dtype=np.float) 
+    vrlatu2 = np.zeros((dy), dtype=np.float) 
+    yprimv = np.zeros((dy), dtype=np.float) 
+    vyprimu1 = np.zeros((dy), dtype=np.float) 
+    vyprimu2 = np.zeros((dy), dtype=np.float) 
+         
+    vrlatv = ffy(np.arange(dy)*1. + 1. + 0.5) 
+    vrlatu1 = ffy(np.arange(dy)*1. + 1. + 0.25) 
+    vrlatu2 = ffy(np.arange(dy)*1. + 1. + 0.75) 
+         
+    vyprimv = fyprim(np.arange(dy)*1. + 1. + 0.5) 
+    vyprimu1 = fyprim(np.arange(dy)*1. + 1. + 0.25) 
+    vyprimu2 = fyprim(np.arange(dy)*1. + 1. + 0.75) 
+         
+#c 
+#c     .....  calcul   des  longitudes et de  x'   ..... 
+#c 
+    vrlonv = np.zeros((dx+1), dtype=np.float) 
+    vrlonu = np.zeros((dx+1), dtype=np.float) 
+    vrlonm025 = np.zeros((dx+1), dtype=np.float) 
+    vrlonp025 = np.zeros((dx+1), dtype=np.float) 
+    vxprimv = np.zeros((dx+1), dtype=np.float) 
+    vxprimu = np.zeros((dx+1), dtype=np.float) 
+    vxprimm025 = np.zeros((dx+1), dtype=np.float) 
+    vxprimo025 = np.zeros((dx+1), dtype=np.float) 
+         
+    vrlonv = ffy(np.arange(dx+1)*1. + 1.) 
+    vrlonu = ffy(np.arange(dx+1)*1. + 1. + 0.5) 
+    vrlonm025 = ffy(np.arange(dx+1)*1. + 1. - 0.25) 
+    vrlonp025 = ffy(np.arange(dx+1)*1. + 1. + 0.25) 
+         
+    vxprimv = fxprim(np.arange(dx+1)*1. + 1.) 
+    vxprimu = fxprim(np.arange(dx+1)*1. + 1. + 0.5) 
+    vxprimm025 = fxprim(np.arange(dx+1)*1. + 1. - 0.25) 
+    vxprimp025 = fxprim(np.arange(dx+1)*1. + 1. + 0.25) 
+         
+    return vrlatu, vyprimu, vrlatv, vyprimv, vrlatu1, vyprimu1, vrlatu2, vyprimu2,   \
+      vrlonu, vxprimu, vrlonv, vxprimv, vrlonm025, vxprimm025, vrlonp025, vxprimp025 
+
+# From grid/fxy_sin.h
+#    ................................................................
+#    ................  Fonctions in line  ...........................
+#    ................................................................
+#
+
+def fy(rj,dy):
+   """
+   """
+   val = np.arcsin(1.+2.*((1.-rj)/np.float(dy)))
+
+   return val
+
+def fyprim(rj,dy):
+    """
+    """
+    val = 1./np.sqrt((rj-1.)*(dy+1.-rj))
+
+    return val
+
+def fx(ri,dx):
+    """
+    """
+    val = 2.*np.pi/np.float(dx) * ( ri - 0.5*np.float(dx) - 1. )
+
+    return val
+
+def fxprim(ri,dx):
+    """
+    """
+    val = 2.*np.pi/np.float(dx)
+
+    return val
+
+#
+#
+#    La valeur de pi est passee par le common/const/ou /const2/ .
+#    Sinon, il faut la calculer avant d'appeler ces fonctions .
+#
+#   ----------------------------------------------------------------
+#     Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
+#   -----------------------------------------------------------------
+#
+#    .....  ici, on a l'application particuliere suivante   ........
+#
+#                **************************************
+#                **     x = 2. * pi/iim *  X         **
+#                **     y =      pi/jjm *  Y         **
+#                **************************************
+#
+#   ..................................................................
+#   ..................................................................
+#
+#
+#
+#-----------------------------------------------------------------------
+
+def fxysinus (ddy,ddx):
+    """c
+       c     Calcul  des longitudes et des latitudes  pour une fonction f(x,y)
+       c            avec y = Asin( j )  .
+       c
+       c     Auteur  :  P. Le Van
+       c
+       c
+      from: dyn3d/fxysinus.F
+    """
+    fname = 'fxysinus'
+
+#    ......  calcul  des  latitudes  et de y'   .....
+#
+    rlatu = np.zeros((ddy+1), dtype=np.float)
+    yprimu = np.zeros((ddy+1), dtype=np.float)
+
+    for j in range(ddy+1):
+        rlatu[j] = fy( np.float(j+1), ddy)
+        yprimu[j] = fyprim( np.float(j+1), ddy)
+
+    rlatv = np.zeros((ddy), dtype=np.float)
+    rlatu1 = np.zeros((ddy), dtype=np.float)
+    rlatu2 = np.zeros((ddy), dtype=np.float)
+    yprimv = np.zeros((ddy), dtype=np.float)
+    yprimu1 = np.zeros((ddy), dtype=np.float)
+    yprimu2 = np.zeros((ddy), dtype=np.float)
+
+    for j in range(ddy):
+        rlatv[j] = fy( np.float(j) + 0.5, ddy)
+        rlatu1[j] = fy( np.float(j) + 0.25, ddy) 
+        rlatu2[j] = fy( np.float(j) + 0.75, ddy) 
+
+        yprimv[j] = fyprim( np.float(j) + 0.5, ddy) 
+        yprimu1[j] = fyprim( np.float(j) + 0.25, ddy)
+        yprimu2[j] = fyprim( np.float(j) + 0.75, ddy)
+
+#
+#     .....  calcul   des  longitudes et de  x'   .....
+#
+    rlonv = np.zeros((ddx+1), dtype=np.float)
+    rlonu = np.zeros((ddx+1), dtype=np.float)
+    rlonm025 = np.zeros((ddx+1), dtype=np.float)
+    rlonp025 = np.zeros((ddx+1), dtype=np.float)
+    xprimv = np.zeros((ddx+1), dtype=np.float)
+    xprimu = np.zeros((ddx+1), dtype=np.float)
+    xprimm025 = np.zeros((ddx+1), dtype=np.float)
+    xprimp025 = np.zeros((ddx+1), dtype=np.float)
+
+    for i in range(ddx + 1):
+        rlonv[i] = fx( np.float(i), ddx )
+        rlonu[i] = fx( np.float(i)+ 0.5, ddx )
+        rlonm025[i] = fx( np.float(i)- 0.25, ddx )
+        rlonp025[i] = fx( np.float(i)+ 0.25, ddx )
+
+        xprimv[i] = fxprim(np.float(i), ddx )
+        xprimu[i] = fxprim(np.float(i)+ 0.5, ddx )
+        xprimm025[i] = fxprim(np.float(i)- 0.25, ddx )
+        xprimp025[i] = fxprim(np.float(i)+ 0.25, ddx )
+
+    return rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu,    \
+      xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025
+
+def fxhyp (dx, dy, xzoomdeg, grossism):
+    """
+       c      Auteur :  P. Le Van 
+
+
+       c    Calcule les longitudes et derivees dans la grille du GCM pour une
+       c     fonction f(x) a tangente  hyperbolique  .
+       c
+       c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
+       c     dzoom  etant  la distance totale de la zone du zoom
+       c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
+       c
+       c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
+       c   ********************************************************************
+    """
+    fname = 'fxhyp'
+
+    nmax = 30000
+    nmax2 = 2*nmax
+
+    scal180 = True
+
+#      scal180 = .TRUE.  si on veut avoir le premier point scalaire pour   
+#      une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
+#      sinon scal180 = .FALSE.
+
+#     ......  arguments  d'entree   .......
+#
+#       REAL xzoomdeg,dzooma,tau,grossism
+
+#    ......   arguments  de  sortie  ......
+    rlonm025 = np.zeros((dx+1), dtype=np.float)
+    xprimm025 = np.zeros((dx+1), dtype=np.float)
+    rlonv = np.zeros((dx+1), dtype=np.float)
+    xprimv = np.zeros((dx+1), dtype=np.float)
+    rlonu = np.zeros((dx+1), dtype=np.float)
+    xprimu = np.zeros((dx+1), dtype=np.float)
+    rlonp025 = np.zeros((dx+1), dtype=np.float)
+    xprimp025 = np.zeros((dx+1), dtype=np.float)
+
+#     .... variables locales  ....
+#
+    xlon = np.zeros((dx+1), dtype=np.float32),
+    xprimm = np.zeros((dx+1), dtype=np.float32)
+    xtild = np.zeros((nmax2), dtype=np.float32)
+    fhyp = np.zeros((nmax2), dtype=np.float32)
+    Xprimt = np.zeros((nmax2), dtype=np.float32)
+    Xf = np.zeros((nmax2), dtype=np.float32)
+    xxpr = np.zeros((nmax2), dtype=np.float32)
+    xvrai = np.zeros((dx+1), dtype=np.float32)
+    xxprim = np.zeros((dx+1), dtype=np.float32) 
+    
+    pi = np.pi
+    depi = 2. * np.pi
+    epsilon = 1.e-3
+    xzoom = xzoomdeg * pi/180. 
+
+    if dx == 1:
+        rlonm025[1] = -pi/2.
+        rlonv[1] = 0.
+        rlonu[1] = pi
+        rlonp025[1] = pi/2.
+        rlonm025[2] = rlonm025[1] + depi
+        rlonv[2] = rlonv[1] + depi
+        rlonu[2] = rlonu[1] + depi
+        rlonp025[2] = rlonp025[1] + depi
+
+        xprimm025 = 1.
+        xprimv = 1.
+        xprimu = 1.
+        xprimp025 = 1.
+        champmin =depi
+        champmax = depi
+        return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu,        \
+         rlonp025, xprimp025, champmin, champmax
+
+    decalx   = .75
+    if grossism == 1. and scal180:
+        decalx   = 1.
+
+    print 'FXHYP scal180,decalx', scal180,decalx
+
+    if dzooma > 1.:
+        dzoom = dzooma * depi
+    elif dzooma < 25.:
+        print erormsg
+        print '  ' + fname +": Le param. dzoomx pour 'fxhyp' est trop petit dzooma:",\
+          dzooma,'!L augmenter et relancer !'
+        quit(-1)
+    else:
+        dzoom = dzooma * pi/180.
+
+
+    print ' xzoom( rad.),grossism,tau,dzoom (radians)'
+    print xzoom,grossism,tau,dzoom
+
+    xtild = - pi + np.range(namx2) * depi /nmax2
+
+    for i in range(nmax,nmax2):
+       fa  = tau*  ( dzoom/2.  - xtild[i] )
+       fb  = xtild[i] *  ( pi - xtild[i] )
+
+       if 200.* fb < - fa:
+           fhyp[i] = -1.
+       elif 200. * fb < fa:
+          fhyp[i] = 1.
+       else:
+           if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13:
+               if 200.*fb + fa < 1.e-10:
+                   fhyp[i] = -1.
+               elif 200.*fb - fa < 1.e-10:
+                   fhyp[i] = 1.
+           else:
+               fhyp[i] = np.tanh(fa/fb)
+       if xtild[i] == 0.: fhyp[i] =  1.
+       if xtild[i] == pi: fhyp[i] = -1.
+
+##  ....  Calcul  de  beta  ....
+
+    ffdx = 0.
+
+    for i in range(nmax+1,nmax2):
+        xmoy = 0.5 * ( xtild[i-1] + xtild[i] )
+        fa = tau*( dzoom/2. - xmoy )
+        fb = xmoy*( pi - xmoy )
+
+        if 200.* fb < -fa:
+            fxm = -1.
+        elif 200. * fb < fa:
+            fxm = 1.
+        else:
+            if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13:
+                if 200.*fb + fa < 1.e-10:
+                    fxm = -1.
+                elif 200.*fb - fa < 1.e-10:
+                    fxm = 1.
+            else:
+                fxm = np.tanh( fa/fb )
+
+        if xmoy == 0.: fxm = 1.
+        if xmoy == pi: fxm = -1.
+
+        ffdx = ffdx + fxm * ( xtild[i] - xtild[i-1] )
+
+    beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
+
+    if 2.*beta - grossism < 0.:
+        print warnmsg
+        print '  ' + fname + ': **  Attention ! La valeur beta calculee dans la ' +   \
+          'routine fxhyp est mauvaise ! '
+        print '    Modifier les valeurs de  grossismx ,tau ou dzoomx et relancer ! ***'
+        quit(-1)
+
+#
+#   .....  calcul  de  Xprimt   .....
+#
+       
+    for i in range(nmax, nmax2):
+        Xprimt[i] = beta  + ( grossism - beta ) * fhyp[i]
+
+    for i in range(nmax+1, nmax2):
+        Xprimt[nmax2-i] = Xprimt[i]
+
+#   .....  Calcul  de  Xf     ........
+
+    Xf[0] = - pi
+
+    for i in range(nmax+1, nmax2):
+        xmoy = 0.5 * ( xtild[i-1] + xtild[i] )
+        fa  = tau*  ( dzoom/2.  - xmoy )
+        fb  = xmoy *  ( pi - xmoy )
+
+        if 200.* fb < -fa:
+            fxm = -1.
+        elif 200. * fb < fa:
+            fxm = 1.
+        else:
+            fxm = np.tanh( fa/fb )
+
+        if xmoy == 0.: fxm =  1.
+        if xmoy == pi: fxm = -1.
+        xxpr[i] = beta + ( grossism - beta ) * fxm
+
+    for i in range(nmax+1, nmax2):
+        xxpr[nmax2-i+1] = xxpr[i]
+
+    for i in range(nmax2):
+         Xf[i]   = Xf[i-1] + xxpr[i] * ( xtild[i] - xtild[i-1] )
+
+#    *****************************************************************
+#
+
+#     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
+#     .....  xuv = 0.5  si  calcul  aux pts      U        ........
+
+    for ik in range(4):
+        if ik == 0:
+            xuv = -0.25
+        elif ik == 1:
+            xuv = 0.
+        elif ik == 2:
+            xuv = 0.50
+        elif ik == 4:
+            xuv = 0.25
+
+        xo1   = 0.
+
+        ii1=1
+        ii2=dx
+        if ik == 0 and grossism == 1.:
+            ii1 = 2 
+            ii2 = dx+1
+
+        for i in range(ii1, ii2):
+            xlon2 = - pi + (np.float(i) + xuv - decalx) * depi / np.float32(dx) 
+
+            Xfi = xlon2
+
+            ended = False
+            for it in range(nmax2,0,-1):
+                if Xfi > Xf[it]:
+                    ended = True
+                    break
+            if not ended: it = 0
+
+#    ......  Calcul de   Xf(xi)    ...... 
+#
+        xi  = xtild[it]
+
+        if it == nmax2:
+            it = nmax2 -1
+            Xf[it+1] = pi
+
+#  .....................................................................
+#
+#   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
+#   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
+#          et (Xf(it+1),xtild(it+1) )
+
+        a0, a1, a2, a3 = coefpoly( Xf[it], Xf[it+1], Xprimt[it], Xprimt[it+1],       \
+          xtild[it], xtild[it+1])
+
+        Xf1 = Xf[it]
+        Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
+
+        for iteri in range(300):
+            xi = xi - ( Xf1 - Xfi )/ Xprimin
+
+            ended = False
+            if np.abs(xi-xo1) < epsilon: 
+                ended = True
+                break
+            xo1 = xi
+            xi2 = xi * xi
+            Xf1 = a0 +  a1 * xi + a2 * xi2  + a3 * xi2 * xi
+            Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
+
+            if not ended:
+                print errmsg
+                print '  ' + fname + ': Pas de solution ***** ',i,xlon2,iteri
+                quit(-1)
+
+        xxprim[i] = depi/ ( np.float32(dx) * Xprimin )
+        xvrai[i]  =  xi + xzoom
+
+    if ik == 0 and grossism == 1:
+        xvrai[0] = xvrai[dx+1] - depi
+        xxprim[0] = xxprim[dx+1]
+
+    xlon[0:dx+1] = xvrai[0:dx+1]
+    xprimm[0:dx+1] = xxprim[0:dx+1]
+
+    for i in range(dx-1):
+        if xvrai[i+1] < xvrai[i]:
+            print errormsg
+            print '  ' + fname + ': PBS. avec rlonu(',i+1,') plus petit que rlonu(', \
+              i,')'
+            quit(-1)
+
+#
+#   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
+#   ........................................................................
+
+    champmin =  1.e12
+    champmax = -1.e12
+    for i in range(dx):
+        champmin = np.min( [champmin,xvrai[i]] )
+        champmax = np.max( [champmax,xvrai[i]] )
+
+#    if champmin >= -pi-0.10 and champmax <= pi+0.10:
+#                GO TO 1600
+
+#
+# HERE -- here
+#
+#      ELSE
+#       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
+#     ,  ' et pi '
+#c
+#        IF( xzoom.LE.0.)  THEN
+#          IF( ik.EQ. 1 )  THEN
+#          DO i = 1, iim
+#           IF( xvrai(i).GE. - pi )  GO TO 80
+#          ENDDO
+#            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
+#            STOP 8
+# 80       CONTINUE
+#          is2 = i
+#          ENDIF
+#
+#          IF( is2.NE. 1 )  THEN
+#            DO ii = is2 , iim
+#             xlon  (ii-is2+1) = xvrai(ii)
+#             xprimm(ii-is2+1) = xxprim(ii)
+#            ENDDO
+#            DO ii = 1 , is2 -1
+#             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
+#             xprimm(ii+iim-is2+1) = xxprim(ii) 
+#            ENDDO
+#          ENDIF
+#        ELSE 
+#          IF( ik.EQ.1 )  THEN
+#           DO i = iim,1,-1
+#             IF( xvrai(i).LE. pi ) GO TO 90
+#           ENDDO
+#             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
+#              STOP 9
+# 90        CONTINUE
+#            is2 = i
+#          ENDIF
+#           idif = iim -is2
+#           DO ii = 1, is2
+#            xlon  (ii+idif) = xvrai(ii)
+#            xprimm(ii+idif) = xxprim(ii)
+#           ENDDO
+#           DO ii = 1, idif
+#            xlon (ii)  = xvrai (ii+is2) - depi
+#            xprimm(ii) = xxprim(ii+is2) 
+#           ENDDO
+#         ENDIF
+#      ENDIF
+#c
+#c     .........   Fin  de la reorganisation   ............................
+
+# 1600    CONTINUE
+
+
+#         xlon  ( iip1)  = xlon(1) + depi
+#         xprimm( iip1 ) = xprimm (1 )
+       
+#         DO i = 1, iim+1
+#         xvrai(i) = xlon(i)*180./pi
+#         ENDDO
+
+#         IF( ik.EQ.1 )  THEN
+#c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
+#c          WRITE(6,18) 
+#c          WRITE(6,68) xvrai
+#c          WRITE(6,*) ' XPRIM k ',ik
+#c          WRITE(6,566)  xprimm
+
+#           DO i = 1,iim +1
+#             rlonm025(i) = xlon( i )
+#            xprimm025(i) = xprimm(i)
+#           ENDDO
+#         ELSE IF( ik.EQ.2 )  THEN
+#c          WRITE(6,18) 
+#c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
+#c          WRITE(6,68) xvrai
+#c          WRITE(6,*) ' XPRIM k ',ik
+#c          WRITE(6,566)  xprimm
+
+#           DO i = 1,iim + 1
+#             rlonv(i) = xlon( i )
+#            xprimv(i) = xprimm(i)
+#           ENDDO
+
+#         ELSE IF( ik.EQ.3)   THEN
+#c          WRITE(6,18) 
+#c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
+#c          WRITE(6,68) xvrai
+#c          WRITE(6,*) ' XPRIM ik ',ik
+#c          WRITE(6,566)  xprimm
+
+#           DO i = 1,iim + 1
+#             rlonu(i) = xlon( i )
+#            xprimu(i) = xprimm(i)
+#           ENDDO
+
+#         ELSE IF( ik.EQ.4 )  THEN
+#c          WRITE(6,18) 
+#c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
+#c          WRITE(6,68) xvrai
+#c          WRITE(6,*) ' XPRIM ik ',ik
+#c          WRITE(6,566)  xprimm
+
+#           DO i = 1,iim + 1
+#             rlonp025(i) = xlon( i )
+#            xprimp025(i) = xprimm(i)
+#           ENDDO
+
+#         ENDIF
+
+#5000    CONTINUE
+#c
+#       WRITE(6,18)
+#c
+#c    ...........  fin  de la boucle  do 5000      ............
+
+#        DO i = 1, iim
+#         xlon(i) = rlonv(i+1) - rlonv(i)
+#        ENDDO
+#        champmin =  1.e12
+#        champmax = -1.e12
+#        DO i = 1, iim
+#         champmin = MIN( champmin, xlon(i) )
+#         champmax = MAX( champmax, xlon(i) )
+#        ENDDO
+#         champmin = champmin * 180./pi
+#         champmax = champmax * 180./pi
+
+#18     FORMAT(/)
+#24     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
+#68     FORMAT(1x,7f9.2)
+#566    FORMAT(1x,7f9.4)
+
+    return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu,           \
+      rlonp025, xprimp025, champmin, champmax
+
+
+def fxyhyper (ddy, ddx, yzoom, grossy, dzoomy, tauy, xzoom, grossx, dzoomx, taux,    \
+     rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, 
+     rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025):
+    """c
+       c      Auteur :  P. Le Van .
+       c
+       c      d'apres  formulations de R. Sadourny .
+       c
+       c
+       c     Ce spg calcule les latitudes( routine fyhyp ) et longitudes( fxhyp )
+       c            par des  fonctions  a tangente hyperbolique .
+       c
+       c     Il y a 3 parametres ,en plus des coordonnees du centre du zoom (xzoom
+       c                      et  yzoom )   :  
+       c
+       c     a) le grossissement du zoom  :  grossy  ( en y ) et grossx ( en x )
+       c     b) l' extension     du zoom  :  dzoomy  ( en y ) et dzoomx ( en x )
+       c     c) la raideur de la transition du zoom  :   taux et tauy   
+       c
+       c  N.B : Il vaut mieux avoir   :   grossx * dzoomx <  pi    ( radians )
+       c ******
+       c                  et              grossy * dzoomy <  pi/2  ( radians )
+       c
+    """
+    fname = 'fxyhyper'
+
+#       CALL fyhyp ( yzoom, grossy, dzoomy,tauy  , 
+#     ,  rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
+#     ,  dymin,dymax                                               )
+
+#       CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
+#     , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax         )
+
+
+    for i in range(ddx+1):
+        if rlonp025[i] < rlonv[i]:
+            print errormsg
+            print '  ' + fname + ' Attention !  rlonp025 < rlonv',i
+            quit(-1)
+
+        if rlonv[i] < rlonm025[i]:
+            print errormsg
+            print '  ' + fname + ' Attention !  rlonm025 > rlonv',i
+            quit(-1)
+
+        if rlonp025[i] > rlonu[i]:
+            print errormsg
+            print '  ' + fname + ' Attention !  rlonp025 > rlonu',i
+            quit(-1)
+
+    print '  *** TEST DE COHERENCE  OK    POUR   FX **** '
+
+
+    for j in range(ddy):
+        if rlatu1[j] <= rlatu2[j]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatu1 < rlatu2',rlatu1[j],rlatu2[j],j
+            quit(-1)
+
+        if rlatu2[j] <= rlatu[j+1]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatu2 < rlatup1',rlatu2[j],rlatu[j+1],j
+            quit(-1)
+
+        if rlatu[j] <= rlatu1[j]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatu < rlatu1',rlatu[j],rlatu1[j],j
+            quit(-1)
+
+        if rlatv[j] <= rlatu2[j]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatv < rlatu2 ',rlatv[j],rlatu2[j],j
+            quit(-1)
+
+        if rlatv[j] >= rlatu1[j]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatv > rlatu1 ',rlatv[j],rlatu1[j],j
+            quit(-1)
+
+        if rlatv[j] >= rlatu[j]:
+            print errormsg
+            print '  ' + fname + 'Attention ! rlatv > rlatu ',rlatv[j],rlatu[j],j
+            quit(-1)
+
+    print '  *** TEST DE COHERENCE  OK    POUR   FY **** '
+
+    print '  Latitudes  '
+    print ' *********** '
+
+    print 'Au centre du zoom, la longueur de la maille est d environ', dymin ,       \
+     'degres alors que la maille en dehors de la zone du zoom est d environ',        \
+     dymax, 'degres'
+    print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\
+      'tau , dzoom pour Y et repasser ! '
+
+    print '  Longitudes  '
+    print ' ************ '
+    print 'Au centre du zoom, la longueur de la maille est d environ', dxmin ,       \
+     'degres alors que la maille en dehors de la zone du zoom est d environ',        \
+     dxmax, 'degres'
+    print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\
+      'tau , dzoom pour Y et repasser ! '
+
+    return
+
+def inigeom(dy,dx):
+    """c
+       c     Auteur :  P. Le Van
+       c
+       c   ............      Version  du 01/04/2001     ........................
+       c
+       c  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-
+       c     endroits que les aires aireij1,..aireij4 .
+
+       c  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
+       c
+       c
+    """
+    fname = 'inigeom'
+
+    cvu = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuv = np.zeros((dy, dx+1), dtype=np.float)
+
+    cuij1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuij2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuij3 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuij4 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvij1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvij2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvij3 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvij4 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aireij1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aireij2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aireij3 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aireij4 = np.zeros((dy+1,dx+1), dtype=np.float)
+
+    aire = np.zeros((dy+1,dx+1), dtype=np.float)
+    aireu = np.zeros((dy+1,dx+1), dtype=np.float)
+    airev = np.zeros((dy+1,dx+1), dtype=np.float)
+    unsaire = np.zeros((dy+1,dx+1), dtype=np.float)
+    unsair_gam1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    unsair_gam2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    airesurg = np.zeros((dy+1,dx+1), dtype=np.float)
+    unsairez = np.zeros((dy+1,dx+1), dtype=np.float)
+    unsairz_gam = np.zeros((dy+1,dx+1), dtype=np.float)
+    fext = np.zeros((dy+1,dx+1), dtype=np.float)
+
+    alpha1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha3 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha4 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha1p2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha1p4 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha2p3 = np.zeros((dy+1,dx+1), dtype=np.float)
+    alpha3p4 = np.zeros((dy+1,dx+1), dtype=np.float)
+
+    rlonvv = np.zeros((dx+1), dtype=np.float)
+    rlatuu = np.zeros((dy+1), dtype=np.float)
+    rlatu1 = np.zeros((dy), dtype=np.float)
+    yprimu1 = np.zeros((dy), dtype=np.float)
+    rlatu2 = np.zeros((dy), dtype=np.float)
+    yprimu2 = np.zeros((dy), dtype=np.float)
+    yprimv = np.zeros((dy), dtype=np.float)
+    yprimu = np.zeros((dy+1), dtype=np.float)
+
+    rlonm025 = np.zeros((dx+1), dtype=np.float)
+    xprimm025 = np.zeros((dx+1), dtype=np.float)
+    rlonp025 = np.zeros((dx+1), dtype=np.float)
+    xprimp025 = np.zeros((dx+1), dtype=np.float)
+
+    cu = np.zeros((dy+1,dx+1), dtype=np.float)
+    cv = np.zeros((dy+1,dx+1), dtype=np.float)
+    unscu2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    unscv2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuvsurcv = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvsurcv = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuvscvgam1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuvscvgam2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvscuvgam = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvusurcu = np.zeros((dy+1,dx+1), dtype=np.float)
+    cusurcvu = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvuscugam1 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cvuscugam2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    cuscvugam = np.zeros((dy+1,dx+1), dtype=np.float)
+    airvscu2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aivscu2gam = np.zeros((dy+1,dx+1), dtype=np.float)
+    airuscv2 = np.zeros((dy+1,dx+1), dtype=np.float)
+    aiuscv2gam = np.zeros((dy+1,dx+1), dtype=np.float)
+
+    constang = np.zeros((dy+1,dx+1), dtype=np.float)
+#
+#
+#   ------------------------------------------------------------------
+#   -                                                                -
+#   -    calcul des coeff. ( cu, cv , 1./cu**2,  1./cv**2  )         -
+#   -                                                                -
+#   ------------------------------------------------------------------
+#
+#      les coef. ( cu, cv ) permettent de passer des vitesses naturelles
+#      aux vitesses covariantes et contravariantes , ou vice-versa ...
+#
+#
+#     on a :  u (covariant) = cu * u (naturel)   , u(contrav)= u(nat)/cu
+#             v (covariant) = cv * v (naturel)   , v(contrav)= v(nat)/cv
+#
+#       on en tire :  u(covariant) = cu * cu * u(contravariant)
+#                     v(covariant) = cv * cv * v(contravariant)
+#
+#
+#     on a l'application (  x(X) , y(Y) )   avec - im/2 +1 <  X  < im/2
+#                                                          =     =
+#                                           et   - jm/2    <  Y  < jm/2
+#                                                          =     =
+#
+#      ...................................................
+#      ...................................................
+#      .  x  est la longitude du point  en radians       .
+#      .  y  est la  latitude du point  en radians       .
+#      .                                                 .
+#      .  on a :  cu(i,j) = rad * COS(y) * dx/dX         .
+#      .          cv( j ) = rad          * dy/dY         .
+#      .        aire(i,j) =  cu(i,j) * cv(j)             .
+#      .                                                 .
+#      . y, dx/dX, dy/dY calcules aux points concernes   .
+#      .                                                 .
+#      ...................................................
+#      ...................................................
+#
+#
+#
+#                                                           ,
+#    cv , bien que dependant de j uniquement,sera ici indice aussi en i
+#          pour un adressage plus facile en  ij  .
+#
+#
+#
+#  **************  aux points  u  et  v ,           *****************
+#      xprimu et xprimv sont respectivement les valeurs de  dx/dX
+#      yprimu et yprimv    .  .  .  .  .  .  .  .  .  .  .  dy/dY
+#      rlatu  et  rlatv    .  .  .  .  .  .  .  .  .  .  .la latitude
+#      cvu    et   cv      .  .  .  .  .  .  .  .  .  .  .    cv
+#
+#  **************  aux points u, v, scalaires, et z  ****************
+#      cu, cuv, cuscal, cuz sont respectiv. les valeurs de    cu
+#
+#
+#
+#         Exemple de distribution de variables sur la grille dans le
+#             domaine de travail ( X,Y ) .
+#     ................................................................
+#                  DX=DY= 1
+#
+#   
+#        +     represente  un  point scalaire ( p.exp  la pression )
+#        >     represente  la composante zonale du  vent
+#        V     represente  la composante meridienne du vent
+#        o     represente  la  vorticite
+#
+#       ----  , car aux poles , les comp.zonales covariantes sont nulles
+#
+#
+#
+#         i ->
+#
+#         1      2      3      4      5      6      7      8
+#  j
+#  v  1   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
+#
+#         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
+#
+#     2   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
+#
+#         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
+#
+#     3   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
+#
+#         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
+#
+#     4   +   >  +   >  +   >  +   >  +   >  +   >  +   >  +  >
+#
+#         V   o  V   o  V   o  V   o  V   o  V   o  V   o  V  o
+#
+#     5   + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
+#
+#
+#      Ci-dessus,  on voit que le nombre de pts.en longitude est egal
+#                 a   IM = 8
+#      De meme ,   le nombre d'intervalles entre les 2 poles est egal
+#                 a   JM = 4
+#
+#      Les points scalaires ( + ) correspondent donc a des valeurs
+#       entieres  de  i ( 1 a IM )   et  de  j ( 1 a  JM +1 )   .
+#
+#      Les vents    U       ( > ) correspondent a des valeurs  semi-
+#       entieres  de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
+#
+#      Les vents    V       ( V ) correspondent a des valeurs entieres
+#       de     i ( 1 a  IM ) et semi-entieres de  j ( 1 +0.5  a JM +0.5)
+#
+#
+#
+
+    if nitergdiv != 2:
+        gamdi_gdiv = coefdis/ ( np.float(nitergdiv) -2. )
+    else:
+        gamdi_gdiv = 0.
+
+    if nitergrot != 2:
+        gamdi_grot = coefdis/ ( np.float(nitergrot) -2. )
+    else:
+        gamdi_grot = 0.
+
+    if niterh != 2:
+        gamdi_h = coefdis/ ( np.float(niterh) -2. )
+    else:
+        gamdi_h = 0.
+
+    print 'gamdi_gd:',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,nitergdiv,nitergrot,niterh
+
+#     ----------------------------------------------------------------
+#
+    if not fxyhypb:
+        if ysinus:
+            print ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
+
+#   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
+
+            rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu,   \
+              xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 =      \
+              fxysinus(dx, dy)
+        else:
+            print '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
+
+#  .... utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ...
+ 
+            pxo   = clon *np.pi /180.
+            pyo   = 2.* clat* np.pi /180.
+#
+#  ....  determination de  transx ( pour le zoom ) par Newton-Raphson ...
+#
+            itmax = 10
+            eps = .1e-7
+
+            xo1 = 0.
+            for iter in range(itmax):
+                x1 = xo1
+                f = x1+ alphax*np.sin(x1-pxo)
+                df = 1.+ alphax*np.cos(x1-pxo)
+                x1 = x1 - f/df
+                xdm = np.abs( x1- xo1 )
+                if xdm > eps: xo1 = x1
+
+            transx = xo1
+
+            itmay = 10
+            eps   = .1e-7
+
+            yo1  = 0.
+            for iter in range(itmay):
+                y1 = yo1
+                f = y1 + alphay*np.sin(y1-pyo)
+                df = 1. + alphay*np.cos(y1-pyo)
+                y1 = y1 -f/df
+                ydm = np.abs(y1-yo1)
+                if ydm > eps: yo1 = y1
+
+            transy = yo1
+
+            rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu,   \
+              xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 =      \
+              fxy(dx,dy)
+    else:
+#
+#   ....  Utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
+#   .....................................................................
+
+        print '*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
+ 
+#!!!!! Lluis
+#!! HERE, how to do all this without zoom!?
+#!
+
+#        CALL fxyhyper( clat, grossismy, dzoomy, tauy    , 
+#     ,                clon, grossismx, dzoomx, taux    ,
+#     , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2  ,
+#     , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
+
+#  -------------------------------------------------------------------
+
+    rlatu[0] = np.arcsin(1.)
+    rlatu[dy] = -rlatu[0]
+
+#   ....  calcul  aux  poles  ....
+
+    yprimu[0] = 0.
+    yprimu[dy] = 0.
+
+    un4rad2 = 0.25 * rad * rad
+
+#
+#   --------------------------------------------------------------------
+#   --------------------------------------------------------------------
+#   -                                                                  -
+#   -  calcul  des aires ( aire,aireu,airev, 1./aire, 1./airez  )      -
+#   -      et de   fext ,  force de coriolis  extensive  .             -
+#   -                                                                  -
+#   --------------------------------------------------------------------
+#   --------------------------------------------------------------------
+#
+#
+#
+#   A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
+#   affectees 4 aires entourant P , calculees respectivement aux points
+#            ( i + 1/4, j - 1/4 )    :    aireij1 (i,j)
+#            ( i + 1/4, j + 1/4 )    :    aireij2 (i,j)
+#            ( i - 1/4, j + 1/4 )    :    aireij3 (i,j)
+#            ( i - 1/4, j - 1/4 )    :    aireij4 (i,j)
+#
+#           ,
+#   Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
+#   Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
+#   des 4 aires  aireij1,aireij2,aireij3,aireij4 qui sont affectees au
+#   point (i,j) .
+#   On definit en outre les coefficients  alpha comme etant egaux a
+#    (aireij / aire), c.a.d par exp.  alpha1(i,j)=aireij1(i,j)/aire(i,j)
+#
+#   De meme, toute aire centree en 1 point U est egale a la somme des
+#   4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
+#    Idem pour  airev, airez .
+#
+#       On a ,pour chaque maille :    dX = dY = 1
+#
+#
+#                             . V
+#
+#                 aireij4 .        . aireij1
+#
+#                   U .       . P      . U
+#
+#                 aireij3 .        . aireij2
+#
+#                             . V
+#
+#
+#
+#
+#
+#   ....................................................................
+#
+#    Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
+#   qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
+#   taires cuij et les 4 elongat. cvij qui sont calculees aux memes 
+#     endroits  que les aireij   .    
+#
+#   ....................................................................
+#
+#     .......  do 35  :   boucle sur les  jjm + 1  latitudes   .....
+#
+#
+    for j in range(dy+1):
+        if j == 1:
+            yprm = yprimu1[j]
+            rlatm = rlatu1[j]
+
+            coslatm = np.cos( rlatm )
+            radclatm = 0.5* rad * coslatm
+
+            for i in range(dx):
+                xprp = xprimp025[i]
+                xprm = xprimm025[i]
+                aireij2[0,i] = un4rad2 * coslatm  * xprp * yprm
+                aireij3[0,i] = un4rad2 * coslatm  * xprm * yprm
+                cuij2[0,i] = radclatm * xprp
+                cuij3[0,i] = radclatm * xprm
+                cvij2[0,i] = 0.5* rad * yprm
+                cvij3[0,i] = cvij2[0,i]
+
+            for i in range(dx):
+                aireij1[0,i] = 0.
+                aireij4[0,i] = 0.
+                cuij1[0,i] = 0.
+                cuij4[0,i] = 0.
+                cvij1[0,i] = 0.
+                cvij4[0,i] = 0.
+
+        elif j == dy:
+            yprp = yprimu2[j-1]
+            rlatp = rlatu2[j-1]
+
+            coslatp = np.cos( rlatp )
+            radclatp = 0.5* rad * coslatp
+
+            for i in range(dx):
+                xprp = xprimp025[i]
+                xprm = xprimm025[i]
+                aireij1[dy,i] = un4rad2 * coslatp  * xprp * yprp
+                aireij4[dy,i] = un4rad2 * coslatp  * xprm * yprp
+                cuij1[dy,i] = radclatp * xprp
+                cuij4[dy,i] = radclatp * xprm
+                cvij1[dy,i] = 0.5 * rad* yprp
+                cvij4[dy,i] = cvij1[dy,i]
+
+            for i in range(dx):
+                aireij2[dy,i] = 0.
+                aireij3[dy,i] = 0.
+                cvij2[dy,i] = 0.
+                cvij3[dy,i] = 0.
+                cuij2[dy,i] = 0.
+                cuij3[dy,i] = 0.
+
+        else:
+
+            rlatp = rlatu2[j-1]
+            yprp = yprimu2[j-1]
+            rlatm = rlatu1[j]
+            yprm = yprimu1[j]
+
+            coslatm = np.cos( rlatm )
+            coslatp = np.cos( rlatp )
+            radclatp = 0.5* rad * coslatp
+            radclatm = 0.5* rad * coslatm
+
+            for i in range(dx):
+                xprp = xprimp025[i]
+                xprm = xprimm025[i]
+      
+                ai14 = un4rad2 * coslatp * yprp
+                ai23 = un4rad2 * coslatm * yprm
+                aireij1[j,i] = ai14 * xprp
+                aireij2[j,i] = ai23 * xprp
+                aireij3[j,i] = ai23 * xprm
+                aireij4[j,i] = ai14 * xprm
+                cuij1[j,i] = radclatp * xprp
+                cuij2[j,i] = radclatm * xprp
+                cuij3[j,i] = radclatm * xprm
+                cuij4[j,i] = radclatp * xprm
+                cvij1[j,i] = 0.5* rad * yprp
+                cvij2[j,i] = 0.5* rad * yprm
+                cvij3[j,i] = cvij2[j,i]
+                cvij4[j,i] = cvij1[j,i]
+
+#
+#    ........       periodicite   ............
+#
+        cvij1[j,dx] = cvij1[j,0]
+        cvij2[j,dx] = cvij2[j,0]
+        cvij3[j,dx] = cvij3[j,0]
+        cvij4[j,dx] = cvij4[j,0]
+        cuij1[j,dx] = cuij1[j,0]
+        cuij2[j,dx] = cuij2[j,0]
+        cuij3[j,dx] = cuij3[j,0]
+        cuij4[j,dx] = cuij4[j,0]
+        aireij1[j,dx] = aireij1[j,0]
+        aireij2[j,dx] = aireij2[j,0]
+        aireij3[j,dx] = aireij3[j,0]
+        aireij4[j,dx] = aireij4[j,0]
+        
+#
+#    ..............................................................
+#
+    for j in range(dy+1):
+        for i in range(dx):
+            aire[j,i] = aireij1[j,i] + aireij2[j,i] + aireij3[j,i] + aireij4[j,i]
+            alpha1[j,i] = aireij1[j,i] / aire[j,i]
+            alpha2[j,i] = aireij2[j,i] / aire[j,i]
+            alpha3[j,i] = aireij3[j,i] / aire[j,i]
+            alpha4[j,i] = aireij4[j,i] / aire[j,i]
+            alpha1p2[j,i] = alpha1 [j,i] + alpha2 [j,i]
+            alpha1p4[j,i] = alpha1 [j,i] + alpha4 [j,i]
+            alpha2p3[j,i] = alpha2 [j,i] + alpha3 [j,i]
+            alpha3p4[j,i] = alpha3 [j,i] + alpha4 [j,i]
+
+        aire[j,dx] = aire[j,0]
+        alpha1[j,dx] = alpha1[j,0]
+        alpha2[j,dx] = alpha2[j,0]
+        alpha3[j,dx] = alpha3[j,0]
+        alpha4[j,dx] = alpha4[j,0]
+        alpha1p2[j,dx] = alpha1p2[j,0]
+        alpha1p4[j,dx] = alpha1p4[j,0]
+        alpha2p3[j,dx] = alpha2p3[j,0]
+        alpha3p4[j,dx] = alpha3p4[j,0]
+
+    for j in range(dy+1):
+        for i in range(dx):
+            aireu[j,i] = aireij1[j,i] + aireij2[j,i] + aireij4[j,i+1] + aireij3[j,i+1]
+            unsaire[j,i] = 1./ aire[j,i]
+            unsair_gam1[j,i] = unsaire[j,i]** ( - gamdi_gdiv )
+            unsair_gam2[j,i] = unsaire[j,i]** ( - gamdi_h    )
+            airesurg[j,i] = aire[j,i]/ g
+
+        aireu[j,dx] = aireu[j,0]
+        unsaire[j,dx] = unsaire[j,0]
+        unsair_gam1[j,dx] = unsair_gam1[j,0]
+        unsair_gam2[j,dx] = unsair_gam2[j,0]
+        airesurg[j,dx] = airesurg[j,0]
+
+
+    for j in range(dy):
+        for i in range(dx):
+            airev[j,i] = aireij2[j,i]+ aireij3[j,i]+ aireij1[j,i] + aireij4[j+1,i]
+
+        for i in range(dx):
+            airez = aireij2[j,i]+aireij1[j+1,i]+aireij3[j,i+1] + aireij4[j+1,i+1]
+            unsairez[j,i] = 1./ airez
+            unsairz_gam[j,i]= unsairez[j,i]** ( - gamdi_grot )
+            fext[j,i] = airez * np.sin(rlatv[j])* 2.* omeg
+
+        airev[j,dx] = airev[j,0]
+        unsairez[j,dx] = unsairez[j,0]
+        fext[j,dx] = fext[j,0]
+        unsairz_gam[j,dx] = unsairz_gam[j,0]
+
+
+#
+#    .....      Calcul  des elongations cu,cv, cvu     .........
+#
+    for j in range(dy):
+        for i in range(dx):
+            cv[j,i] = 0.5 *( cvij2[j,i]+cvij3[j,i]+cvij1[j+1,i]+cvij4[j+1,i])
+            cvu[j,i]= 0.5 *( cvij1[j,i]+cvij4[j,i]+cvij2[j,i]+cvij3[j,i] )
+            cuv[j,i]= 0.5 *( cuij2[j,i]+cuij3[j,i]+cuij1[j+1,i]+cuij4[j+1,i])
+            unscv2[j,i] = 1./ ( cv[j,i]*cv[j,i] )
+
+        for i in range(dx):
+            cuvsurcv [j,i] = airev[j,i] * unscv2[j,i]
+            cvsurcuv [j,i] = 1./cuvsurcv[j,i]
+            cuvscvgam1[j,i] = cuvsurcv [j,i] ** ( - gamdi_gdiv )
+            cuvscvgam2[j,i] = cuvsurcv [j,i] ** ( - gamdi_h )
+            cvscuvgam[j,i] = cvsurcuv [j,i] ** ( - gamdi_grot )
+
+        cv[j,dx] = cv[j,0]
+        cvu[j,dx] = cvu[j,0]
+        unscv2[j,dx] = unscv2[j,0]
+        cuv[j,dx] = cuv[j,0]
+        cuvsurcv[j,dx] = cuvsurcv[j,0]
+        cvsurcuv[j,dx] = cvsurcuv[j,0]
+        cuvscvgam1[j,dx] = cuvscvgam1[j,0]
+        cuvscvgam2[j,dx] = cuvscvgam2[j,0]
+        cvscuvgam[j,dx] = cvscuvgam[j,0]
+
+
+    for j in range(1,dy):
+        for i in range(dx):
+            cu[j,i] = 0.5*(cuij1[j,i]+cuij4[j,i+1]+cuij2[j,i]+cuij3[j,i+1])
+            unscu2[j,i] = 1./ ( cu[j,i] * cu[j,i] )
+            cvusurcu[j,i] = aireu[j,i] * unscu2[j,i]
+            cusurcvu[j,i] = 1./ cvusurcu[j,i]
+            cvuscugam1[j,i] = cvusurcu[j,i] ** ( - gamdi_gdiv ) 
+            cvuscugam2[j,i] = cvusurcu[j,i] ** ( - gamdi_h    ) 
+            cuscvugam[j,i] = cusurcvu[j,i] ** ( - gamdi_grot )
+
+        cu[j,dx]  = cu[j,0]
+        unscu2[j,dx] = unscu2[j,0]
+        cvusurcu[j,dx] = cvusurcu[j,0]
+        cusurcvu[j,dx] = cusurcvu[j,0]
+        cvuscugam1[j,dx] = cvuscugam1[j,0]
+        cvuscugam2[j,dx] = cvuscugam2[j,0]
+        cuscvugam[j,dx] = cuscvugam[j,0]
+
+#
+#   ....  calcul aux  poles  ....
+#
+    for i in range(dx+1):
+        cu[0, i] = 0.
+        unscu2[0, i] = 0.
+        cvu[0, i] = 0.
+
+        cu[dy,i] = 0.
+        unscu2[dy,i] = 0.
+        cvu[dy,i] = 0.
+
+#
+#    ..............................................................
+#
+    for j in range(dy):
+        for i in range(dx):
+           airvscu2[j,i] = airev[j,i]/ ( cuv[j,i] * cuv[j,i] )
+           aivscu2gam[j,i] = airvscu2[j,i]** ( - gamdi_grot )
+
+        airvscu2[j,dx] = airvscu2[j,0]
+        aivscu2gam[j,dx] = aivscu2gam[j,0]
+
+    for j in range(dy):
+        for i in range(dx):
+            airuscv2[j,i] = aireu[j,i]/ ( cvu[j,i] * cvu[j,i] )
+            aiuscv2gam[j,i] = airuscv2[j,i]** ( - gamdi_grot ) 
+
+        airuscv2[j,dx]  = airuscv2[j,0]
+        aiuscv2gam[j,dx]  = aiuscv2gam[j,0]
+
+#
+#   calcul des aires aux  poles :
+#   -----------------------------
+#
+    apoln = SSUM(dx,aire[0,0],1)
+    apols = SSUM(dx,aire[dy,0],1)
+    unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
+    unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
+    unsapolnga2 = 1./ ( apoln ** ( - gamdi_h    ) )
+    unsapolsga2 = 1./ ( apols ** ( - gamdi_h    ) )
+
+#
+#-----------------------------------------------------------------------
+#     gtitre='Coriolis version ancienne'
+#     gfichier='fext1'
+#     CALL writestd(fext,iip1*jjm)
+#
+#   changement F. Hourdin calcul conservatif pour fext
+#   constang contient le produit a * cos ( latitude ) * omega
+#
+    for i in range(dx):
+        constang[0,i] = 0.
+
+    for j in range(dy-1):
+        for i in range(dx):
+            constang[j+1,i] = rad*omeg*cu[j+1,i]*np.cos(rlatu[j+1])
+
+    for i in range(dx):
+        constang[dy,i] = 0.
+
+#
+#   periodicite en longitude
+#
+    for j in range(dy):
+        fext[j,dx] = fext[j,0]
+
+    for j in range(dy+1):
+        constang[j,dx] = constang[j,0]
+
+# fin du changement
+
+#
+#-----------------------------------------------------------------------
+#
+    print '   ***  Coordonnees de la grille  *** '
+    print '   LONGITUDES  aux pts.   V  ( degres )  '
+
+    for i in range(dx+1):
+        rlonvv[i] = rlonv[i]*180./np.pi
+
+    print rlonvv
+
+    print '   LATITUDES   aux pts.   V  ( degres )  '
+
+    for i in range(dy):
+        rlatuu[i] = rlatv[i]*180./np.pi
+
+    print rlatuu[dy]
+
+    for i in range(dx+1):
+        rlonvv[i]=rlonu[i]*180./np.pi
+
+    print '   LONGITUDES  aux pts.   U  ( degres )  '
+    print rlonvv
+
+
+    print '   LATITUDES   aux pts.   U  ( degres )  '
+    for i in range(dy+1):
+        rlatuu[i]=rlatu[i]*180./np.pi
+
+    print rlatuu[0:dy+2]
+
+    return aire
 
 def SSUM(n,sx,incx):
@@ -30 +1406 @@
     ix = 0
 
-    for i in range(n):
-        ssumv=ssumv+sx[ix]
-        ix=ix+incx
+    if len(sx.shape) != 0:
+        for i in range(n):
+            ssumv=ssumv+sx[ix]
+            ix=ix+incx
+    else:
+        ssumv=ssumv+sx
 
     return ssumv
@@ -50 +1429 @@
     return sy
 
-def exner_hyb (dx, dy, dz, psv, pv, alphav, betav):
+def exner_hyb (dx, dy, dz, psv, pv, aire):
     """c
        c     Auteurs :  P.Le Van  , Fr. Hourdin  .
@@ -83 +1462 @@
     pk = np.zeros((dz, dy, dx), dtype=np.float)
     pkf = np.zeros((dz, dy, dx), dtype=np.float)
+
+    ppn = np.zeros((dy, dx), dtype=np.float)
+    pps = np.zeros((dy, dx), dtype=np.float)
+
+    ip1jm = (dx+1)*dy
 
     if dz == 1:       
@@ -101 +1485 @@
     unpl2k = 1.+ 2.* kappa
 
-    for ij in range(ngrid):
-        pks[ij] = cpp * ( ps[ij]/preff ) ** kappa
-
-    for ij in range(iim):
-        ppn[ij] = aire[ij] * pks[ij]
-        pps[ij] = aire[ij+ip1jm] * pks[ij+ip1jm]
+    for j in range(dy):
+        for i in range(dx):
+            pks[j,i] = cpp * ( ps[j,i]/preff ) ** kappa
+
+    for j in range(dy):
+        for i in range(dx):
+            ppn[j,i] = aire[j,i] * pks[j,i]
+            pps[j,i] = aire[j,i+ip1jm] * pks[ij+ip1jm]
 
     xpn = SSUM(iim,ppn,1) /apoln
@@ -152 +1538 @@
 #      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
 
-    return pksv, pkv, pkfv
+    return pksv, pkv, pkfv, aplhav, betav
 
 def exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ):
@@ -187 +1573 @@
     ppn = np.zeros((iim), dtype=np.float)
     ppn = np.zeros((iim), dtype=np.float)
+
+    ip1jm = (dx+1)*dy
 
     firstcall = True
@@ -204 +1592 @@
                 quit(-1)
 
-
         firstcall = False
-
 
 #### Specific behaviour for Shallow Water (1 vertical layer) case:
     if llm == 1:
-      
+     
 #         Compute pks(:),pk(:),pkf(:)
         
@@ -216 +1602 @@
             pks[ij] = (cpp/preff) * ps[ij]
             pk[ij,1] = .5*pks[ij]
-
           
         pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
@@ -291 +1676 @@
     fname = 'pression'
       
-    press = np.zeros((dz, dy, dx), dtype=np.float)
-
+    press = np.zeros((dz+1, dy, dx), dtype=np.float)
+
+    print 'shape psv:',psv.shape,'press:',press.shape,'ap:',apv.shape,'bp:',bpv.shape
     for l in range(dz+1):
-        press[l,:,:] = apv[l] + bpv[l]*psv[:,:]
+        press[l,:,:] = apv[l] + bpv[l]*psv
    
     return press
@@ -368 +1754 @@
     dsig = np.zeros((dz), dtype=np.float)
     dpres = np.zeros((dz), dtype=np.float)
-    ap = np.zeros((dz+1), dtype=np.float)
-    bp = np.zeros((dz+1), dtype=np.float)
+    apv = np.zeros((dz+1), dtype=np.float)
+    bpv = np.zeros((dz+1), dtype=np.float)
 
 # default scaleheight is 8km for earth
@@ -412 +1798 @@
         sig[dz-1]=0.
 
-        bp[0:dz] = np.exp(1.-1./sig[0:dz]**2)
-        bp[llmp1-1] = 0.
-
-        ap = pa * (sig - bp)
+        bpv[0:dz] = np.exp(1.-1./sig[0:dz]**2)
+        bpv[llmp1-1] = 0.
+
+        apv = pa * (sig - bp)
 
     elif vert_sampling  == 'tropo':
@@ -427 +1813 @@
             sig[l] = sig[l+1] + dsig[l]
 
-        bp[0]=1.
-        bp[1:dz] = np.exp(1.-1./sig[1:dz]**2)
-        bp[llmp1-1] = 0.
-
-        ap[0] = 0.
-        ap[1:dz+1] = pa*(sig[1:dz+1]-bp[1:dz+1])
+        bpv[0]=1.
+        bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
+        bpv[llmp1-1] = 0.
+
+        apv[0] = 0.
+        apv[1:dz+1] = pa*(sig[1:dz+1]-bpv[1:dz+1])
 
     elif vert_sampling  == 'strato':
@@ -455 +1841 @@
             sig[l] = sig[l+1] + dsig[l]
 
-        bp[0] = 1.
-        bp[1:dz] = np.exp(1.-1./sig[1:dz]**2)
-        bp[llmp1-1] = 0.
-
-        ap[0] = 0.
-        ap[1:dz+1] = pa*(sig[1:dz+1] - bp[1:dz+1])
+        bpv[0] = 1.
+        bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
+        bpv[llmp1-1] = 0.
+
+        apv[0] = 0.
+        apv[1:dz+1] = pa*(sig[1:dz+1] - bpv[1:dz+1])
 
     elif vert_sampling  == 'read':
@@ -475 +1861 @@
         for l in range(dz+1):
             values = ncvar.readuce_space(sfobj.readline())
-            ap[l] = np.float(values[0])
-            bp[l] = np.float(values[0])
+            apv[l] = np.float(values[0])
+            bpv[l] = np.float(values[0])
 
         sfobj.close()
-        if ap[0] == 0. or ap[dz+1] == 0. or bp[0] == 1. or bp[dz+1] == 0.:
+        if apv[0] == 0. or apv[dz+1] == 0. or bpv[0] == 1. or bpv[dz+1] == 0.:
             print errormsg
             print '  ' + fname + ': bad ap or bp values !!'
             print '    k ap bp ___________'
             for k in range(dz+1):
-                print k, ap[k], bp[k]
+                print k, apv[k], bpv[k]
           
     else:
@@ -497 +1883 @@
     print '  ' + fname + ': k BP AP ________'
     for k in range(dz+1):
-        print k, bp[k], ap[k]
+        print k, bpv[k], apv[k]
 
     print 'Niveaux de pressions approximatifs aux centres des'
@@ -505 +1891 @@
 
     for l in range(dz):
-        dpres[l] = bp[l] - bp[l+1]
-        pnivs[l] = 0.5*( ap[l]+bp[l]*preff + ap[l+1]+bp[l+1]*preff )
+        dpres[l] = bpv[l] - bpv[l+1]
+        pnivs[l] = 0.5*( apv[l]+bpv[l]*preff + apv[l+1]+bpv[l+1]*preff )
         print '    PRESNIVS(', l, ')=', pnivs[l], ' Z ~ ',                           \
           np.log(preff/pnivs[l])*scaleheight, ' DZ ~ ',                              \
-          scaleheight*np.log((ap[l]+bp[l]*preff)/                                    \
-          np.max([ap[l+1]+bp[l+1]*preff, 1.e-10]))
+          scaleheight*np.log((apv[l]+bpv[l]*preff)/                                  \
+          np.max([apv[l+1]+bpv[l+1]*preff, 1.e-10]))
 
     print '  ' + fname + ': PRESNIVS [Pa]:', pnivs
 
-    return pnivs
+    return pnivs, apv, bpv
 
 
@@ -536 +1922 @@
     s = np.zeros((dz), dtype=np.float)
     sig = np.zeros((dz+1), dtype=np.float)
-    ap = np.zeros((dz+1), dtype=np.float)
-    bp = np.zeros((dz+1), dtype=np.float)
+    apv = np.zeros((dz+1), dtype=np.float)
+    bpv = np.zeros((dz+1), dtype=np.float)
 
 # Ouverture possible de fichiers typiquement E.T.
@@ -635 +2021 @@
         for l in range(dz):
             newsig = sig_hybrid(sig[l],pa,preff)
-            bp[l] = np.exp(1.-1./(newsig**2))
-            ap[l] = pa * (newsig - bp[l] )
+            bpv[l] = np.exp(1.-1./(newsig**2))
+            apv[l] = pa * (newsig - bp[l] )
 
         bp[llmp1-1] = 0.
@@ -647 +2033 @@
 #       Pour ne pas passer en coordonnees hybrides
         for l in range(dz):
-            ap[l] = 0.
-            bp[l] = sig[l]
-
-        ap[llmp1-1] = 0.
+            apv[l] = 0.
+            bpv[l] = sig[l]
+
+        apv[llmp1-1] = 0.
 
     bp[llmp1-1] = 0.
@@ -666 +2052 @@
 
     for l in range(dz-1):
-        aps[0:dz-1] = 0.5*( ap[0:dz-1]+ap[1:dz]) 
-        bps[0:dz-1] = 0.5*( bp[0:dz-1]+bp[1:dz]) 
+        aps[0:dz-1] = 0.5*( apv[0:dz-1]+apv[1:dz]) 
+        bps[0:dz-1] = 0.5*( bpv[0:dz-1]+bpv[1:dz]) 
 
     if hybrid:
         aps[dz-1] = aps[dz-2]**2 / aps[dz-3] 
-        bps[dz-1] = 0.5*(bp[dz-1] + bp[dz])
+        bps[dz-1] = 0.5*(bpv[dz-1] + bpv[dz])
     else:
         bps[dz-1] = bps[dz-2]**2 / bps[dz-3] 
@@ -1063 +2449 @@
 # For extra-terrestrial planets
 #presnivs, pseudoalt = presnivs_calc_noterre(opts.znivs, dimz)
-presnivs = presnivs_calc(opts.znivs, dimz)
+presnivs, ap, bp = presnivs_calc(opts.znivs, dimz)
 lon, lat = global_lonlat(dimy,dimx)
 lonu, latu = global_lonlat(dimy,dimx+1)
@@ -1119 +2505 @@
 # 3. Initialize fields (if necessary)
 # surface pressure
-ps =  np.zeros((dimy, dimx), dtype=np.float)
 
 if iflag_phys > 2:
@@ -1125 +2510 @@
 # "Specify the initial dry mass to be equivalent to
 #  a global mean surface pressure (101325 minus 245) Pa."
-    ps = 101080.  
+    ps = np.ones((dimy, dimx), dtype=np.float)*101080.  
 else:
 # use reference surface pressure
-    ps = preff
+    ps = np.ones((dimy, dimx), dtype=np.float)*preff
         
 # ground geopotential
@@ -1135 +2520 @@
 p = pression(dimx, dimy, dimz, ap, bp, ps)
 
+aire = inigeom(dimx, dimy)
+
 if opts.pexner == 'hybdrid':
-    pks, pk, pkf = exner_hyb(ip1jmp1, ps, p, alpha, beta)
+    pks, pk, pkf, alpha, beta = exner_hyb(dimx, dimy, dimz, ps, p, aire)
 else:
-    pks, pk, pkf = exner_milieu(ip1jmp1, ps, p, beta)
+    pks, pk, pkf = exner_milieu(ip1jmp1, dimx, dimy, dimz, ps, p, beta)
 
 masse = massdair(p)

trunk/tools/lmdz_const.py

@@ -66 +66 @@
 a = 6371229.
 r1sa = np.float(np.float64(1.)/np.float64(a))
+#rad = 6371220 
+rad = a*1.
+omeg = 7.272205e-05
+
 #
 #     -----------------------------------------------------------------
@@ -92 +96 @@
 kappa = rd/cpd
 etv = rv/rd-1.
+cpp = cpd*1.
 #
 #     ----------------------------------------------------------------
@@ -153 +158 @@
 type_unsddu=3
 type_unsddv=4
+
+# Zoom related
+#
+fxyhypb = False
+nitergdiv = 1
+nitergrot = 2
+niterh = 2
+coefdis = 0.
+ysinus = True