Index: LMDZ5/trunk/libf/dyn3d/comvert.h
===================================================================
--- LMDZ5/trunk/libf/dyn3d/comvert.h	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/comvert.h	(revision 1520)
@@ -5,8 +5,28 @@
 !   INCLUDE 'comvert.h'
 
-      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
-     &               pa,preff,nivsigs(llm),nivsig(llm+1)
+      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
+     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
+     &               aps(llm),bps(llm),scaleheight
 
-      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
+      common/comverti/disvert_type
+
+      real ap     ! hybrid pressure contribution at interlayers
+      real bp     ! hybrid sigma contribution at interlayer
+      real presnivs ! (reference) pressure at mid-layers
+      real dpres
+      real pa     ! reference pressure (Pa) at which hybrid coordinates
+                  ! become purely pressure
+      real preff  ! reference surface pressure (Pa)
+      real nivsigs
+      real nivsig
+      real aps    ! hybrid pressure contribution at mid-layers
+      real bps    ! hybrid sigma contribution at mid-layers
+      real scaleheight ! atmospheric (reference) scale height (km)
+
+      integer disvert_type ! type of vertical discretization:
+                           ! 1: Earth (default for planet_type==earth),
+                           !     automatic generation
+                           ! 2: Planets (default for planet_type!=earth),
+                           !     using 'z2sig.def' (or 'esasig.def) file
 
  !-----------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dyn3d/disvert.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/disvert.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/disvert.F90	(revision 1520)
@@ -1,5 +1,5 @@
 ! $Id$
 
-SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
@@ -17,7 +17,9 @@
   ! ds(l) : distance entre les couches l et l-1 en coord.s
 
-  REAL pa, preff
-  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
-  REAL presnivs(llm)
+  real,intent(in) :: pa, preff
+  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
+  real,intent(out) :: presnivs(llm)
+  real,intent(out) :: scaleheight
 
   REAL sig(llm+1), dsig(llm)
@@ -26,11 +28,13 @@
   INTEGER l
   REAL dsigmin
-  REAL alpha, beta, deltaz, h
+  REAL alpha, beta, deltaz,scaleheight
   INTEGER iostat
-  REAL pi, x
+  REAL x
+  character(len=*),parameter :: modname="disvert"
 
   !-----------------------------------------------------------------------
 
-  pi = 2 * ASIN(1.)
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
 
   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
@@ -38,5 +42,5 @@
   IF (iostat == 0) THEN
      ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
      READ(99, *) beta ! facteur d'acroissement en haut 1.3
@@ -44,6 +48,7 @@
      READ(99, *) k1 ! nombre de couches dans la transition haute
      CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'h, alpha, k0, k1, beta'
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
 
      alpha=deltaz/tanh(1./k0)*2.
@@ -51,7 +56,8 @@
      sig(1)=1.
      do l=1, llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
-             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
+             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
+                  *beta**(l-(llm-k1))/log(beta))
+        zk=-scaleheight*log(sig(l+1))
 
         dzk1=alpha*tanh(l/k0)
@@ -73,8 +79,10 @@
            dsigmin=1.
         else
-           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119,23 +127,23 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)' BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)' AP '
+  write(lunout, *)  trim(modname),': AP '
   write(lunout, *) ap
 
   write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
   write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
-  write(lunout, *)'8km'
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  write(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*8. &
-          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
           max(ap(l+1)+bp(l+1)*preff, 1.e-10))
   ENDDO
 
-  write(lunout, *) 'PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs
 
