Changeset 212 in lmdz_wrf for trunk

Timestamp:

Jan 9, 2015, 3:13:49 PM (10 years ago)

Author:

lfita

Message:

First final version without checking if everithing is done fine

Location:

trunk/tools

Files:

• iniaqua.py (modified) (10 diffs)
• lmdz_const.py (modified) (2 diffs)

trunk/tools/iniaqua.py

@@ -19 +19 @@
 filekinds = ['CF', 'WRF']
 
-## g.e. # iniaqua.py -d 32,32,39 -o WRF -t 19791201000000 -y true
+## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z param 
+
+# from dyn3d/fxy_new.h
+#c
+#c....stretching in x...
+#c
+def fripx(ri,dx):
+    fname = 'fripx'
+
+    ripx = (ri-1.0) *2.*np.pi/np.float(dx)
+
+    return ripx
+
+def ffx(ri):
+    fname = 'ffx'
+
+    fx = fripx(ri) + transx  + alphax * np.sin(fripx(ri)+transx-pxo ) - np.pi
+
+    return fx
+
+def ffxprim(ri,dx):
+    fname = 'ffxprim'
+
+    fxprim = 2.*np.pi/REAL(dx)*( 1.+ alphax*np.cos( fripx(ri)+transx-pxo ) )
+
+    return fxprim
+
+#c....stretching in y...
+#c
+def fbigy(rj,jjp1,jjm):
+    fname = 'fbigy'
+
+    bigy = 2.*(np.float(jjp1)-rj)*np.pi/jjm
+
+    return bigy
+
+def fy(rj):
+    fname = 'fy'
+
+    fy = ( fbigy(rj) + transy + alphay * np.sin(fbigy(rj)+transy-pyo) )/2. - np.pi/2.
+
+    return fy
+
+def ffyprim(rj,jjm):
+    fname = 'ffyprim'
+
+    fyprim = ( np.pi/jjm ) * ( 1. + alphay * np.cos( fbigy(rj)+transy-pyo ) )
+
+    return fyprim
+
+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
+
+def jacobi(A,N,NP):
+    """!
+       ! $Id: jacobi.F90 1289 2009-12-18 14:51:22Z emillour $
+       !
+    """
+    fname = 'jacobi'
+
+    D = np.zeros((NP), dtype=np.float)
+    V = np.zeros((NP,NP), dtype=np.float)
+    B = np.zeros((NP), dtype=np.float)
+    Z = np.zeros((NP), dtype=np.float)
+
+    for IP in range(N):
+        V[IP,IP]=1.
+        B[IP] = A[IP,IP]
+
+    D = B
+    NROT = 0
+    for I in range(50):
+#      DO I=1,50 ! 50? I suspect this should be NP
+#                !     but convergence is fast enough anyway
+        SM = 0.
+        for IP in range(N-1):
+            for IQ in range(IP+1,N):
+                SM = SM + np.abs(A[IP,IQ])
+
+        if SM  == 0.: return A, D, V, NROT
+
+        if I < 4:
+          TRESH = 0.2*SM/N**2
+        else:
+          TRESH = 0.
+
+        for IP in range(N-1):
+            for IQ in range(IP+1,N):
+                G = 100.*np.abs(A[IP,IQ])
+                if I > 4 and np.abs(D[I]) + G == np.abs(D[IP]) and                   \
+                  np.abs(D[IQ])+G == np.abs(D[IQ]):
+                    A[IP,IQ] = 0.
+                elif np.abs(A[IP,IQ]) > TRESH:
+                    H = D[IQ]-D[IP]
+                    if np.abs(H)+G == np.abs(H):
+                        T = A[IP,IQ]/H
+                    else:
+                        THETA = 0.5*H/A[IP,IQ]
+                        T = 1./(np.abs(THETA)+np.sqrt(1.+THETA**2))
+                        if THETA < 0.: T = -T
+
+                    C=1./SQRT(1+T**2)
+                    S=T*C
+                    TAU=S/(1.+C)
+                    H=T*A[IP,IQ]
+                    Z[IP]=Z[IP]-H
+                    Z[IQ]=Z[IQ]+H
+                    D[IP]=D[IP]-H
+                    D[IQ]=D[IQ]+H
+                    A[IP,IQ]=0.
+                    for J in range(IP-1):
+                        G=A[J,IP]
+                        H=A[J,IQ]
+                        A[J,IP]=G-S*(H+G*TAU)
+                        A[J,IQ]=H+S*(G-H*TAU)
+
+                    for J in range(IP+1,IQ-1):
+                        G = A[IP,J]
+                        H = A[J,IQ]
+                        A[IP,J] = G-S*(H+G*TAU)
+                        A[J,IQ] = H+S*(G-H*TAU)
+
+                    for J in range(IQ+1,N):
+                        G = A[IP,J]
+                        H = A[IQ,J]
+                        A[IP,J] = G-S*(H+G*TAU)
+                        A[IQ,J] = H+S*(G-H*TAU)
+
+                    for J in range(N):
+                        G = V[J,IP]
+                        H = V[J,IQ]
+                        V[J,IP] = G-S*(H+G*TAU)
+                        V[J,IQ] = H+S*(G-H*TAU)
+
+                    NROT=NROT+1
+
+          B = B + Z
+          D = B
+          Z = 0.
+
+    print errormsg
+    print '  ' + fname + ': 50 iterations should never happen !!'
+    quit(-1)
+
+    return
+
+def acc(vec,im):
+    """!
+       ! $Header$
+       !
+    """
+    fname ='acc'
+
+    d = np.zeros((im), dtype=np.float)
+
+    for j in range(im):
+        for i in range(im):
+            d[i] = vec[i,j]*vec[i,j]
+
+        vsum = np.sqrt(np.sum(d))
+        for i in range(im):
+            vec[i,j] = vec[i,j] / vsum
+
+    return vec, d
+
+def eigen_sort(d,n,np):
+    """!
+       ! $Header$
+       !
+    """
+    fname = 'eigen_sort'
+
+    v = np.zeros((np,np), dtype=np.float)
+
+    for i in range(n-1):     
+        k=i
+        p=d[i]
+        for j in range(i+1,n):
+            if d[j] >= p:
+                k=j
+                p=d[j]
+          
+        if k != i:
+            d[k]=d[i]
+            d[i]=p
+            for j in range(n):
+                p=v[j,i]
+                v[j,i]=v[j,k]
+                v[j,k]=p
+
+    return d, v
+
+def inifgn(dx, dy, dv):
+    """! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004-05-19 12:53:09 lmdzadmin Exp $
+       !
+       c  
+       c    ...  H.Upadyaya , O.Sharma  ... 
+       c
+    """
+    fname = 'inifgn'
+
+    imm1  = dx - 1
+
+    rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, xprimu,   \
+      rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = fxy(dx,dy)
+
+    dlonu = xprimu
+    dlonv = xprimv
+
+    sddv = np.sqrt(dlonv)
+    sddu = np.sqrt(dlonu)
+    unsddu = 1./sddu
+    unsddv = 1./sddv
+
+    vec = np.zeros((dx,dy), dtype=np.float)
+    vec1 = np.zeros((dx,dy), dtype=np.float)
+    eigenfnv = np.zeros((dx,dy), dtype=np.float)
+    eigenfnu = np.zeros((dx,dy), dtype=np.float)
+c
+c
+    eignfnv[0,0] = -1.
+    eignfnv[dx-1,0] = 1.
+
+    for i in range(imm1):
+        eignfnv[i+1,i+1]= -1.
+        eignfnv[i,i+1] = 1.
+
+    for j in range(dx):
+        for i in range(dx):
+            eignfnv[i,j] = eignfnv[i,j]/(sddu[i]*sddv[j])
+
+    for j in range(dx):
+        for i in range(dx):
+            eignfnu[i,j] = -eignfnv[j,i]
+
+    for j in range(dx):
+        for i in range(dx):
+            vec[i,j] = 0.0
+            vec1[i,j] = 0.0
+            for k in range(dx):
+                vec[i,j] = vec[i,j]  + eignfnu[i,k]*eignfnv[k,j]
+                vec1[i,j] = vec1[i,j] + eignfnv[i,k]*eignfnu[k,j]
+
+    vec, dv, eignfnv, nrot = jacobi(vec,iim,iim)
+    eignfnv, d = acc(eignfnv,iim)
+    dv, eignfnv = eigen_sort(dv,iim,iim)
+c
+    vec1, du, eignfnu, nrot = jacobi(vec1,iim,iim)
+    eignfnu, d = acc(eignfnu,iim)
+    du, eignfnu = eigen_sort(du,iim,iim)
+
+    return sddu, unsddu, sddv, unsddv
+
+
+def inifilr(iim, xprimu):
+    """ from: filtrez/filtreg_mod.F90
+        !    ... H. Upadhyaya, O.Sharma   ...
+        !
+        !     version 3 .....
+
+        !     Correction  le 28/10/97    P. Le Van .
+        !  -------------------------------------------------------------------
+        !
+        ! ------------------------------------------------------------
+        !   This routine computes the eigenfunctions of the laplacien
+        !   on the stretched grid, and the filtering coefficients
+        !      
+        !  We designate:
+        !   eignfn   eigenfunctions of the discrete laplacien
+        !   eigenvl  eigenvalues
+        !   jfiltn   indexof the last scalar line filtered in NH
+        !   jfilts   index of the first line filtered in SH
+        !   modfrst  index of the mode from WHERE modes are filtered
+        !   modemax  maximum number of modes ( im )
+        !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
+        !   sdd      SQRT( dx )
+        !      
+        !     the modes are filtered from modfrst to modemax
+        !      
+        !-----------------------------------------------------------
+        !
+    """
+    fname = 'inifilr'
+
+    dlonu = xprimu
+
+    !
+    CALL inifgn(eignvl)
+    !
+    PRINT *,'inifilr: EIGNVL '
+    PRINT 250,eignvl
+250 FORMAT( 1x,5e14.6)
+    !
+    ! compute eigenvalues and eigenfunctions
+    !
+    !
+    !.................................................................
+    !
+    !  compute the filtering coefficients for scalar lines and 
+    !  meridional wind v-lines
+    !
+    !  we filter all those latitude lines WHERE coefil < 1
+    !  NO FILTERING AT POLES
+    !
+    !  colat0 is to be used  when alpha (stretching coefficient)
+    !  is set equal to zero for the regular grid CASE 
+    !
+    !    .......   Calcul  de  colat0   .........
+    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
+    !
+    !
+    DO j = 1,jjm
+       dlatu( j ) = rlatu( j ) - rlatu( j+1 )
+    ENDDO
+    !
+#ifdef CRAY
+    iymin   = ISMIN( jjm, dlatu, 1 )
+    ixmineq = ISMIN( iim, dlonu, 1 )
+    dymin   = dlatu( iymin )
+    dxmin   = dlonu( ixmineq )
+#else
+    dxmin   =  dlonu(1)
+    DO  i  = 2, iim
+       dxmin = MIN( dxmin,dlonu(i) )
+    ENDDO
+    dymin  = dlatu(1)
+    DO j  = 2, jjm
+       dymin = MIN( dymin,dlatu(j) )
+    ENDDO
+#endif
+    !
+    ! For a regular grid, we want the filter to start at latitudes
+    ! corresponding to lengths dx of the same size as dy (in terms
+    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
+    !  <=> latitude=60 degrees).
+    ! Same idea for the zoomed grid: start filtering polewards as soon
+    ! as length dx becomes of the same size as dy 
+    !
+    colat0  =  MIN( 0.5, dymin/dxmin )
+    !
+    IF( .NOT.fxyhypb.AND.ysinus )  THEN
+       colat0 = 0.6
+       !         ...... a revoir  pour  ysinus !   .......
+       alphax = 0.
+    ENDIF
+    !
+    PRINT 50, colat0,alphax
+50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
+    !
+    IF(alphax.EQ.1. )  THEN
+       PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
+       STOP
+    ENDIF
+    !
+    lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
+
+    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
+    !
+    DO i = 2,iim
+       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
+    ENDDO
+    !
+
+    DO j = 1,jjm
+       DO i = 1,iim
+          coefilu( i,j )  = 0.0
+          coefilv( i,j )  = 0.0
+          coefilu2( i,j ) = 0.0
+          coefilv2( i,j ) = 0.0
+       ENDDO
+    ENDDO
+
+    !
+    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
+    !    .........................................................
+    !
+    modemax = iim
+
+!!!!    imx = modemax - 4 * (modemax/iim)
+
+    imx  = iim
+    !
+    PRINT *,'inifilr: TRUNCATION AT ',imx
+    !
+! Ehouarn: set up some defaults
+    jfiltnu=2 ! avoid north pole
+    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
+    jfiltnv=1 ! NB: no poles on the V grid
+    jfiltsv=jjm
+
+    DO j = 2, jjm/2+1
+       cof = COS( rlatu(j) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
+            jfiltnu= j
+          ENDIF
+       ENDIF
+
+       cof = COS( rlatu(jjp1-j+1) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
+               jfiltsu= jjp1-j+1
+          ENDIF
+       ENDIF
+    ENDDO
+    !
+    DO j = 1, jjm/2
+       cof = COS( rlatv(j) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
+            jfiltnv= j
+          ENDIF
+       ENDIF
+
+       cof = COS( rlatv(jjm-j+1) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
+               jfiltsv= jjm-j+1
+          ENDIF
+       ENDIF
+    ENDDO
+    !                                 
+
+    IF( jfiltnu.GT. jjm/2 +1 )  THEN
+       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
+       STOP
+    ENDIF
+
+    IF( jfiltsu.GT.  jjm  +1 )  THEN
+       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
+       STOP
+    ENDIF
+
+    IF( jfiltnv.GT. jjm/2    )  THEN
+       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
+       STOP
+    ENDIF
+
+    IF( jfiltsv.GT.     jjm  )  THEN
+       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
+       STOP
+    ENDIF
+
+    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
+         jfiltnv,jfiltsv,jfiltnu,jfiltsu
+
+    IF(first_call_inifilr) THEN
+       ALLOCATE(matriceun(iim,iim,jfiltnu))
+       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
+       ALLOCATE(matricevn(iim,iim,jfiltnv))
+       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
+       ALLOCATE( matrinvn(iim,iim,jfiltnu))
+       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
+       first_call_inifilr = .FALSE.
+    ENDIF
+
+    !                                 
+    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
+    !................................................................
+    !
+    !
+    DO j = 1,jjm
+    !default initialization: all modes are retained (i.e. no filtering)
+       modfrstu( j ) = iim
+       modfrstv( j ) = iim
+    ENDDO
+    !
+    DO j = 2,jfiltnu
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          IF ( cof .LT. 1. ) GOTO 82
+       ENDDO
+       GOTO 84
+82     modfrstu( j ) = k
+       !
+       kf = modfrstu( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          coefilu(k,j) = cof - 1.
+          coefilu2(k,j) = cof*cof - 1.
+       ENDDO
+84     CONTINUE
+    ENDDO
+    !                                 
+    !
+    DO j = 1,jfiltnv
+       !
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          IF ( cof .LT. 1. ) GOTO 87
+       ENDDO
+       GOTO 89
+87     modfrstv( j ) = k
+       !
+       kf = modfrstv( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          coefilv(k,j) = cof - 1.
+          coefilv2(k,j) = cof*cof - 1.
+       ENDDO
+89     CONTINUE
+    ENDDO
+    !
+    DO j = jfiltsu,jjm
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          IF ( cof .LT. 1. ) GOTO 92
+       ENDDO
+       GOTO 94
+92     modfrstu( j ) = k
+       !
+       kf = modfrstu( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          coefilu(k,j) = cof - 1.
+          coefilu2(k,j) = cof*cof - 1.
+       ENDDO
+94     CONTINUE
+    ENDDO
+    !                                 
+    DO j = jfiltsv,jjm
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          IF ( cof .LT. 1. ) GOTO 97
+       ENDDO
+       GOTO 99
+97     modfrstv( j ) = k
+       !
+       kf = modfrstv( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          coefilv(k,j) = cof - 1.
+          coefilv2(k,j) = cof*cof - 1.
+       ENDDO
+99     CONTINUE
+    ENDDO
+    !
+
+    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
+! Ehouarn: and what are these for??? Trying to handle a limit case
+!          where filters extend to and meet at the equator?
+       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
+       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
+
+       PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
+            jfiltnv,jfiltsv,jfiltnu,jfiltsu
+    ENDIF
+
+    PRINT *,'   Modes premiers  v  '
+    PRINT 334,modfrstv
+    PRINT *,'   Modes premiers  u  '
+    PRINT 334,modfrstu
+
+    !  
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
+    !                       sur la grille scalaire                 ........
+    !   ...................................................................
+    !
+    DO j = 2, jfiltnu
+
+       DO i=1,iim
+          coff = coefilu(i,j)
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO ! of DO i=1,iim
+#ifdef CRAY
+       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matriceun(i,k,j) = 0.0
+             DO ii = 1, iim
+                matriceun(i,k,j) = matriceun(i,k,j) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO ! of DO k = 1, iim
+#endif
+#endif
+
+    ENDDO ! of DO j = 2, jfiltnu
+
+    DO j = jfiltsu, jjm
+
+       DO i=1,iim
+          coff = coefilu(i,j)
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO ! of DO i=1,iim
+#ifdef CRAY
+       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, &
+            matriceus(1,1,j-jfiltsu+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matriceus(i,k,j-jfiltsu+1) = 0.0
+             DO ii = 1, iim
+                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO ! of DO k = 1, iim
+#endif
+#endif
+
+    ENDDO ! of DO j = jfiltsu, jjm
+
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
+    !                       sur la grille   de V ou de Z           ........
+    !   ...................................................................
+    !
+    DO j = 1, jfiltnv
+
+       DO i = 1, iim
+          coff = coefilv(i,j)
+          IF( i.LT.modfrstv(j) ) coff = 0.
+          DO k = 1, iim
+             eignft(i,k) = eignfnu(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matricevn(i,k,j) = 0.0
+             DO ii = 1, iim
+                matricevn(i,k,j) = matricevn(i,k,j) &
+                     + eignfnu(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO ! of DO j = 1, jfiltnv
+
+    DO j = jfiltsv, jjm
+
+       DO i = 1, iim
+          coff = coefilv(i,j)
+          IF( i.LT.modfrstv(j) ) coff = 0.
+          DO k = 1, iim
+             eignft(i,k) = eignfnu(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnu, iim, eignft, iim, 0.0, & 
+            matricevs(1,1,j-jfiltsv+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matricevs(i,k,j-jfiltsv+1) = 0.0
+             DO ii = 1, iim
+                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
+                     + eignfnu(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO ! of DO j = jfiltsv, jjm
+
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
+    !              sur la grille scalaire , pour le filtre inverse ........
+    !   ...................................................................
+    !
+    DO j = 2, jfiltnu
+
+       DO i = 1,iim
+          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matrinvn(i,k,j) = 0.0
+             DO ii = 1, iim
+                matrinvn(i,k,j) = matrinvn(i,k,j) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO ! of DO j = 2, jfiltnu
+
+    DO j = jfiltsu, jjm
+
+       DO i = 1,iim
+          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matrinvs(i,k,j-jfiltsu+1) = 0.0
+             DO ii = 1, iim
+                matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO ! of DO j = jfiltsu, jjm
+
+    IF (use_filtre_fft) THEN
+       CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
+                           coefilv,modfrstv,jfiltnv,jfiltsv)
+       CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
+                           coefilv,modfrstv,jfiltnv,jfiltsv)
+    ENDIF
+
+    !   ...................................................................
+
+    !
+334 FORMAT(1x,24i3)
+755 FORMAT(1x,6f10.3,i3)
+
+    RETURN
+  END SUBROUTINE inifilr
+
+
+def filtreg(dx, dy, champorig, nlat, nbniv, ifiltre, iaire, griscal ,iterv):
+    """c=======================================================================
+       c
+       c   Auteur: P. Le Van        07/10/97
+       c   ------
+       c
+       c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
+       c                     pour l'operateur  Filtre    .
+       c   ------
+       c
+       c   Arguments:
+       c   ----------
+       c
+       c      nblat                 nombre de latitudes a filtrer
+       c      nbniv                 nombre de niveaux verticaux a filtrer
+       c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
+       c                            en sortie : champ filtre
+       c      ifiltre               +1  Transformee directe
+       c                            -1  Transformee inverse
+       c                            +2  Filtre directe
+       c                            -2  Filtre inverse
+       c
+       c      iaire                 1   si champ intensif
+       c                            2   si champ extensif (pondere par les aires)
+       c
+       c      iterv                 1   filtre simple
+       c
+       c=======================================================================
+       c
+       c
+       c                      Variable Intensive
+       c                ifiltre = 1     filtre directe
+       c                ifiltre =-1     filtre inverse
+       c
+       c                      Variable Extensive
+       c                ifiltre = 2     filtre directe
+       c                ifiltre =-2     filtre inverse
+       c
+    """
+
+    fname = 'filtreg'
+    first = True
+
+    sdd12 = np.zeros((dx,4), dtype = np.float)
+
+    champ = fortran_mat(champorig)
+
+    sddu, unsddu, sddv, unsddv = inifgn(dx, dy)
+
+    if first:
+         sdd12[:,type_sddu] = sddu
+         sdd12[:,type_sddv] = sddv
+         sdd12[:,type_unsddu] = unsddu
+         sdd12[:,type_unsddv] = unsddv
+
+         first=False
+
+    if ifiltre  == 1 or ifiltre == -1:
+        print '  ' + fname + ': Pas de transformee simple dans cette version'
+        quit()
+      
+    if iter == 2:
+        print errormsg
+        print '  ' + fname + ': Pas d iteration du filtre dans cette version !'
+        print ' Utiliser old_filtreg et repasser !'
+        quit()
+      
+    if ifiltre == -2 and not griscal:
+        print errormsg
+        print '  ' + fname + ' Cette routine ne calcule le filtre inverse que ' +    \
+          ' sur la grille des scalaires !'
+        quit(-1)
+      
+    if ifiltre != 2 and ifiltre != - 2:
+        print errormsg
+        print '  ' + fname + ': Probleme dans filtreg car ifiltre NE 2 et NE -2' +   \
+          ' corriger et repasser !'
+        quit(-1)
+      
+    iim = dx
+    jjm = dy
+
+    iim2 = iim*iim
+    immjm = iim*jjm
+
+    if griscal:
+        if nlat != jjp1:
+            print errormsg
+            print '  ' + fname  + ': 1111'
+            quit(-1)
+        else
+            if iaire == 1:
+                sdd1_type = type_sddv
+                sdd2_type = type_unsddv
+            else:
+                sdd1_type = type_unsddv
+                sdd2_type = type_sddv
+
+            jdfil1 = 2
+            jffil1 = jfiltnu
+            jdfil2 = jfiltsu
+            jffil2 = jjm
+    else:
+        if nlat != jjm:
+            print errormsg
+            print '  ' + fname  + ': 2222'
+            quit(-1)
+        else:           
+            if iaire == 1:
+                sdd1_type = type_sddu
+                sdd2_type = type_unsddu
+            else:
+                sdd1_type = type_unsddu
+                sdd2_type = type_sddu
+
+            
+            jdfil1 = 1
+            jffil1 = jfiltnv
+            jdfil2 = jfiltsv
+            jffil2 = jjm
+      
+    for hemisph in range(1,3):
+        if (hemisph == 1 ):
+            jdfil = jdfil1
+            jffil = jffil1
+        else:
+            jdfil = jdfil2
+            jffil = jffil2
+ 
+        for l in range(nbniv):
+            for j in rangE(jdfil,jffil):
+                for i in range(iim):
+                    champ[i,j,l] = champ[i,j,l]*sdd12[i,sdd1_type] # sdd1(i)
+         
+        if hemisph == 1:
+            if ifiltre == -2:
+               for j in range(jdfil,jffil):
+                  eignq[:,j-jdfil+1,:] = np.matmul(matrinvn[:,:,j], champ[:,j,:])
+#endif
+               END DO
+               
+            ELSE IF ( griscal )     THEN
+               
+               DO j = jdfil,jffil
+#ifdef BLAS
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matriceun(1,1,j),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+#else
+                  eignq(:,j-jdfil+1,:)
+     $                 = matmul(matriceun(:,:,j), champ(:iim,j,:))
+#endif
+               END DO
+               
+            ELSE 
+               
+               DO j = jdfil,jffil
+#ifdef BLAS
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matricevn(1,1,j),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+#else
+                  eignq(:,j-jdfil+1,:)
+     $                 = matmul(matricevn(:,:,j), champ(:iim,j,:))
+#endif
+               END DO
+               
+            ENDIF
+            
+         ELSE
+            
+            IF( ifiltre. EQ. -2 )   THEN
+               
+               DO j = jdfil,jffil
+#ifdef BLAS
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matrinvs(1,1,j-jfiltsu+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+#else
+                  eignq(:,j-jdfil+1,:)
+     $                 = matmul(matrinvs(:,:,j-jfiltsu+1),
+     $                 champ(:iim,j,:))
+#endif
+               END DO
+               
+               
+            ELSE IF ( griscal )     THEN
+               
+               DO j = jdfil,jffil
+#ifdef BLAS
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matriceus(1,1,j-jfiltsu+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+#else
+                  eignq(:,j-jdfil+1,:)
+     $                 = matmul(matriceus(:,:,j-jfiltsu+1),
+     $                 champ(:iim,j,:))
+#endif
+               END DO
+                              
+            ELSE 
+               
+               DO j = jdfil,jffil
+#ifdef BLAS
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matricevs(1,1,j-jfiltsv+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+#else
+                  eignq(:,j-jdfil+1,:)
+     $                 = matmul(matricevs(:,:,j-jfiltsv+1),
+     $                 champ(:iim,j,:))
+#endif
+               END DO
+                              
+            ENDIF
+            
+         ENDIF
+         
+         IF( ifiltre.EQ. 2 )  THEN
+            
+            DO l = 1, nbniv
+               DO j = jdfil,jffil
+                  DO i = 1, iim
+                     champ( i,j,l ) = 
+     &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
+     &                    * sdd12(i,sdd2_type) ! sdd2(i)
+                  END DO
+               END DO
+            END DO
+
+         ELSE
+
+            DO l = 1, nbniv
+               DO j = jdfil,jffil
+                  DO i = 1, iim
+                     champ( i,j,l ) = 
+     &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
+     &                    * sdd12(i,sdd2_type) ! sdd2(i)
+                  END DO
+               END DO
+            END DO
+
+         ENDIF
+
+         DO l = 1, nbniv
+            DO j = jdfil,jffil
+               champ( iip1,j,l ) = champ( 1,j,l )
+            END DO
+         END DO
+
+     
+      ENDDO
+
+1111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a 
+     &     filtrer, sur la grille des scalaires'/)
+2222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
+     &     ltrer, sur la grille de V ou de Z'/)
+      RETURN
+      END
+
+def SSUM(n,sx,incx):
+    """ Obsolete version of sum for non Fortan 90 code
+      from dyn3d/cray.F
+    """
+
+    ssumv = 0.
+
+    ix = 0
+
+    for i in range(n)
+        ssumv=ssumv+sx[ix]
+        ix=ix+incx
+
+    return ssumv
+
+def SCOPY(n,sx,incx,sy,incy):
+    """ Obsolete function to copy matrix values
+      from dyn3d/cray.F
+    """
+    iy = 1
+    ix = 1
+
+    for i in range(n):
+        sy[iy] = sx[ix]
+        ix = ix+incx
+        iy = iy+incy
+
+    return sy
+
+def exner_hyb (dx, dy, dz, psv, pv, alphav, betav)
+    """c
+       c     Auteurs :  P.Le Van  , Fr. Hourdin  .
+       c    ..........
+       c
+       c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+       c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+       c
+       c   ************************************************************************
+       c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
+       c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+       c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+       c   ************************************************************************
+       c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+       c    la pression et la fonction d'Exner  au  sol  .
+       c
+       c                                 -------- z                                   
+       c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+       c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+       c    ( voir note de Fr.Hourdin )  ,
+       c
+       c    on determine successivement , du haut vers le bas des couches, les 
+       c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+       c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+       c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+       c
+    """
+
+    fname = 'exner_hyb'
+
+    pks = np.zeros((dy, dx), dtype=np.float)
+    pk = np.zeros((dz, dy, dx), dtype=np.float)
+    pkf = np.zeros((dz, dy, dx), dtype=np.float)
+
+    if dz == 1:       
+# Compute pks(:),pk(:),pkf(:)
+        pks = (cpp/preff)*ps
+        pk[0,:,:] = 0.5*pks
+#        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 ) is the same as the next line?
+        pkf = pk
+        
+#  No filtering... not necessary on aquaplanet
+#        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+# our work is done, exit routine
+        return pksv, pkv, pkfv
+
+#!!!! General case:
+     
+    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]
+
+    xpn = SSUM(iim,ppn,1) /apoln
+    xps = SSUM(iim,pps,1) /apols
+
+    for ij in range(iip1):
+        pks[ij] = xpn
+        pks[ij+ip1jm] = xps
+#
+#
+#    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+#
+    for ij in range(ngrid):
+        alpha[ij,llm] = 0.
+        beta [ij,llm] = 1./ unpl2k
+
+#
+#     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+#
+    for l in range (llm-1,1,-1):
+c
+        for ij in range(ngrid):
+            dellta = p[ij,l]* unpl2k + p[ij,l+1]* ( beta[ij,l+1]-unpl2k )
+            alpha[ij,l] = -p[ij,l+1] / dellta * alpha[ij,l+1]
+            beta[ij,l] = p[ij,l] / dellta   
+
+#  ***********************************************************************
+#     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+#
+    for ij in range(ngrid):
+        pk[ij,0] = ( p[ij,0]*pks[ij] - 0.5*alpha[ij,1]*p[ij,1] )                     \
+          *(  p[ij,0]* (1.+kappa) + 0.5*( beta[ij,1]-unpl2k )* p[ij,1] )
+
+#
+#    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+#
+    for l in range(1,llm):
+        for ij in range(ngrid):
+            pk[ij,l] = alpha[ij,l] + beta[ij,l] * pk[ij,l-1]
+#
+#
+    pkfv = SCOPY ( ngrid * llm, pk, 1, pkfv, 1 )
+
+# We do not filter for iniaqua
+#      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+
+    return pksv, pkv, pkfv
+
+def exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ):
+    """c
+       c     Auteurs :  F. Forget , Y. Wanherdrick
+       c P.Le Van  , Fr. Hourdin  .
+       c    ..........
+       c
+       c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+       c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+       c
+       c   ************************************************************************
+       c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+       c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+       c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+       c   ************************************************************************
+       c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+       c    la pression et la fonction d'Exner  au  sol  .
+       c
+       c     WARNING : CECI est une version speciale de exner_hyb originale
+       c               Utilise dans la version martienne pour pouvoir 
+       c               tourner avec des coordonnees verticales complexe
+       c              => Il ne verifie PAS la condition la proportionalite en 
+       c              energie totale/ interne / potentielle (F.Forget 2001)
+       c    ( voir note de Fr.Hourdin )  ,
+       c
+    """
+    fname = 'exner_milieu'
+
+    pks = np.zeros((dy, dx), dtype=np.float)
+    pk = np.zeros((dz, dy, dx), dtype=np.float)
+    pkf = np.zeros((dz, dy, dx), dtype=np.float)
+
+    ppn = np.zeros((iim), dtype=np.float))
+    ppn = np.zeros((iim), dtype=np.float))
+
+    firstcall = True
+    modname = 'exner_milieu'
+
+# Sanity check
+    if firstcall:
+#       sanity checks for Shallow Water case (1 vertical layer)
+        if llm == 1:
+            if kappa != 1:
+                print errormsg
+                print '  ' + fname+ ': kappa!=1 , but running in Shallow Water mode!!'
+                quit(-1)
+            if cpp != r:
+                print errormsg
+                print '  ' + fname+ ': cpp!=r , but running in Shallow Water mode!!'
+                quit(-1)
+
+
+        firstcall = False
+
+
+#### Specific behaviour for Shallow Water (1 vertical layer) case:
+      if llm == 1:
+      
+#         Compute pks(:),pk(:),pkf(:)
+        
+          for ij in range(ngrid):
+              pks[ij] = (cpp/preff) * ps[ij]
+              pk[ij,1] = .5*pks[ij]
+
+          
+          pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
+# We do not filter for iniaqua
+#        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+#         our work is done, exit routine
+          return pksv, pkv, pkfv
+
+
+#### General case:
+
+#     -------------
+#     Calcul de pks
+#     -------------
+   
+    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]
+
+    xpn = SSUM(iim,ppn,1) /apoln
+    xps = SSUM(iim,pps,1) /apols
+
+    for ij in range(iip1)
+        pks[ij] = xpn
+        pks[ij+ip1jm] = xps
+
+#
+#
+#    .... Calcul de pk  pour la couche l 
+#    --------------------------------------------
+#
+    dum1 = cpp * (2*preff)**(-kappa) 
+    for l in range(llm-1):
+        for ij in range(ngrid):
+            pk[ij,l] = dum1 * (p[ij,l] + p[ij,l+1])**kappa
+
+#    .... Calcul de pk  pour la couche l = llm ..
+#    (on met la meme distance (en log pression)  entre Pk(llm)
+#    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+    for ij in range(ngrid):
+        pk[ij,llm] = pk[ij,llm-1]**2 / pk[ij,llm-2]
+
+#    calcul de pkf
+#    -------------
+    pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
+
+# We do not filter iniaqua
+#      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+      
+#    EST-CE UTILE ?? : calcul de beta
+#    --------------------------------
+    for l in range(1, llm):
+        for ij in range(ngrid):
+            beta[ij,l] = pk[ij,l] / pk[ij,l-1]
+
+    return pksv, pkv, pkfv
+
+def pression(dx, dy, dz, apv, bpv, psv)
+    """c
+
+       c      Auteurs : P. Le Van , Fr.Hourdin  .
+
+      c  ************************************************************************
+      c     Calcule la pression p(l) aux differents niveaux l = 1 ( niveau du
+      c     sol) a l = llm +1 ,ces niveaux correspondant aux interfaces des (llm) 
+      c     couches , avec  p(ij,llm +1) = 0.  et p(ij,1) = ps(ij)  .      
+      c  ************************************************************************
+      c
+    """ 
+    fname = 'pression'
+      
+    press = np.zeros((dz, dy, dx), dtype=np.float)
+
+    for l in range(dz+1):
+        press[l,:,:] = apv[l] + bpv[l]*psv[:,:]
+   
+    return press
 
 def sig_hybrid(sig,pa,preff):
@@ -76 +1449 @@
     return nsig
 
-def presnivs_calc(hybrid, dz):
+def presnivs_calc(vert_sampling, dz):
+    """ From dyn3d/disvert.F calculation of vertical pressure levels
+      vert_sampling= which kind of vertical sampling is desired: "param", "tropo", 
+        "strato" and "read"
+      dz= numbef of vertical layers
+    """
+
+    fname = 'presnivs_calc'
+
+    llmp1 = dz + 1
+    pnivs = np.zeros((dz), dtype=np.float)
+    s = np.zeros((dz), dtype=np.float)
+    sig = np.zeros((dz+1), dtype=np.float)
+    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)
+
+# default scaleheight is 8km for earth
+    scaleheight = 8.
+
+#    vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
+
+    if vert_sampling == 'param':
+# On lit les options dans sigma.def:
+        if not os.path.isfile('easig.def'):
+            print errormsg
+            print '  ' + fname  + ": parameters file 'easig.def' does not exist!!"
+            quit(-1)
+
+        sfobj = open('sigma.def', 'r')
+        scaleheight = np.float( ncvar.reduce_spaces(fobj.readline()))
+        deltaz = np.float( ncvar.reduce_spaces(fobj.readline()))
+        beta = np.float( ncvar.reduce_spaces(fobj.readline()))
+        k0 = np.float( ncvar.reduce_spaces(fobj.readline()))
+        k1 = np.float( ncvar.reduce_spaces(fobj.readline()))
+        sfobj.close()
+
+        alpha = deltaz/(dz*scaleheight)
+        print ':scaleheight, alpha, k0, k1, beta', scaleheight, alpha, k0, k1, beta
+
+        alpha=deltaz/np.tanh(1./k0)*2.
+        zkm1=0.
+        sig[0]=1.
+        for l in range(dz):
+            sig[l+1]=(np.cosh(l/k0))**(-alpha*k0/scaleheight)                        \
+               *exp(-alpha/scaleheight*np.tanh((llm-k1)/k0)                          \
+               *beta**(l-(llm-k1))/np.log(beta))
+
+            zk=-scaleheight*np.log(sig[l+1])
+
+            dzk1=alpha*np.tanh(l/k0)
+            dzk2=alpha*np.tanh((llm-k1)/k0)*beta**(l-(llm-k1))/np.log(beta)
+
+            print l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
+            zkm1=zk
+
+        sig[dz-1]=0.
+
+        bp[0:dz] = np.exp(1.-1./sig[0:dz]**2)
+        bp[llmp1-1] = 0.
+
+        ap = pa * (sig - bp)
+
+    elif vert_sampling  == 'tropo':
+        for l in range(dz):
+            x = 2.*np.arcsin(1.)*(l-0.5)/(dz+1.)
+            dsig[l] = 1.0+7.0*np.sin(x)**2
+
+        dsig = dsig/np.sum(dsig)
+        sig[dz] = 0.
+        for l in range(dz-1,0,-1):
+            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])
+
+    elif vert_sampling  == 'strato':
+        if dz == 39:
+            dsigmin = 0.3
+        elif dz == 50:
+            dsigmin = 1.
+        else:
+            print ' ATTENTION discretisation z a ajuster'
+            dsigmin = 1.
+
+        print 'Discretisation verticale DSIGMIN=',dsigmin
+
+        for l in range(dz):
+            x = 2.*np.arcsin(1.)*(l - 0.5)/(dz+1)
+            dsig[l] =(dsigmin+7.*np.sin(x)**2)                                       \
+             *(0.5*(1.-np.tanh(1.*(x-np.arcsin(1.))/np.arcsin(1.))))**2
+
+        dsig = dsig/np.sum(dsig)
+        sig[dz] = 0.
+        for l in range(dz-1,0,-1):
+            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])
+
+    elif vert_sampling  == 'read':
+# Read "ap" and "bp". First line is skipped (title line). "ap"
+# should be in Pa. First couple of values should correspond to
+# the surface, that is : "bp" should be in descending order.
+        if not os.path.isfile('hybrid.txt'):
+            print errormsg
+            print '  ' + fname  + ": parameters file 'hybrid.txt' does not exist!!"
+            quit(-1)
+        sfobj = ('hybrid.txt', 'r')
+# skip title line
+        title = sfobj.readline()
+        for l in range(dz+1):
+            values = ncvar.readuce_space(sfobj.readline())
+            ap[l] = np.float(values[0])
+            bp[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.:
+            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]
+          
+    else:
+        print errormsg
+        print '  ' + fname + ': wrong value for vert_sampling:', vert_sampling
+        quit(-1)
+
+
+    nivsigs = np.arange(dz)*1.+1.
+    nivsig = np.arange(llmp1)*1.+1.
+
+    print '  ' + fname + ': k BP AP ________'
+    for k in range(dz+1):
+        print k, bp[k], ap[k]
+
+    print 'Niveaux de pressions approximatifs aux centres des'
+    print 'couches calcules pour une pression de surface =', preff
+    print 'et altitudes equivalentes pour une hauteur d echelle de '
+    print scaleheight,' km'
+
+    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 )
+        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]))
+
+    print '  ' + fname + ': PRESNIVS [Pa]:', pnivs
+
+    return pnivs
+
+
+def presnivs_calc_noterre(hybrid, dz):
     """ From dyn3d/disvert_noterre.F calculation of vertical pressure levels
       hybrid= whether hydbrid coordinates have to be used
@@ -82 +1619 @@
     """
 
-    fname = 'presnivs_calc'
+    fname = 'presnivs_calc_noterre'
 
 #-----------------------------------------------------------------------
@@ -272 +1809 @@
     return longitude, latitude
 
+def massdair(p):
+    """c
+       c *********************************************************************
+       c       ....  Calcule la masse d'air  dans chaque maille   ....
+       c *********************************************************************
+       c
+       c    Auteurs : P. Le Van , Fr. Hourdin  .
+       c   ..........
+       c
+       c  ..    p                      est  un argum. d'entree pour le s-pg ...
+       c  ..  masse                    est un  argum.de sortie pour le s-pg ...
+       c     
+       c  ....  p est defini aux interfaces des llm couches   .....
+       c
+    """
+    fname = 'massdair'
+
+    masse = np.zeros((ip1jmp1,llm), dtype=np.float)
+#
+#
+#   Methode pour calculer massebx et masseby .
+#   ----------------------------------------
+#
+#    A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
+#       alpha1(i,j)  calcule  au point ( i+1/4,j-1/4 )
+#       alpha2(i,j)  calcule  au point ( i+1/4,j+1/4 )
+#       alpha3(i,j)  calcule  au point ( i-1/4,j+1/4 )
+#       alpha4(i,j)  calcule  au point ( i-1/4,j-1/4 )
+#
+#    Avec  alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)        
+#
+#    N.B .  Pour plus de details, voir s-pg  ...  iniconst ...
+#
+#
+#
+#   alpha4 .         . alpha1    . alpha4
+#    (i,j)             (i,j)       (i+1,j)
+#
+#             P .        U .          . P
+#           (i,j)       (i,j)         (i+1,j)
+#
+#   alpha3 .         . alpha2    .alpha3 
+#    (i,j)              (i,j)     (i+1,j)
+#
+#             V .        Z .          . V
+#           (i,j)
+#
+#   alpha4 .         . alpha1    .alpha4
+#   (i,j+1)            (i,j+1)   (i+1,j+1) 
+#
+#             P .        U .          . P
+#          (i,j+1)                    (i+1,j+1)
+#
+#
+#
+#                       On  a :
+#
+#    massebx(i,j) = masse(i  ,j) * ( alpha1(i  ,j) + alpha2(i,j))   +
+#                   masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
+#     localise  au point  ... U (i,j) ...
+#
+#    masseby(i,j) = masse(i,j  ) * ( alpha2(i,j  ) + alpha3(i,j  )  +
+#                   masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)  
+#     localise  au point  ... V (i,j) ...
+#
+#
+#=======================================================================
+
+    for l in range (llm):
+        for ij in range(ip1jmp1):
+            masse[ij,l] = airesurg[ij] * ( p[ij,l] - p[ij,l+1] )
+
+        for ij in range(ip1jmp1,iip1)
+            masse[ij+iim,l] = masse[ij,l]
+       
+    return masse
+
+def geopot(ngrid, teta, pk, pks)
+    """c=======================================================================
+       c
+       c   Auteur:  P. Le Van
+       c   -------
+       c
+       c   Objet:
+       c   ------
+       c
+       c    *******************************************************************
+       c    ....   calcul du geopotentiel aux milieux des couches    .....
+       c    *******************************************************************
+       c
+       c     ....   l'integration se fait de bas en haut  ....
+       c
+       c     .. ngrid,teta,pk,pks,phis sont des argum. d'entree pour le s-pg ..
+       c              phi               est un  argum. de sortie pour le s-pg .
+       c
+       c=======================================================================
+    """
+    fname = 'geopot'
+
+    phis = np.zeros((ngrid), dtype=np.float)
+    phi = np.zeros((ngrid,llm), dtype=np.float)
+
+#-----------------------------------------------------------------------
+#     calcul de phi au niveau 1 pres du sol  .....
+
+    for ij in range(ngrid):
+        phi[ij,1] = phis[ij] + teta[ij,0] * ( pks[ij] - pk[ij,0] )
+
+#     calcul de phi aux niveaux superieurs  .......
+
+    for l in range(1,llm):
+        for ij in range(ngrid):
+            phi[ij,l] = phi[ij,l-1] + 0.5 * ( teta[ij,l]+teta[ij,l-1] ) *            \
+              ( pk[ij,l-1]-pk[ij,l] )
+
+    return phis, phi 
+
+def dump2d(im,jm,nom_z):
+    """ Function to create a dump 2d variable
+      from dyn3d/dump2d.F
+    """
+    fname = 'dump2d'
+
+    z = np.zeros((im,jm), dtype=np.float)
+
+    zmin = z[0,0]
+    zllm = z[0,0]
+    imin = 0
+    illm = 0
+    jmin = 0
+    jllm = 0
+
+    for j in range(jm):
+        for i in range(im):
+            if z[i,j] > zllm:
+                illm=i
+                jllm=j
+                zllm=z[i,j]
+
+            if z[i,j] < zmin:
+                imin=i
+                jmin=j
+                zmin=z[i,j]
+
+    print 'MIN:',zmin
+    print 'MAX:',zllm
+
+    if zllm > zmin:
+        for j in range(jm):
+            print int(10.*(z[:,j]-zmin)/(zllm-zmin))
+
+    return z
+
+def ugeostr(phi):
+    """! Calcul du vent covariant geostrophique a partir du champ de
+       ! geopotentiel.
+       ! We actually compute: (1 - cos^8 \phi) u_g
+       ! to have a wind going smoothly to 0 at the equator.
+       ! We assume that the surface pressure is uniform so that model
+       ! levels are pressure levels.
+    """
+    fname = 'ugeostr'
+    ucov = np.zeros((iip1,jjp1,llm), dtype=np.float)
+    um = np.zeros((jjm,llm), dtype=np.float)
+    u = np.zeros((iip1,jjm,llm), dtype=np.float)
+
+
+    for j in range(jjm)
+
+        if np.abs(np.sin(rlatv[j])) < 1.e-4:
+            zlat = 1.e-4
+        else:
+            zlat=rlatv(j)
+
+        fact = np.cos(zlat)
+        fact = fact*fact
+        fact = fact*fact
+        fact = fact*fact
+        fact = (1.-fact)/ (2.*omeg*sin(zlat)*(rlatu[j+1]-rlatu[j]))
+        fact = -fact/rad
+        for l in range(llm):
+            for i in range(iim):
+                u[i,j,l] = fact*(phi[i,j+1,l]-phi[i,j,l])
+                um[j,l]=um[j,l]+u[i,j,l]/np.float(iim)
+
+    um = dump2d(jjm,llm,'Vent-u geostrophique')
+
+# calcul des champ de vent:
+
+    for l in range(llm):
+        for i in range(iip1):
+            ucov[i,1,l]=0.
+            ucov[i,jjp1,l]=0.
+        for j in range(1,jjm):
+            for i in range(iim):
+                ucov[i,j,l] = 0.5*(u[i,j,l]+u[i,j-1,l])*cu[i,j]
+
+        ucov[iip1,j,l]=ucov[0,j,l]
+
+    return ucov
+
+def RAN1(IDUM, Nvals):
+    """ Function to generate Nvals random numbers
+      from dyn3d/ran1.F
+      IDUM= Random Seed 
+      Nvals= number of values
+    """
+    fname = 'RAN1'
+    Nvals = 97
+
+    R = np.zeros((Nvals), dtype=np.float)
+
+    M1 = 259200
+    IA1 = 7141
+    IC1 = 54773 
+    RM1 = 3.8580247E-6
+    M2 = 134456 
+    IA2 = 8121
+    IC2 = 28411
+    RM2 = 7.4373773E-6
+    M3 = 243000 
+    IA3 = 4561
+    IC3 = 51349
+
+    if IDUM < 0 or IFF == 0:
+        IFF = 1
+        IX1 = np.mod(IC1-IDUM,M1)
+        IX1 = np.mod(IA1*IX1+IC1,M1)
+        IX2 = np.mod(IX1,M2)
+        IX1 = np.mod(IA1*IX1+IC1,M1)
+        IX3 = np.mod(IX1,M3)
+        for J in range(Nvals):
+            IX1 = np.mod(IA1*IX1+IC1,M1)
+            IX2 = np.mod(IA2*IX2+IC2,M2)
+            R[J] = (np.float(IX1)+np.float(IX2)*RM2)*RM1
+
+        IDUM=1
+
+    IX1 = np.mod(IA1*IX1+IC1,M1)
+    IX2 = np.mod(IA2*IX2+IC2,M2)
+    IX3 = np.mod(IA3*IX3+IC3,M3)
+    J = 1+(Nvals*IX3)/M3
+    if J > Nvals or J lt 1: quit()
+    ran1=R[J]
+    R[J]=(np.float(IX1)+np.float(IX2)*RM2)*RM1
+
+    return ran1
+
+
 ####### ###### ##### #### ### ## #
 
@@ -281 +2067 @@
 parser.add_option("-d", "--dimensions", dest="dims", 
   help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES")
+parser.add_option("-p", "--pressure_exner", dest="pexner, 
+  help="how as to b computed Exner pressure ('hybrid': hybrid coordinates, 'middle': middle layer)",  
+  metavar="VALUE")
 parser.add_option("-t", "--time", dest="time", 
   help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE")
-parser.add_option("-y", "--hybrid", dest="hybrid",
-  help="whether vertical levels have to compute hydbid or not (true/false)", metavar="VAR")
+parser.add_option("-z", "--z_levels", type='choice', dest="znivs",
+  choices=['param', 'tropo', 'strato', 'read'], 
+  help="which kind of vertical levels have to be computed ('param', 'tropo', 'strato', 'read')", 
+  metavar="VAR")
 
 (opts, args) = parser.parse_args()
@@ -317 +2108 @@
 dimy = int(opts.dims.split(',')[1])
 dimz = int(opts.dims.split(',')[2])
-
-if opts.hybrid == 'true':
-    hybrid_calc = True
-else:
-    hybrid_calc = False
 
 ofile = 'iniaqua.nc'
@@ -373 +2159 @@
 gam_pv=4.
 
-presnivs, pseudoalt = presnivs_calc(hybrid_calc, dimz)
-lon, lat = global_lonlat(dx,dy)
-lonu, latu = global_lonlat(dx+1,dy)
-lonv, latv = global_lonlat(dx,dy+1)
+# For extra-terrestrial planets
+#presnivs, pseudoalt = presnivs_calc_noterre(opts.znivs, dimz)
+presnivs = presnivs_calc(opts.znivs, dimz)
+lon, lat = global_lonlat(dimy,dimx)
+lonu, latu = global_lonlat(dimy,dimx+1)
+lonv, latv = global_lonlat(dimy+1,dimx)
 
 # 2. Initialize fields towards which to relax
@@ -383 +2171 @@
 knewt_t = np.zeros((dimz), dtype=np.float)
 kfrict = np.zeros((dimz), dtype=np.float)
-clat4 =  np.zeros((dimy+1, dimx), dtype=np.float)
+clat4 =  np.zeros((dimy+1, dimx+1), dtype=np.float)
 
 # Friction
@@ -393 +2181 @@
 
     for j in 1,range(dimy+1):
-        clat4[j,:]=np.cos(rlatu[j])**4
+        clat4[j,:]=np.cos(latu[j,0])**4
 
 # Vertical profile
-tetajl =  np.zeros((dimy+1, dimx), dtype=np.float)
+tetajl =  np.zeros((dimz, dimy+1, dimx), dtype=np.float)
 teta = np.zeros((dimy+1, dimx+1), dtype=np.float)
-tetarappel = np.zeros((dimy+1, dimx+1), dtype=np.float)
-
-for l in range (dz):
+tetarappel = np.zeros((dimz, dimy+1, dimx+1), dtype=np.float)
+
+for l in range (dimz):
     zsig = presnivs[l]/preff
     tetastrat = ttp*zsig**(-kappa)
@@ -411 +2199 @@
     if planet_type == 'giant':
         tetajl[l,:,:] = teta0+(delt_y*((np.sin(latu*3.14159*eps+0.0001))**2) /       \
-          ((rlatu*3.14159*eps+0.0001)**2))-delt_z*np.log(zsig)
+          ((latu*3.14159*eps+0.0001)**2))-delt_z*np.log(zsig)
 
 # Profile stratospheric isotherm (+vortex)
-    w_pv=(1.-np.tanh((rlatu-phi_pv)/dphi_pv))/2.
-    tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
-    for iy in range(dy):
-        for ix in range(dx):
-            tetajl[l,iy,ix]=np.max([tetajl(l,ix,iy),tetastrat])
-
-for iz in range(dimz):
-    for iy in range(dimy+1):
-        tetarappel[iz,iy,0:dimx+1] = tetajl[iz,iy,:]
-
-    tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]
-
+    for iy in range(dimy):
+        for ix in range(dimx):
+            w_pv=(1.-np.tanh((latu[iy,ix]-phi_pv)/dphi_pv))/2.
+            tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+            tetajl[l,iy,ix]=np.max([tetajl[l,iy,ix],tetastrat])
+
+#for iz in range(dimz):
+#    for iy in range(dimy+1):
+#        tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,:]
+#
+#    tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]
+tetarappel = tetajl
+
+# 3. Initialize fields (if necessary)
+# surface pressure
+ps =  np.zeros((dimy, dimx), dtype=np.float)
+
+if iflag_phys > 2:
+# specific value for CMIP5 aqua/terra planets
+# "Specify the initial dry mass to be equivalent to
+#  a global mean surface pressure (101325 minus 245) Pa."
+    ps = 101080.  
+else:
+# use reference surface pressure
+    ps = preff
+        
+# ground geopotential
+phiss =  np.zeros((dimy, dimx), dtype=np.float)
+
+p = pression(dimx, dimy, dimz, ap, bp, ps)
+
+if opts.pexner == 'hybdrid':
+    pks, pk, pkf = exner_hyb(ip1jmp1, ps, p, alpha, beta)
+else:
+    pks, pk, pkf = exner_milieu(ip1jmp1, ps, p, beta)
+
+masse = massdair(p)
+
+# bulk initialization of temperature
+teta = tetarappel.copy()
+
+# geopotential
+phis =  np.zeros((dimy, dimx), dtype=np.float)
+phi =  np.zeros((dimz, dimy, dimx), dtype=np.float)
+
+phis, phi = geopot(ip1jmp1,teta,pk,pks)
+
+# winds
+ucov =  np.zeros((dimz, dimy, dimx), dtype=np.float)
+vcov =  np.zeros((dimz, dimy, dimx), dtype=np.float)
+
+
+if ok_geost:
+    ucov = ugeostr(phi)
+
+# bulk initialization of tracers
+q = np.zeros((dimz, dimy, dimx, nqtot), dtype=np.float)
+
+if planet_type == 'earth':
+# Earth: first two tracers will be water
+    for i in range(nqtot):
+        if i == 1: q[:,:,i] = 1.e-10
+        if i == 2: q[:,:,i] = 1.e-15
+        if i > 2: q[:,:,i] = 0.
+
+# add random perturbation to temperature
+idum  = -1
+zz = RAN1(idum)
+idum  = 0
+for l in range(llm):
+    for ij in range(iip2,ip1jm):
+        teta[ij,l] = teta[ij,l]*(1.+0.005*RAN1(idum))
+
+# maintain periodicity in longitude
+for l in range(llm):
+    for ij in range(0,ip1jmp1,iip1):
+        teta[ij+iim,l]=teta[ij,l]
 
 ncf = NetCDFFile(ofile, 'w')

trunk/tools/lmdz_const.py

@@ -29 +29 @@
 # ihf: (W/m2) Intrinsic heat flux (for giant planets)
 
+planet_type = 'terre'
+iflag_phys = 2
 daysec = 86400.
 preff = 1013250.
@@ -144 +146 @@
 r5les = r3les*(tt-r4les)
 r5ies = r3ies*(tt-r4ies)
+
+# For filtreg
+#
+type_sddu=1
+type_sddv=2
+type_unsddu=3
+type_unsddv=4