Index: LMDZ5/trunk/libf/dyn3d/disvert_noterre.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/disvert_noterre.F	(revision 1520)
+++ LMDZ5/trunk/libf/dyn3d/disvert_noterre.F	(revision 1520)
@@ -0,0 +1,331 @@
+! $Id: $
+      SUBROUTINE disvert_noterre
+
+c    Auteur :  F. Forget Y. Wanherdrick, P. Levan
+c    Nouvelle version 100% Mars !!
+c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
+
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
+
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "logic.h"
+#include "iniprint.h"
+c
+c=======================================================================
+c    Discretisation verticale en coordonnée hybride (ou sigma)
+c
+c=======================================================================
+c
+c   declarations:
+c   -------------
+c
+c
+      INTEGER l,ll
+      REAL snorm
+      REAL alpha,beta,gama,delta,deltaz
+      real quoi,quand
+      REAL zsig(llm),sig(llm+1)
+      INTEGER np,ierr
+      integer :: ierr1,ierr2,ierr3,ierr4
+      REAL x
+
+      REAL SSUM
+      EXTERNAL SSUM
+      real newsig 
+      REAL dz0,dz1,nhaut,sig1,esig,csig,zz
+      real tt,rr,gg, prevz
+      real s(llm),dsig(llm) 
+      real pseudoalt(llm)
+
+      integer iz 
+      real z, ps,p
+      character(len=*),parameter :: modname="disvert_noterre"
+
+c
+c-----------------------------------------------------------------------
+c
+! Initializations:
+!      pi=2.*ASIN(1.) ! already done in iniconst
+      
+      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
+      CALL getin('hybrid',hybrid)
+      write(lunout,*) trim(modname),': hybrid=',hybrid
+
+! Ouverture possible de fichiers typiquement E.T.
+
+         open(99,file="esasig.def",status='old',form='formatted',
+     s   iostat=ierr2)
+         if(ierr2.ne.0) then
+              close(99)
+              open(99,file="z2sig.def",status='old',form='formatted',
+     s        iostat=ierr4)
+         endif
+
+c-----------------------------------------------------------------------
+c   cas 1 on lit les options dans esasig.def:
+c   ----------------------------------------
+
+      IF(ierr2.eq.0) then
+
+c        Lecture de esasig.def :
+c        Systeme peu souple, mais qui respecte en theorie
+c        La conservation de l'energie (conversion Energie potentielle
+c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
+
+         write(lunout,*)'*****************************'
+         write(lunout,*)'WARNING reading esasig.def'
+         write(lunout,*)'*****************************'
+         READ(99,*) scaleheight
+         READ(99,*) dz0
+         READ(99,*) dz1
+         READ(99,*) nhaut
+         CLOSE(99)
+
+         dz0=dz0/scaleheight
+         dz1=dz1/scaleheight
+
+         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
+
+         esig=1.
+
+         do l=1,20
+            esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
+         enddo
+         csig=(1./sig1-1.)/(exp(esig)-1.)
+
+         DO L = 2, llm
+            zz=csig*(exp(esig*(l-1.))-1.)
+            sig(l) =1./(1.+zz)
+     &      * tanh(.5*(llm+1-l)/nhaut)
+         ENDDO
+         sig(1)=1.
+         sig(llm+1)=0.
+         quoi      = 1. + 2.* kappa
+         s( llm )  = 1.
+         s(llm-1) = quoi
+         IF( llm.gt.2 )  THEN
+            DO  ll = 2, llm-1
+               l         = llm+1 - ll
+               quand     = sig(l+1)/ sig(l)
+               s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
+            ENDDO
+         END IF
+c
+         snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
+         DO l = 1, llm
+            s(l)    = s(l)/ snorm
+         ENDDO
+
+c-----------------------------------------------------------------------
+c   cas 2 on lit les options dans z2sig.def:
+c   ----------------------------------------
+
+      ELSE IF(ierr4.eq.0) then
+         write(lunout,*)'****************************'
+         write(lunout,*)'Reading z2sig.def'
+         write(lunout,*)'****************************'
+
+         READ(99,*) scaleheight
+         do l=1,llm
+            read(99,*) zsig(l)
+         end do
+         CLOSE(99)
+
+         sig(1) =1
+         do l=2,llm
+           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + 
+     &                      exp(-zsig(l-1)/scaleheight) )
+         end do
+         sig(llm+1) =0
+
+c-----------------------------------------------------------------------
+      ELSE
+         write(lunout,*) 'didn t you forget something ??? '
+         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
+         stop
+      ENDIF
+c-----------------------------------------------------------------------
+
+      DO l=1,llm
+        nivsigs(l) = REAL(l)
+      ENDDO
+
+      DO l=1,llmp1
+        nivsig(l)= REAL(l)
+      ENDDO
+
+ 
+c-----------------------------------------------------------------------
+c    ....  Calculs  de ap(l) et de bp(l)  ....
+c    .........................................
+c
+c   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
+c-----------------------------------------------------------------------
+c
+
+      if (hybrid) then  ! use hybrid coordinates
+         write(lunout,*) "*********************************"
+         write(lunout,*) "Using hybrid vertical coordinates"
+         write(lunout,*) 
+c        Coordonnees hybrides avec mod
+         DO l = 1, llm
+
+         call sig_hybrid(sig(l),pa,preff,newsig)
+            bp(l) = EXP( 1. - 1./(newsig**2)  )
+            ap(l) = pa * (newsig - bp(l) )
+         enddo
+         bp(llmp1) = 0.
+         ap(llmp1) = 0.
+      else ! use sigma coordinates
+         write(lunout,*) "********************************"
+         write(lunout,*) "Using sigma vertical coordinates"
+         write(lunout,*) 
+c        Pour ne pas passer en coordonnees hybrides
+         DO l = 1, llm
+            ap(l) = 0.
+            bp(l) = sig(l)
+         ENDDO
+         ap(llmp1) = 0.
+      endif
+
+      bp(llmp1) =   0.
+
+      write(lunout,*) trim(modname),': BP '
+      write(lunout,*)  bp
+      write(lunout,*) trim(modname),': AP '
+      write(lunout,*)  ap
+
+c     Calcul au milieu des couches :
+c     WARNING : le choix de placer le milieu des couches au niveau de
+c     pression intermédiaire est arbitraire et pourrait etre modifié.
+c     Le calcul du niveau pour la derniere couche 
+c     (on met la meme distance (en log pression)  entre P(llm)
+c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
+c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
+
+      DO l = 1, llm-1
+       aps(l) =  0.5 *( ap(l) +ap(l+1)) 
+       bps(l) =  0.5 *( bp(l) +bp(l+1)) 
+      ENDDO
+     
+      if (hybrid) then
+         aps(llm) = aps(llm-1)**2 / aps(llm-2) 
+         bps(llm) = 0.5*(bp(llm) + bp(llm+1))
+      else
+         bps(llm) = bps(llm-1)**2 / bps(llm-2) 
+         aps(llm) = 0.
+      end if
+
+      write(lunout,*) trim(modname),': BPs '
+      write(lunout,*)  bps
+      write(lunout,*) trim(modname),': APs'
+      write(lunout,*)  aps
+
+      DO l = 1, llm
+       presnivs(l) = aps(l)+bps(l)*preff
+       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
+      ENDDO
+
+      write(lunout,*)trim(modname),' : PRESNIVS' 
+      write(lunout,*)presnivs 
+      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
+     &                'height of ',scaleheight,' km)' 
+      write(lunout,*)pseudoalt
+
+c     --------------------------------------------------
+c     This can be used to plot the vertical discretization
+c     (> xmgrace -nxy testhybrid.tab )
+c     --------------------------------------------------
+c     open (53,file='testhybrid.tab')
+c     scaleheight=15.5
+c     do iz=0,34
+c       z = -5 + min(iz,34-iz)
+c     approximation of scale height for Venus
+c       scaleheight = 15.5 - z/55.*10.
+c       ps = preff*exp(-z/scaleheight)
+c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
+c       do l=2,llm
+c     approximation of scale height for Venus
+c          if (zsig(l-1).le.55.) then
+c             scaleheight = 15.5 - zsig(l-1)/55.*10.
+c          else
+c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
+c          endif
+c          zsig(l)= zsig(l-1)-scaleheight*
+c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
+c       end do
+c       write(53,'(I3,50F10.5)') iz, zsig
+c      end do
+c      close(53)
+c     --------------------------------------------------
+
+
+      RETURN
+      END
+
+c ************************************************************
+      subroutine sig_hybrid(sig,pa,preff,newsig)
+c     ----------------------------------------------
+c     Subroutine utilisee pour calculer des valeurs de sigma modifie
+c     pour conserver les coordonnees verticales decrites dans
+c     esasig.def/z2sig.def lors du passage en coordonnees hybrides
+c     F. Forget 2002
+c     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
+c     L'objectif est de calculer newsig telle que
+c       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
+c     Cela ne se résoud pas analytiquement: 
+c     => on résoud par iterration bourrine 
+c     ----------------------------------------------
+c     Information  : where exp(1-1./x**2) become << x
+c           x      exp(1-1./x**2) /x
+c           1           1
+c           0.68       0.5
+c           0.5        1.E-1
+c           0.391      1.E-2
+c           0.333      1.E-3
+c           0.295      1.E-4
+c           0.269      1.E-5
+c           0.248      1.E-6
+c        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
+
+
+      implicit none
+      real x1, x2, sig,pa,preff, newsig, F
+      integer j
+
+      newsig = sig
+      x1=0
+      x2=1
+      if (sig.ge.1) then
+            newsig= sig
+      else if (sig*preff/pa.ge.0.25) then
+        DO J=1,9999  ! nombre d''iteration max
+          F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
+c         write(0,*) J, ' newsig =', newsig, ' F= ', F
+          if (F.gt.1) then
+              X2 = newsig
+              newsig=(X1+newsig)*0.5
+          else
+              X1 = newsig
+              newsig=(X2+newsig)*0.5
+          end if
+c         Test : on arete lorsque on approxime sig à moins de 0.01 m près 
+c         (en pseudo altitude) :
+          IF(abs(10.*log(F)).LT.1.E-5) goto 999
+        END DO
+       else   !    if (sig*preff/pa.le.0.25) then
+             newsig= sig*preff/pa
+       end if
+ 999   continue
+       Return
+      END
Index: LMDZ5/trunk/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/etat0_netcdf.F90	(revision 1520)
@@ -251,5 +251,9 @@
 !*******************************************************************************
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  if (disvert_type.eq.1) then
+    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  else ! we assume that we are in the disvert_type==2 case
+    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
+  endif
   pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
 !  WRITE(lunout,*) 'P3D :', p3d(10,20,:)
Index: LMDZ5/trunk/libf/dyn3d/exner_hyb.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/exner_hyb.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/exner_hyb.F	(revision 1520)
@@ -51,17 +51,32 @@
       REAL SSUM
 c
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_hyb"
+      
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
 
       if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
-        endif
         
         ! Compute pks(:),pk(:),pkf(:)
@@ -77,6 +92,8 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
+!!!! General case:
      
       unpl2k    = 1.+ 2.* kappa
Index: LMDZ5/trunk/libf/dyn3d/exner_milieu.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/exner_milieu.F	(revision 1520)
+++ LMDZ5/trunk/libf/dyn3d/exner_milieu.F	(revision 1520)
@@ -0,0 +1,151 @@
+!
+! $Id $
+!
+      SUBROUTINE  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
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "comvert.h"
+#include "serre.h"
+
+      INTEGER  ngrid
+      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
+      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
+
+c    .... variables locales   ...
+
+      INTEGER l, ij
+      REAL dum1
+
+      REAL ppn(iim),pps(iim)
+      REAL xpn, xps
+      REAL SSUM
+      EXTERNAL SSUM
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_milieu"
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.2) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==2",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
+      if (llm.eq.1) then
+      
+        ! Compute pks(:),pk(:),pkf(:)
+        
+        DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij) 
+          pk(ij,1) = .5*pks(ij)
+        ENDDO
+        
+        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+        ! our work is done, exit routine
+        return
+
+      endif ! of if (llm.eq.1)
+
+!!!! General case:
+
+c     -------------
+c     Calcul de pks
+c     -------------
+   
+      DO   ij  = 1, ngrid
+        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+      ENDDO
+
+      DO  ij   = 1, iim
+        ppn(ij) = aire(   ij   ) * pks(  ij     )
+        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
+      ENDDO
+      xpn      = SSUM(iim,ppn,1) /apoln
+      xps      = SSUM(iim,pps,1) /apols
+
+      DO ij   = 1, iip1
+        pks(   ij     )  =  xpn
+        pks( ij+ip1jm )  =  xps
+      ENDDO
+c
+c
+c    .... Calcul de pk  pour la couche l 
+c    --------------------------------------------
+c
+      dum1 = cpp * (2*preff)**(-kappa) 
+      DO l = 1, llm-1
+        DO   ij   = 1, ngrid
+         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+        ENDDO
+      ENDDO
+
+c    .... Calcul de pk  pour la couche l = llm ..
+c    (on met la meme distance (en log pression)  entre Pk(llm)
+c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+      DO   ij   = 1, ngrid
+         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+      ENDDO
+
+
+c    calcul de pkf
+c    -------------
+      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+      
+c    EST-CE UTILE ?? : calcul de beta
+c    --------------------------------
+      DO l = 2, llm
+        DO   ij   = 1, ngrid
+          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
+        ENDDO
+      ENDDO
+
+      RETURN
+      END
Index: LMDZ5/trunk/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/guide_mod.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/guide_mod.F90	(revision 1520)
@@ -644,6 +644,9 @@
 ! -----------------------------------------------------------------
     CALL pression( ip1jmp1, ap, bp, psi, p )
-    CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
-
+    if (disvert_type==1) then
+      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+    else ! we assume that we are in the disvert_type==2 case
+      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
+    endif
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
     unskap=1./kappa
Index: LMDZ5/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 1520)
@@ -71,4 +71,7 @@
 
   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  
+  character(len=*),parameter :: modname="iniacademic"
+  character(len=80) :: abort_message
 
   !-----------------------------------------------------------------------
@@ -210,5 +213,13 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        if (disvert_type.eq.1) then
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        elseif (disvert_type.eq.2) then
+          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        else
+          write(abort_message,*) "Wrong value for disvert_type: ", &
+                              disvert_type
+          call abort_gcm(modname,abort_message,0)
+        endif
         CALL massdair(p,masse)
 
Index: LMDZ5/trunk/libf/dyn3d/iniconst.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/iniconst.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/iniconst.F	(revision 1520)
@@ -5,4 +5,10 @@
 
       USE control_mod
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
 
       IMPLICIT NONE
@@ -21,5 +27,6 @@
 #include "iniprint.h"
 
-
+      character(len=*),parameter :: modname="iniconst"
+      character(len=80) :: abort_message
 c
 c
@@ -48,11 +55,30 @@
       r       = cpp * kappa
 
-      write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
+      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
 c
 c-----------------------------------------------------------------------
 
-       CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-c
-c
-       RETURN
-       END
+! vertical discretization: default behavior depends on planet_type flag
+      if (planet_type=="earth") then
+        disvert_type=1
+      else
+        disvert_type=2
+      endif
+      ! but user can also specify using one or the other in run.def:
+      call getin('disvert_type',disvert_type)
+      write(lunout,*) trim(modname),': disvert_type=',disvert_type
+      
+      if (disvert_type==1) then
+       ! standard case for Earth (automatic generation of levels)
+       call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     &              scaleheight)
+      else if (disvert_type==2) then
+        ! standard case for planets (levels generated using z2sig.def file)
+        call disvert_noterre
+      else
+        write(abort_message,*) "Wrong value for disvert_type: ",
+     &                        disvert_type
+        call abort_gcm(modname,abort_message,0)
+      endif
+
+      END
Index: LMDZ5/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 1520)
@@ -172,5 +172,5 @@
 
       character*80 dynhist_file, dynhistave_file
-      character(len=20) :: modname
+      character(len=*),parameter :: modname="leapfrog"
       character*80 abort_message
 
@@ -193,5 +193,4 @@
       itaufin   = nday*day_step
       itaufinp1 = itaufin +1
-      modname="leapfrog"
       
 
@@ -211,5 +210,9 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+      if (disvert_type==1) then
+        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+      else ! we assume that we are in the disvert_type==2 case
+        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+      endif
 
 c-----------------------------------------------------------------------
@@ -359,5 +362,9 @@
 
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         if (disvert_type==1) then
+           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         else ! we assume that we are in the disvert_type==2 case
+           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+         endif
 
 !           rdaym_ini  = itau * dtvr / daysec
@@ -463,5 +470,9 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        if (disvert_type==1) then
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        else ! we assume that we are in the disvert_type==2 case
+          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        endif
 
 
Index: LMDZ5/trunk/libf/dyn3d/limy.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/limy.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/limy.F	(revision 1520)
@@ -1,5 +1,5 @@
-!
-! $Header$
-!
+c
+c $Id$
+c
       SUBROUTINE limy(s0,sy,sm,pente_max)
 c
@@ -40,5 +40,5 @@
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
+      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
       Logical extremum,first
       save first
@@ -117,14 +117,14 @@
 
 c     print*,dyqv(iip1+1)
-c     apn=abs(dyq(1)/dyqv(iip1+1))
+c     appn=abs(dyq(1)/dyqv(iip1+1))
 c     print*,dyq(ip1jm+1)
 c     print*,dyqv(ip1jm-iip1+1)
-c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 c     do ij=2,iim
-c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 c     enddo
-c     apn=min(pente_max/apn,1.)
-c     aps=min(pente_max/aps,1.)
+c     appn=min(pente_max/appn,1.)
+c     apps=min(pente_max/apps,1.)
 
 
@@ -132,13 +132,13 @@
 
 c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-c    &   apn=0.
+c    &   appn=0.
 c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-c    &   aps=0.
+c    &   apps=0.
 
 c   limitation des pentes aux poles
 c     do ij=1,iip1
-c        dyq(ij)=apn*dyq(ij)
-c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+c        dyq(ij)=appn*dyq(ij)
+c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 c     enddo
 
Index: LMDZ5/trunk/libf/dyn3d/logic.h
===================================================================
--- LMDZ5/trunk/libf/dyn3d/logic.h	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/logic.h	(revision 1520)
@@ -3,18 +3,23 @@
 !
 !
-!
+! NB: keep items of different kinds in seperate common blocs to avoid
+!     "misaligned commons" issues
 !-----------------------------------------------------------------------
 ! INCLUDE 'logic.h'
 
-      COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
+      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
 
+      COMMON/logici/ iflag_phys,iflag_trac
+      
       LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
      &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
+                     ! (only used if disvert_type==2)
 
-      INTEGER iflag_phys
+      integer iflag_phys,iflag_trac
 !-----------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dyn3d/vlsplt.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/vlsplt.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/vlsplt.F	(revision 1520)
@@ -1,6 +1,4 @@
-!
-! $Header$
-!
-c
+c
+c $Id$
 c
 
@@ -478,5 +476,5 @@
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
 c     REAL newq,oldmasse
       Logical extremum,first,testcpu
@@ -602,14 +600,14 @@
 C     PRINT*,dyq(1)
 C     PRINT*,dyqv(iip1+1)
-C     apn=abs(dyq(1)/dyqv(iip1+1))
+C     appn=abs(dyq(1)/dyqv(iip1+1))
 C     PRINT*,dyq(ip1jm+1)
 C     PRINT*,dyqv(ip1jm-iip1+1)
-C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 C     DO ij=2,iim
-C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 C     ENDDO
-C     apn=min(pente_max/apn,1.)
-C     aps=min(pente_max/aps,1.)
+C     appn=min(pente_max/appn,1.)
+C     apps=min(pente_max/apps,1.)
 C
 C
@@ -617,13 +615,13 @@
 C
 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-C    &   apn=0.
+C    &   appn=0.
 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-C    &   aps=0.
+C    &   apps=0.
 C
 C   limitation des pentes aux poles
 C     DO ij=1,iip1
-C        dyq(ij)=apn*dyq(ij)
-C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+C        dyq(ij)=appn*dyq(ij)
+C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 C     ENDDO
 C
Index: LMDZ5/trunk/libf/dyn3d/vlspltqs.F
===================================================================
--- LMDZ5/trunk/libf/dyn3d/vlspltqs.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3d/vlspltqs.F	(revision 1520)
@@ -1,5 +1,5 @@
-!
-! $Header$
-!
+c
+c $Id$
+c
        SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
      ,                                  p,pk,teta                 )
@@ -634,14 +634,14 @@
 C     PRINT*,dyq(1)
 C     PRINT*,dyqv(iip1+1)
-C     apn=abs(dyq(1)/dyqv(iip1+1))
+C     appn=abs(dyq(1)/dyqv(iip1+1))
 C     PRINT*,dyq(ip1jm+1)
 C     PRINT*,dyqv(ip1jm-iip1+1)
-C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 C     DO ij=2,iim
-C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 C     ENDDO
-C     apn=min(pente_max/apn,1.)
-C     aps=min(pente_max/aps,1.)
+C     appn=min(pente_max/appn,1.)
+C     apps=min(pente_max/apps,1.)
 C
 C
@@ -649,13 +649,13 @@
 C
 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-C    &   apn=0.
+C    &   appn=0.
 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-C    &   aps=0.
+C    &   apps=0.
 C
 C   limitation des pentes aux poles
 C     DO ij=1,iip1
-C        dyq(ij)=apn*dyq(ij)
-C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+C        dyq(ij)=appn*dyq(ij)
+C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 C     ENDDO
 C
Index: LMDZ5/trunk/libf/dyn3dpar/comvert.h
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/comvert.h	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/comvert.h	(revision 1520)
@@ -1,12 +1,32 @@
 !
-! $Header$
+! $Id$
 !
 !-----------------------------------------------------------------------
 !   INCLUDE 'comvert.h'
 
-      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),       &
-     &               pa,preff,nivsigs(llm),nivsig(llm+1)
+      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
+     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
+     &               aps(llm),bps(llm),scaleheight
 
-      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
+      common/comverti/disvert_type
 
-!-----------------------------------------------------------------------
+      real ap     ! hybrid pressure contribution at interlayers
+      real bp     ! hybrid sigma contribution at interlayer
+      real presnivs ! (reference) pressure at mid-layers
+      real dpres
+      real pa     ! reference pressure (Pa) at which hybrid coordinates
+                  ! become purely pressure
+      real preff  ! reference surface pressure (Pa)
+      real nivsigs
+      real nivsig
+      real aps    ! hybrid pressure contribution at mid-layers
+      real bps    ! hybrid sigma contribution at mid-layers
+      real scaleheight ! atmospheric (reference) scale height (km)
+
+      integer disvert_type ! type of vertical discretization:
+                           ! 1: Earth (default for planet_type==earth),
+                           !     automatic generation
+                           ! 2: Planets (default for planet_type!=earth),
+                           !     using 'z2sig.def' (or 'esasig.def) file
+
+ !-----------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dyn3dpar/disvert.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/disvert.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/disvert.F90	(revision 1520)
@@ -1,5 +1,5 @@
 ! $Id$
 
-SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
@@ -17,7 +17,9 @@
   ! ds(l) : distance entre les couches l et l-1 en coord.s
 
-  REAL pa, preff
-  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
-  REAL presnivs(llm)
+  real,intent(in) :: pa, preff
+  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
+  real,intent(out) :: presnivs(llm)
+  real,intent(out) :: scaleheight
 
   REAL sig(llm+1), dsig(llm)
@@ -26,11 +28,13 @@
   INTEGER l
   REAL dsigmin
-  REAL alpha, beta, deltaz, h
+  REAL alpha, beta, deltaz,scaleheight
   INTEGER iostat
-  REAL pi, x
+  REAL x
+  character(len=*),parameter :: modname="disvert"
 
   !-----------------------------------------------------------------------
 
-  pi = 2 * ASIN(1.)
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
 
   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
@@ -38,5 +42,5 @@
   IF (iostat == 0) THEN
      ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
      READ(99, *) beta ! facteur d'acroissement en haut 1.3
@@ -44,6 +48,7 @@
      READ(99, *) k1 ! nombre de couches dans la transition haute
      CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'h, alpha, k0, k1, beta'
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
 
      alpha=deltaz/tanh(1./k0)*2.
@@ -51,7 +56,8 @@
      sig(1)=1.
      do l=1, llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
-             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
+             *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
+                  *beta**(l-(llm-k1))/log(beta))
+        zk=-scaleheight*log(sig(l+1))
 
         dzk1=alpha*tanh(l/k0)
@@ -73,8 +79,10 @@
            dsigmin=1.
         else
-           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
            dsigmin=1.
         endif
-        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
      endif
 
@@ -119,23 +127,23 @@
   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
-  write(lunout, *)' BP '
+  write(lunout, *)  trim(modname),': BP '
   write(lunout, *) bp
-  write(lunout, *)' AP '
+  write(lunout, *)  trim(modname),': AP '
   write(lunout, *) ap
 
   write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
   write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
-  write(lunout, *)'8km'
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  write(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*8. &
-          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
           max(ap(l+1)+bp(l+1)*preff, 1.e-10))
   ENDDO
 
-  write(lunout, *) 'PRESNIVS '
+  write(lunout, *) trim(modname),': PRESNIVS '
   write(lunout, *) presnivs
 
Index: LMDZ5/trunk/libf/dyn3dpar/disvert_noterre.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/disvert_noterre.F	(revision 1520)
+++ LMDZ5/trunk/libf/dyn3dpar/disvert_noterre.F	(revision 1520)
@@ -0,0 +1,331 @@
+! $Id: $
+      SUBROUTINE disvert_noterre
+
+c    Auteur :  F. Forget Y. Wanherdrick, P. Levan
+c    Nouvelle version 100% Mars !!
+c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
+
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
+
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "logic.h"
+#include "iniprint.h"
+c
+c=======================================================================
+c    Discretisation verticale en coordonnée hybride (ou sigma)
+c
+c=======================================================================
+c
+c   declarations:
+c   -------------
+c
+c
+      INTEGER l,ll
+      REAL snorm
+      REAL alpha,beta,gama,delta,deltaz
+      real quoi,quand
+      REAL zsig(llm),sig(llm+1)
+      INTEGER np,ierr
+      integer :: ierr1,ierr2,ierr3,ierr4
+      REAL x
+
+      REAL SSUM
+      EXTERNAL SSUM
+      real newsig 
+      REAL dz0,dz1,nhaut,sig1,esig,csig,zz
+      real tt,rr,gg, prevz
+      real s(llm),dsig(llm) 
+      real pseudoalt(llm)
+
+      integer iz 
+      real z, ps,p
+      character(len=*),parameter :: modname="disvert_noterre"
+
+c
+c-----------------------------------------------------------------------
+c
+! Initializations:
+!      pi=2.*ASIN(1.) ! already done in iniconst
+      
+      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
+      CALL getin('hybrid',hybrid)
+      write(lunout,*) trim(modname),': hybrid=',hybrid
+
+! Ouverture possible de fichiers typiquement E.T.
+
+         open(99,file="esasig.def",status='old',form='formatted',
+     s   iostat=ierr2)
+         if(ierr2.ne.0) then
+              close(99)
+              open(99,file="z2sig.def",status='old',form='formatted',
+     s        iostat=ierr4)
+         endif
+
+c-----------------------------------------------------------------------
+c   cas 1 on lit les options dans esasig.def:
+c   ----------------------------------------
+
+      IF(ierr2.eq.0) then
+
+c        Lecture de esasig.def :
+c        Systeme peu souple, mais qui respecte en theorie
+c        La conservation de l'energie (conversion Energie potentielle
+c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
+
+         write(lunout,*)'*****************************'
+         write(lunout,*)'WARNING reading esasig.def'
+         write(lunout,*)'*****************************'
+         READ(99,*) scaleheight
+         READ(99,*) dz0
+         READ(99,*) dz1
+         READ(99,*) nhaut
+         CLOSE(99)
+
+         dz0=dz0/scaleheight
+         dz1=dz1/scaleheight
+
+         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
+
+         esig=1.
+
+         do l=1,20
+            esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.)
+         enddo
+         csig=(1./sig1-1.)/(exp(esig)-1.)
+
+         DO L = 2, llm
+            zz=csig*(exp(esig*(l-1.))-1.)
+            sig(l) =1./(1.+zz)
+     &      * tanh(.5*(llm+1-l)/nhaut)
+         ENDDO
+         sig(1)=1.
+         sig(llm+1)=0.
+         quoi      = 1. + 2.* kappa
+         s( llm )  = 1.
+         s(llm-1) = quoi
+         IF( llm.gt.2 )  THEN
+            DO  ll = 2, llm-1
+               l         = llm+1 - ll
+               quand     = sig(l+1)/ sig(l)
+               s(l-1)    = quoi * (1.-quand) * s(l)  + quand * s(l+1)
+            ENDDO
+         END IF
+c
+         snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2)
+         DO l = 1, llm
+            s(l)    = s(l)/ snorm
+         ENDDO
+
+c-----------------------------------------------------------------------
+c   cas 2 on lit les options dans z2sig.def:
+c   ----------------------------------------
+
+      ELSE IF(ierr4.eq.0) then
+         write(lunout,*)'****************************'
+         write(lunout,*)'Reading z2sig.def'
+         write(lunout,*)'****************************'
+
+         READ(99,*) scaleheight
+         do l=1,llm
+            read(99,*) zsig(l)
+         end do
+         CLOSE(99)
+
+         sig(1) =1
+         do l=2,llm
+           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + 
+     &                      exp(-zsig(l-1)/scaleheight) )
+         end do
+         sig(llm+1) =0
+
+c-----------------------------------------------------------------------
+      ELSE
+         write(lunout,*) 'didn t you forget something ??? '
+         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
+         stop
+      ENDIF
+c-----------------------------------------------------------------------
+
+      DO l=1,llm
+        nivsigs(l) = REAL(l)
+      ENDDO
+
+      DO l=1,llmp1
+        nivsig(l)= REAL(l)
+      ENDDO
+
+ 
+c-----------------------------------------------------------------------
+c    ....  Calculs  de ap(l) et de bp(l)  ....
+c    .........................................
+c
+c   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
+c-----------------------------------------------------------------------
+c
+
+      if (hybrid) then  ! use hybrid coordinates
+         write(lunout,*) "*********************************"
+         write(lunout,*) "Using hybrid vertical coordinates"
+         write(lunout,*) 
+c        Coordonnees hybrides avec mod
+         DO l = 1, llm
+
+         call sig_hybrid(sig(l),pa,preff,newsig)
+            bp(l) = EXP( 1. - 1./(newsig**2)  )
+            ap(l) = pa * (newsig - bp(l) )
+         enddo
+         bp(llmp1) = 0.
+         ap(llmp1) = 0.
+      else ! use sigma coordinates
+         write(lunout,*) "********************************"
+         write(lunout,*) "Using sigma vertical coordinates"
+         write(lunout,*) 
+c        Pour ne pas passer en coordonnees hybrides
+         DO l = 1, llm
+            ap(l) = 0.
+            bp(l) = sig(l)
+         ENDDO
+         ap(llmp1) = 0.
+      endif
+
+      bp(llmp1) =   0.
+
+      write(lunout,*) trim(modname),': BP '
+      write(lunout,*)  bp
+      write(lunout,*) trim(modname),': AP '
+      write(lunout,*)  ap
+
+c     Calcul au milieu des couches :
+c     WARNING : le choix de placer le milieu des couches au niveau de
+c     pression intermédiaire est arbitraire et pourrait etre modifié.
+c     Le calcul du niveau pour la derniere couche 
+c     (on met la meme distance (en log pression)  entre P(llm)
+c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
+c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
+
+      DO l = 1, llm-1
+       aps(l) =  0.5 *( ap(l) +ap(l+1)) 
+       bps(l) =  0.5 *( bp(l) +bp(l+1)) 
+      ENDDO
+     
+      if (hybrid) then
+         aps(llm) = aps(llm-1)**2 / aps(llm-2) 
+         bps(llm) = 0.5*(bp(llm) + bp(llm+1))
+      else
+         bps(llm) = bps(llm-1)**2 / bps(llm-2) 
+         aps(llm) = 0.
+      end if
+
+      write(lunout,*) trim(modname),': BPs '
+      write(lunout,*)  bps
+      write(lunout,*) trim(modname),': APs'
+      write(lunout,*)  aps
+
+      DO l = 1, llm
+       presnivs(l) = aps(l)+bps(l)*preff
+       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
+      ENDDO
+
+      write(lunout,*)trim(modname),' : PRESNIVS' 
+      write(lunout,*)presnivs 
+      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
+     &                'height of ',scaleheight,' km)' 
+      write(lunout,*)pseudoalt
+
+c     --------------------------------------------------
+c     This can be used to plot the vertical discretization
+c     (> xmgrace -nxy testhybrid.tab )
+c     --------------------------------------------------
+c     open (53,file='testhybrid.tab')
+c     scaleheight=15.5
+c     do iz=0,34
+c       z = -5 + min(iz,34-iz)
+c     approximation of scale height for Venus
+c       scaleheight = 15.5 - z/55.*10.
+c       ps = preff*exp(-z/scaleheight)
+c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
+c       do l=2,llm
+c     approximation of scale height for Venus
+c          if (zsig(l-1).le.55.) then
+c             scaleheight = 15.5 - zsig(l-1)/55.*10.
+c          else
+c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
+c          endif
+c          zsig(l)= zsig(l-1)-scaleheight*
+c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
+c       end do
+c       write(53,'(I3,50F10.5)') iz, zsig
+c      end do
+c      close(53)
+c     --------------------------------------------------
+
+
+      RETURN
+      END
+
+c ************************************************************
+      subroutine sig_hybrid(sig,pa,preff,newsig)
+c     ----------------------------------------------
+c     Subroutine utilisee pour calculer des valeurs de sigma modifie
+c     pour conserver les coordonnees verticales decrites dans
+c     esasig.def/z2sig.def lors du passage en coordonnees hybrides
+c     F. Forget 2002
+c     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
+c     L'objectif est de calculer newsig telle que
+c       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
+c     Cela ne se résoud pas analytiquement: 
+c     => on résoud par iterration bourrine 
+c     ----------------------------------------------
+c     Information  : where exp(1-1./x**2) become << x
+c           x      exp(1-1./x**2) /x
+c           1           1
+c           0.68       0.5
+c           0.5        1.E-1
+c           0.391      1.E-2
+c           0.333      1.E-3
+c           0.295      1.E-4
+c           0.269      1.E-5
+c           0.248      1.E-6
+c        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
+
+
+      implicit none
+      real x1, x2, sig,pa,preff, newsig, F
+      integer j
+
+      newsig = sig
+      x1=0
+      x2=1
+      if (sig.ge.1) then
+            newsig= sig
+      else if (sig*preff/pa.ge.0.25) then
+        DO J=1,9999  ! nombre d''iteration max
+          F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
+c         write(0,*) J, ' newsig =', newsig, ' F= ', F
+          if (F.gt.1) then
+              X2 = newsig
+              newsig=(X1+newsig)*0.5
+          else
+              X1 = newsig
+              newsig=(X2+newsig)*0.5
+          end if
+c         Test : on arete lorsque on approxime sig à moins de 0.01 m près 
+c         (en pseudo altitude) :
+          IF(abs(10.*log(F)).LT.1.E-5) goto 999
+        END DO
+       else   !    if (sig*preff/pa.le.0.25) then
+             newsig= sig*preff/pa
+       end if
+ 999   continue
+       Return
+      END
Index: LMDZ5/trunk/libf/dyn3dpar/etat0_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/etat0_netcdf.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/etat0_netcdf.F90	(revision 1520)
@@ -251,5 +251,9 @@
 !*******************************************************************************
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  if (disvert_type.eq.1) then
+    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
+  else ! we assume that we are in the disvert_type==2 case
+    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
+  endif
   pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
 !  WRITE(lunout,*) 'P3D :', p3d(10,20,:)
Index: LMDZ5/trunk/libf/dyn3dpar/exner_hyb.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_hyb.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_hyb.F	(revision 1520)
@@ -51,17 +51,32 @@
       REAL SSUM
 c
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_hyb"
+      
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
 
       if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
-        endif
         
         ! Compute pks(:),pk(:),pkf(:)
@@ -77,6 +92,8 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
+!!!! General case:
      
       unpl2k    = 1.+ 2.* kappa
Index: LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p.F	(revision 1520)
@@ -53,20 +53,38 @@
       EXTERNAL SSUM
       INTEGER ije,ijb,jje,jjb
-c
-c$OMP BARRIER
-
-      if (llm.eq.1) then
-        ! Specific behaviour for Shallow Water (1 vertical layer) case
-      
-        ! Sanity checks
-        if (kappa.ne.1) then
-          call abort_gcm("exner_hyb",
-     &    "kappa!=1 , but running in Shallow Water mode!!",42)
-        endif
-        if (cpp.ne.r) then
-        call abort_gcm("exner_hyb",
-     &    "cpp!=r , but running in Shallow Water mode!!",42)
+      logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall) 
+      character(len=*),parameter :: modname="exner_hyb_p"
+c
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
         endif
         
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+c$OMP BARRIER
+
+! Specific behaviour for Shallow Water (1 vertical layer) case
+      if (llm.eq.1) then
+      
         ! Compute pks(:),pk(:),pkf(:)
         ijb=ij_begin
@@ -116,4 +134,5 @@
       endif ! of if (llm.eq.1)
 
+!!!! General case:
 
       unpl2k    = 1.+ 2.* kappa
Index: LMDZ5/trunk/libf/dyn3dpar/exner_milieu.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_milieu.F	(revision 1520)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_milieu.F	(revision 1520)
@@ -0,0 +1,151 @@
+!
+! $Id $
+!
+      SUBROUTINE  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
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "comvert.h"
+#include "serre.h"
+
+      INTEGER  ngrid
+      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
+      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
+
+c    .... variables locales   ...
+
+      INTEGER l, ij
+      REAL dum1
+
+      REAL ppn(iim),pps(iim)
+      REAL xpn, xps
+      REAL SSUM
+      EXTERNAL SSUM
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_milieu"
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.2) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==2",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+
+!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
+      if (llm.eq.1) then
+      
+        ! Compute pks(:),pk(:),pkf(:)
+        
+        DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij) 
+          pk(ij,1) = .5*pks(ij)
+        ENDDO
+        
+        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+        ! our work is done, exit routine
+        return
+
+      endif ! of if (llm.eq.1)
+
+!!!! General case:
+
+c     -------------
+c     Calcul de pks
+c     -------------
+   
+      DO   ij  = 1, ngrid
+        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+      ENDDO
+
+      DO  ij   = 1, iim
+        ppn(ij) = aire(   ij   ) * pks(  ij     )
+        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
+      ENDDO
+      xpn      = SSUM(iim,ppn,1) /apoln
+      xps      = SSUM(iim,pps,1) /apols
+
+      DO ij   = 1, iip1
+        pks(   ij     )  =  xpn
+        pks( ij+ip1jm )  =  xps
+      ENDDO
+c
+c
+c    .... Calcul de pk  pour la couche l 
+c    --------------------------------------------
+c
+      dum1 = cpp * (2*preff)**(-kappa) 
+      DO l = 1, llm-1
+        DO   ij   = 1, ngrid
+         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+        ENDDO
+      ENDDO
+
+c    .... Calcul de pk  pour la couche l = llm ..
+c    (on met la meme distance (en log pression)  entre Pk(llm)
+c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+      DO   ij   = 1, ngrid
+         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+      ENDDO
+
+
+c    calcul de pkf
+c    -------------
+      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+      
+c    EST-CE UTILE ?? : calcul de beta
+c    --------------------------------
+      DO l = 2, llm
+        DO   ij   = 1, ngrid
+          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
+        ENDDO
+      ENDDO
+
+      RETURN
+      END
Index: LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p.F	(revision 1520)
+++ LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p.F	(revision 1520)
@@ -0,0 +1,224 @@
+!
+! $Id $
+!
+      SUBROUTINE  exner_milieu_p ( 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
+      USE parallel
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "comvert.h"
+#include "serre.h"
+
+      INTEGER  ngrid
+      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
+      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
+
+c    .... variables locales   ...
+
+      INTEGER l, ij
+      REAL dum1
+
+      REAL ppn(iim),pps(iim)
+      REAL xpn, xps
+      REAL SSUM
+      EXTERNAL SSUM
+      INTEGER ije,ijb,jje,jjb
+      logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall) 
+      character(len=*),parameter :: modname="exner_milieu_p"
+
+      ! Sanity check
+      if (firstcall) then
+        ! check that vertical discretization is compatible
+        ! with this routine
+        if (disvert_type.ne.1) then
+          call abort_gcm(modname,
+     &     "this routine should only be called if disvert_type==1",42)
+        endif
+        
+        ! sanity checks for Shallow Water case (1 vertical layer)
+        if (llm.eq.1) then
+          if (kappa.ne.1) then
+            call abort_gcm(modname,
+     &      "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+            call abort_gcm(modname,
+     &      "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+        endif ! of if (llm.eq.1)
+
+        firstcall=.false.
+      endif ! of if (firstcall)
+      
+c$OMP BARRIER
+
+! Specific behaviour for Shallow Water (1 vertical layer) case
+      if (llm.eq.1) then
+              
+        ! Compute pks(:),pk(:),pkf(:)
+        ijb=ij_begin
+        ije=ij_end
+!$OMP DO SCHEDULE(STATIC)
+        DO ij=ijb, ije
+          pks(ij)=(cpp/preff)*ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          pkf(ij,1)=pk(ij,1)
+        ENDDO
+!$OMP ENDDO
+
+!$OMP MASTER
+      if (pole_nord) then
+        DO  ij   = 1, iim
+          ppn(ij) = aire(   ij   ) * pks(  ij     )
+        ENDDO
+        xpn      = SSUM(iim,ppn,1) /apoln
+  
+        DO ij   = 1, iip1
+          pks(   ij     )  =  xpn
+          pk(ij,1) = .5*pks(ij)
+          pkf(ij,1)=pk(ij,1)
+        ENDDO
+      endif
+      
+      if (pole_sud) then
+        DO  ij   = 1, iim
+          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
+        ENDDO
+        xps      = SSUM(iim,pps,1) /apols 
+  
+        DO ij   = 1, iip1
+          pks( ij+ip1jm )  =  xps
+          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
+          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
+        ENDDO
+      endif
+!$OMP END MASTER
+
+        jjb=jj_begin
+        jje=jj_end
+        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+
+        ! our work is done, exit routine
+        return
+      endif ! of if (llm.eq.1)
+
+!!!! General case:
+
+c     -------------
+c     Calcul de pks
+c     -------------
+   
+      ijb=ij_begin
+      ije=ij_end
+
+c$OMP DO SCHEDULE(STATIC)
+      DO   ij  = ijb, ije
+        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+      ENDDO
+c$OMP ENDDO
+c Synchro OPENMP ici
+
+c$OMP MASTER
+      if (pole_nord) then
+        DO  ij   = 1, iim
+          ppn(ij) = aire(   ij   ) * pks(  ij     )
+        ENDDO
+        xpn      = SSUM(iim,ppn,1) /apoln
+  
+        DO ij   = 1, iip1
+          pks(   ij     )  =  xpn
+        ENDDO
+      endif
+      
+      if (pole_sud) then
+        DO  ij   = 1, iim
+          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
+        ENDDO
+        xps      = SSUM(iim,pps,1) /apols 
+  
+        DO ij   = 1, iip1
+          pks( ij+ip1jm )  =  xps
+        ENDDO
+      endif
+c$OMP END MASTER
+c
+c
+c    .... Calcul de pk  pour la couche l 
+c    --------------------------------------------
+c
+      dum1 = cpp * (2*preff)**(-kappa) 
+      DO l = 1, llm-1
+c$OMP DO SCHEDULE(STATIC)
+        DO   ij   = ijb, ije
+         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+        ENDDO
+c$OMP ENDDO NOWAIT        
+      ENDDO
+
+c    .... Calcul de pk  pour la couche l = llm ..
+c    (on met la meme distance (en log pression)  entre Pk(llm)
+c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+c$OMP DO SCHEDULE(STATIC)
+      DO   ij   = ijb, ije
+         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+      ENDDO
+c$OMP ENDDO NOWAIT        
+
+
+c    calcul de pkf
+c    -------------
+c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+      DO l = 1, llm
+c$OMP DO SCHEDULE(STATIC)
+         DO   ij   = ijb, ije
+           pkf(ij,l)=pk(ij,l)
+         ENDDO
+c$OMP ENDDO NOWAIT             
+      ENDDO
+
+c$OMP BARRIER
+      
+      jjb=jj_begin
+      jje=jj_end
+      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+      
+c    EST-CE UTILE ?? : calcul de beta
+c    --------------------------------
+      DO l = 2, llm
+c$OMP DO SCHEDULE(STATIC)
+        DO   ij   = ijb, ije
+          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
+        ENDDO
+c$OMP ENDDO NOWAIT             
+      ENDDO
+
+      RETURN
+      END
Index: LMDZ5/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/gcm.F	(revision 1520)
@@ -538,5 +538,5 @@
 c       write(78,*) 'q',q
 
-c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
+c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
       CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
      .              time_0)
Index: LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90	(revision 1520)
@@ -455,6 +455,10 @@
 !       Calcul niveaux pression milieu de couches 
 	CALL pression_p( ip1jmp1, ap, bp, ps, p )
-	CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
-	unskap=1./kappa
+	if (disvert_type==1) then
+          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+	else
+          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        endif
+        unskap=1./kappa
 	DO l = 1, llm
 	    DO j=jjb_u,jje_u
@@ -751,5 +755,9 @@
     ELSE
 	CALL pression_p( ip1jmp1, ap, bp, psi, p )
-	CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+	if (disvert_type==1) then
+          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+        else ! we assume that we are in the disvert_type==2 case
+          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+        endif
 	unskap=1./kappa
 	DO l = 1, llm
Index: LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90	(revision 1520)
@@ -71,4 +71,7 @@
 
   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  
+  character(len=*),parameter :: modname="iniacademic"
+  character(len=80) :: abort_message
 
   !-----------------------------------------------------------------------
@@ -210,5 +213,13 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        if (disvert_type.eq.1) then
+          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        elseif (disvert_type.eq.2) then
+          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        else
+          write(abort_message,*) "Wrong value for disvert_type: ", &
+                              disvert_type
+          call abort_gcm(modname,abort_message,0)
+        endif
         CALL massdair(p,masse)
 
Index: LMDZ5/trunk/libf/dyn3dpar/iniconst.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/iniconst.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/iniconst.F	(revision 1520)
@@ -5,4 +5,10 @@
 
       USE control_mod
+#ifdef CPP_IOIPSL
+      use IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      use ioipsl_getincom
+#endif
 
       IMPLICIT NONE
@@ -21,5 +27,6 @@
 #include "iniprint.h"
 
-
+      character(len=*),parameter :: modname="iniconst"
+      character(len=80) :: abort_message
 c
 c
@@ -48,11 +55,30 @@
       r       = cpp * kappa
 
-      write(lunout,*)'iniconst: R  CP  Kappa ',  r , cpp,  kappa
+      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
 c
 c-----------------------------------------------------------------------
 
-       CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-c
-c
-       RETURN
-       END
+! vertical discretization: default behavior depends on planet_type flag
+      if (planet_type=="earth") then
+        disvert_type=1
+      else
+        disvert_type=2
+      endif
+      ! but user can also specify using one or the other in run.def:
+      call getin('disvert_type',disvert_type)
+      write(lunout,*) trim(modname),': disvert_type=',disvert_type
+      
+      if (disvert_type==1) then
+       ! standard case for Earth (automatic generation of levels)
+       call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     &              scaleheight)
+      else if (disvert_type==2) then
+        ! standard case for planets (levels generated using z2sig.def file)
+        call disvert_noterre
+      else
+        write(abort_message,*) "Wrong value for disvert_type: ",
+     &                        disvert_type
+        call abort_gcm(modname,abort_message,0)
+      endif
+
+      END
Index: LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1520)
@@ -164,5 +164,5 @@
 
       character*80 dynhist_file, dynhistave_file
-      character*20 modname
+      character(len=*),parameter :: modname="leapfrog"
       character*80 abort_message
 
@@ -208,5 +208,4 @@
       itaufin   = nday*day_step
       itaufinp1 = itaufin +1
-      modname="leapfrog_p"
 
       itau = 0
@@ -236,5 +235,9 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+      if (disvert_type==1) then
+        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+      else ! we assume that we are in the disvert_type==2 case
+        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+      endif
 c$OMP END MASTER
 c-----------------------------------------------------------------------
@@ -691,5 +694,9 @@
 
 c$OMP BARRIER
-         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         if (disvert_type==1) then
+           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+         else ! we assume that we are in the disvert_type==2 case
+           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+         endif
 c$OMP BARRIER
            jD_cur = jD_ref + day_ini - day_ref
@@ -1030,5 +1037,9 @@
         CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
 c$OMP BARRIER
-        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        if (disvert_type==1) then
+          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        else ! we assume that we are in the disvert_type==2 case
+          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        endif
 c$OMP BARRIER
 
Index: LMDZ5/trunk/libf/dyn3dpar/limy.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/limy.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/limy.F	(revision 1520)
@@ -1,5 +1,5 @@
-!
-! $Header$
-!
+c
+c $Id$
+c
       SUBROUTINE limy(s0,sy,sm,pente_max)
 c
@@ -40,5 +40,5 @@
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
+      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
       Logical extremum,first
       save first
@@ -52,5 +52,6 @@
       REAL      SSUM
       integer ismax,ismin
-      EXTERNAL  SSUM, ismin,ismax
+      EXTERNAL  SSUM, convflu,ismin,ismax
+      EXTERNAL filtreg
 
       data first/.true./
@@ -116,14 +117,14 @@
 
 c     print*,dyqv(iip1+1)
-c     apn=abs(dyq(1)/dyqv(iip1+1))
+c     appn=abs(dyq(1)/dyqv(iip1+1))
 c     print*,dyq(ip1jm+1)
 c     print*,dyqv(ip1jm-iip1+1)
-c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 c     do ij=2,iim
-c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 c     enddo
-c     apn=min(pente_max/apn,1.)
-c     aps=min(pente_max/aps,1.)
+c     appn=min(pente_max/appn,1.)
+c     apps=min(pente_max/apps,1.)
 
 
@@ -131,13 +132,13 @@
 
 c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-c    &   apn=0.
+c    &   appn=0.
 c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-c    &   aps=0.
+c    &   apps=0.
 
 c   limitation des pentes aux poles
 c     do ij=1,iip1
-c        dyq(ij)=apn*dyq(ij)
-c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+c        dyq(ij)=appn*dyq(ij)
+c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 c     enddo
 
Index: LMDZ5/trunk/libf/dyn3dpar/logic.h
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/logic.h	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/logic.h	(revision 1520)
@@ -3,19 +3,25 @@
 !
 !
-!
+! NB: keep items of different kinds in seperate common blocs to avoid
+!     "misaligned commons" issues
 !-----------------------------------------------------------------------
 ! INCLUDE 'logic.h'
 
-      COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
+      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
 
+      COMMON/logici/ iflag_phys,iflag_trac
+      
       LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
      &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
+                     ! (only used if disvert_type==2)
 
-      INTEGER iflag_phys
-!$OMP THREADPRIVATE(/logic/)
+      integer iflag_phys,iflag_trac
+!$OMP THREADPRIVATE(/logicl/)
+!$OMP THREADPRIVATE(/logici/)
 !-----------------------------------------------------------------------
Index: LMDZ5/trunk/libf/dyn3dpar/vlsplt_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/vlsplt_p.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/vlsplt_p.F	(revision 1520)
@@ -1,2 +1,6 @@
+c
+c $Id$
+c
+
       SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
 c
@@ -561,5 +565,5 @@
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
 c     REAL newq,oldmasse
       Logical extremum,first,testcpu
@@ -732,14 +736,14 @@
 C     PRINT*,dyq(1)
 C     PRINT*,dyqv(iip1+1)
-C     apn=abs(dyq(1)/dyqv(iip1+1))
+C     appn=abs(dyq(1)/dyqv(iip1+1))
 C     PRINT*,dyq(ip1jm+1)
 C     PRINT*,dyqv(ip1jm-iip1+1)
-C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 C     DO ij=2,iim
-C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 C     ENDDO
-C     apn=min(pente_max/apn,1.)
-C     aps=min(pente_max/aps,1.)
+C     appn=min(pente_max/appn,1.)
+C     apps=min(pente_max/apps,1.)
 C
 C
@@ -747,13 +751,13 @@
 C
 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-C    &   apn=0.
+C    &   appn=0.
 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-C    &   aps=0.
+C    &   apps=0.
 C
 C   limitation des pentes aux poles
 C     DO ij=1,iip1
-C        dyq(ij)=apn*dyq(ij)
-C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+C        dyq(ij)=appn*dyq(ij)
+C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 C     ENDDO
 C
Index: LMDZ5/trunk/libf/dyn3dpar/vlspltqs_p.F
===================================================================
--- LMDZ5/trunk/libf/dyn3dpar/vlspltqs_p.F	(revision 1519)
+++ LMDZ5/trunk/libf/dyn3dpar/vlspltqs_p.F	(revision 1520)
@@ -1,5 +1,5 @@
-!
-! $Header$
-!
+c
+c $Id$
+c
        SUBROUTINE vlspltqs_p ( q,pente_max,masse,w,pbaru,pbarv,pdt,
      ,                                  p,pk,teta                 )
@@ -772,14 +772,14 @@
 C     PRINT*,dyq(1)
 C     PRINT*,dyqv(iip1+1)
-C     apn=abs(dyq(1)/dyqv(iip1+1))
+C     appn=abs(dyq(1)/dyqv(iip1+1))
 C     PRINT*,dyq(ip1jm+1)
 C     PRINT*,dyqv(ip1jm-iip1+1)
-C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
 C     DO ij=2,iim
-C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
-C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
+C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
 C     ENDDO
-C     apn=min(pente_max/apn,1.)
-C     aps=min(pente_max/aps,1.)
+C     appn=min(pente_max/appn,1.)
+C     apps=min(pente_max/apps,1.)
 C
 C
@@ -787,13 +787,13 @@
 C
 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-C    &   apn=0.
+C    &   appn=0.
 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
 C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-C    &   aps=0.
+C    &   apps=0.
 C
 C   limitation des pentes aux poles
 C     DO ij=1,iip1
-C        dyq(ij)=apn*dyq(ij)
-C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
+C        dyq(ij)=appn*dyq(ij)
+C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
 C     ENDDO
 C
