Index: trunk/libf/dyn3d/comvert.h
===================================================================
--- trunk/libf/dyn3d/comvert.h	(revision 124)
+++ trunk/libf/dyn3d/comvert.h	(revision 127)
@@ -1,14 +1,32 @@
 !
-! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: comvert.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !-----------------------------------------------------------------------
 !   INCLUDE 'comvert.h'
 
-      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
+      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,aps,bps
+      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: trunk/libf/dyn3d/disvert.F90
===================================================================
--- trunk/libf/dyn3d/disvert.F90	(revision 127)
+++ trunk/libf/dyn3d/disvert.F90	(revision 127)
@@ -0,0 +1,150 @@
+! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
+
+SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
+
+  ! Auteur : P. Le Van
+
+  IMPLICIT NONE
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  ! s = sigma ** kappa : coordonnee verticale
+  ! dsig(l) : epaisseur de la couche l ds la coord. s
+  ! sig(l) : sigma a l'interface des couches l et l-1
+  ! ds(l) : distance entre les couches l et l-1 en coord.s
+
+  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)
+  real zk, zkm1, dzk1, dzk2, k0, k1
+
+  INTEGER l
+  REAL dsigmin
+  REAL alpha, beta, deltaz,scaleheight
+  INTEGER iostat
+  REAL x
+  character(len=*),parameter :: modname="disvert"
+
+  !-----------------------------------------------------------------------
+
+  ! default scaleheight is 8km for earth
+  scaleheight=8.
+
+  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
+
+  IF (iostat == 0) THEN
+     ! cas 1 on lit les options dans sigma.def:
+     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
+     READ(99, *) k0 ! nombre de couches dans la transition surf
+     READ(99, *) k1 ! nombre de couches dans la transition haute
+     CLOSE(99)
+     alpha=deltaz/(llm*scaleheight)
+     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+                               scaleheight, alpha, k0, k1, beta
+
+     alpha=deltaz/tanh(1./k0)*2.
+     zkm1=0.
+     sig(1)=1.
+     do l=1, llm
+        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)
+        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
+        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
+        zkm1=zk
+     enddo
+
+     sig(llm+1)=0.
+
+     DO l = 1, llm
+        dsig(l) = sig(l)-sig(l+1)
+     end DO
+  ELSE
+     if (ok_strato) then
+        if (llm==39) then
+           dsigmin=0.3
+        else if (llm==50) then
+           dsigmin=1.
+        else
+           write(lunout,*) trim(modname), &
+           ' ATTENTION discretisation z a ajuster'
+           dsigmin=1.
+        endif
+        write(lunout,*) trim(modname), &
+        ' Discretisation verticale DSIGMIN=',dsigmin
+     endif
+
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+
+        IF (ok_strato) THEN
+           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
+        ELSE
+           dsig(l) = 1.0 + 7.0 * SIN(x)**2
+        ENDIF
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig(l) = sig(l+1) + dsig(l)
+     ENDDO
+  ENDIF
+
+  DO l=1, llm
+     nivsigs(l) = REAL(l)
+  ENDDO
+
+  DO l=1, llmp1
+     nivsig(l)= REAL(l)
+  ENDDO
+
+  ! .... Calculs de ap(l) et de bp(l) ....
+  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
+
+  bp(llmp1) = 0.
+
+  DO l = 1, llm
+     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+     ap(l) = pa * ( sig(l) - bp(l) )
+  ENDDO
+
+  bp(1)=1.
+  ap(1)=0.
+
+  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
+
+  write(lunout, *)  trim(modname),': BP '
+  write(lunout, *) bp
+  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, *) 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))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
+          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
+  ENDDO
+
+  write(lunout, *) trim(modname),': PRESNIVS '
+  write(lunout, *) presnivs
+
+END SUBROUTINE disvert
Index: trunk/libf/dyn3d/disvert_noterre.F
===================================================================
--- trunk/libf/dyn3d/disvert_noterre.F	(revision 124)
+++ trunk/libf/dyn3d/disvert_noterre.F	(revision 127)
@@ -1,2 +1,3 @@
+! $Id: $
       SUBROUTINE disvert_noterre
 
@@ -22,5 +23,5 @@
 c
 c=======================================================================
-c    Discretisation verticale en coordonnée hybride 
+c    Discretisation verticale en coordonnée hybride (ou sigma)
 c
 c=======================================================================
@@ -49,4 +50,5 @@
       integer iz 
       real z, ps,p
+      character(len=*),parameter :: modname="disvert_noterre"
 
 c
@@ -54,9 +56,9 @@
 c
 ! Initializations:
-      pi=2.*ASIN(1.)
+!      pi=2.*ASIN(1.) ! already done in iniconst
       
       hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
       CALL getin('hybrid',hybrid)
-      write(lunout,*)'disvert_noterre: hybrid=',hybrid
+      write(lunout,*) trim(modname),': hybrid=',hybrid
 
 ! Ouverture possible de fichiers typiquement E.T.
@@ -156,9 +158,9 @@
 
       DO l=1,llm
-        nivsigs(l) = FLOAT(l)
+        nivsigs(l) = REAL(l)
       ENDDO
 
       DO l=1,llmp1
-        nivsig(l)= FLOAT(l)
+        nivsig(l)= REAL(l)
       ENDDO
 
@@ -199,7 +201,7 @@
       bp(llmp1) =   0.
 
-      write(lunout,*)' BP '
+      write(lunout,*) trim(modname),': BP '
       write(lunout,*)  bp
-      write(lunout,*)' AP '
+      write(lunout,*) trim(modname),': AP '
       write(lunout,*)  ap
 
@@ -225,7 +227,7 @@
       end if
 
-      write(lunout,*)' BPs '
+      write(lunout,*) trim(modname),': BPs '
       write(lunout,*)  bps
-      write(lunout,*)' APs'
+      write(lunout,*) trim(modname),': APs'
       write(lunout,*)  aps
 
@@ -235,5 +237,5 @@
       ENDDO
 
-      write(lunout,*)' PRESNIVS' 
+      write(lunout,*)trim(modname),' : PRESNIVS' 
       write(lunout,*)presnivs 
       write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
Index: trunk/libf/dyn3d/disvert_terre.F90
===================================================================
--- trunk/libf/dyn3d/disvert_terre.F90	(revision 124)
+++ 	(revision )
@@ -1,142 +1,0 @@
-! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
-
-SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
-
-  ! Auteur : P. Le Van
-
-  IMPLICIT NONE
-
-  include "dimensions.h"
-  include "paramet.h"
-  include "iniprint.h"
-  include "logic.h"
-
-  ! s = sigma ** kappa : coordonnee verticale
-  ! dsig(l) : epaisseur de la couche l ds la coord. s
-  ! sig(l) : sigma a l'interface des couches l et l-1
-  ! 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 sig(llm+1), dsig(llm)
-  real zk, zkm1, dzk1, dzk2, k0, k1
-
-  INTEGER l
-  REAL dsigmin
-  REAL alpha, beta, deltaz, h
-  INTEGER iostat
-  REAL pi, x
-
-  !-----------------------------------------------------------------------
-
-  pi = 2 * ASIN(1.)
-
-  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
-
-  IF (iostat == 0) THEN
-     ! cas 1 on lit les options dans sigma.def:
-     READ(99, *) h ! hauteur d'echelle 8.
-     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
-     READ(99, *) beta ! facteur d'acroissement en haut 1.3
-     READ(99, *) k0 ! nombre de couches dans la transition surf
-     READ(99, *) k1 ! nombre de couches dans la transition haute
-     CLOSE(99)
-     alpha=deltaz/(llm*h)
-     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
-
-     alpha=deltaz/tanh(1./k0)*2.
-     zkm1=0.
-     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))
-
-        dzk1=alpha*tanh(l/k0)
-        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
-        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
-        zkm1=zk
-     enddo
-
-     sig(llm+1)=0.
-
-     DO l = 1, llm
-        dsig(l) = sig(l)-sig(l+1)
-     end DO
-  ELSE
-     if (ok_strato) then
-        if (llm==39) then
-           dsigmin=0.3
-        else if (llm==50) then
-           dsigmin=1.
-        else
-           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
-           dsigmin=1.
-        endif
-        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
-     endif
-
-     DO l = 1, llm
-        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
-
-        IF (ok_strato) THEN
-           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
-                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
-        ELSE
-           dsig(l) = 1.0 + 7.0 * SIN(x)**2
-        ENDIF
-     ENDDO
-     dsig = dsig / sum(dsig)
-     sig(llm+1) = 0.
-     DO l = llm, 1, -1
-        sig(l) = sig(l+1) + dsig(l)
-     ENDDO
-  ENDIF
-
-  DO l=1, llm
-     nivsigs(l) = REAL(l)
-  ENDDO
-
-  DO l=1, llmp1
-     nivsig(l)= REAL(l)
-  ENDDO
-
-  ! .... Calculs de ap(l) et de bp(l) ....
-  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
-
-  bp(llmp1) = 0.
-
-  DO l = 1, llm
-     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-     ap(l) = pa * ( sig(l) - bp(l) )
-  ENDDO
-
-  bp(1)=1.
-  ap(1)=0.
-
-  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
-
-  write(lunout, *)'disvert: BP '
-  write(lunout, *) bp
-  write(lunout, *)'disvert: 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'
-  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)/ &
-          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
-  ENDDO
-
-  write(lunout, *) 'disvert: PRESNIVS '
-  write(lunout, *) presnivs
-
-END SUBROUTINE disvert_terre
Index: trunk/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- trunk/libf/dyn3d/etat0_netcdf.F90	(revision 124)
+++ trunk/libf/dyn3d/etat0_netcdf.F90	(revision 127)
@@ -1,8 +1,8 @@
 !
-! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
+! $Id: etat0_netcdf.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE etat0_netcdf(ib, masque, letat0)
+SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
 !
 !-------------------------------------------------------------------------------
@@ -37,4 +37,5 @@
   LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
   REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
+  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
   LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
 #ifndef CPP_EARTH
@@ -51,5 +52,5 @@
   REAL,    DIMENSION(klon)                 :: tsol, qsol
   REAL,    DIMENSION(klon)                 :: sn, rugmer, run_off_lic_0
-  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol, phis
+  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol
   REAL,    DIMENSION(iip1,jjp1,llm+1)      :: p3d
   REAL,    DIMENSION(iip1,jjp1,llm)        :: uvent, t3d, tpot, qsat, qd
@@ -138,7 +139,5 @@
                    flag_aerosol, new_aod,                               &
                    bl95_b0, bl95_b1,                                    &
-                   iflag_thermals,nsplit_thermals,tau_thermals,         &
-                   iflag_thermals_ed,iflag_thermals_optflux,            &
-                   iflag_coupl,iflag_clos,iflag_wake, read_climoz,      &
+                   read_climoz,                                         &
                    alp_offset)
 
@@ -252,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: trunk/libf/dyn3d/exner_hyb.F
===================================================================
--- trunk/libf/dyn3d/exner_hyb.F	(revision 124)
+++ trunk/libf/dyn3d/exner_hyb.F	(revision 127)
@@ -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: trunk/libf/dyn3d/exner_milieu.F
===================================================================
--- trunk/libf/dyn3d/exner_milieu.F	(revision 124)
+++ trunk/libf/dyn3d/exner_milieu.F	(revision 127)
@@ -48,18 +48,34 @@
       REAL SSUM
       EXTERNAL SSUM
+      logical,save :: firstcall=.true.
+      character(len=*),parameter :: modname="exner_milieu"
 
-      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)
+      ! 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(:)
         
@@ -74,7 +90,9 @@
         ! our work is done, exit routine
         return
+
       endif ! of if (llm.eq.1)
 
-      
+!!!! General case:
+
 c     -------------
 c     Calcul de pks
Index: trunk/libf/dyn3d/guide_mod.F90
===================================================================
--- trunk/libf/dyn3d/guide_mod.F90	(revision 124)
+++ trunk/libf/dyn3d/guide_mod.F90	(revision 127)
@@ -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: trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- trunk/libf/dyn3d/iniacademic.F90	(revision 124)
+++ trunk/libf/dyn3d/iniacademic.F90	(revision 127)
@@ -1,4 +1,4 @@
 !
-! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
+! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -71,4 +71,7 @@
 
   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+  
+  character(len=*),parameter :: modname="iniacademic"
+  character(len=80) :: abort_message
 
   !-----------------------------------------------------------------------
@@ -76,7 +79,4 @@
   ! --------------------------------------
   !
-  
-  write(lunout,*) 'Iniacademic'
-  
   ! initialize planet radius, rotation rate,...
   call conf_planete
@@ -136,6 +136,4 @@
      teta0=315.     ! mean Teta (S.H. 315K)
      CALL getin('teta0',teta0)
-     write(lunout,*) 'Iniacademic - teta0 ',teta0
-     write(lunout,*) 'Iniacademic - rad ',rad
      ttp=200.       ! Tropopause temperature (S.H. 200K)
      CALL getin('ttp',ttp)
@@ -183,5 +181,4 @@
            tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
                 -delt_z*(1.-ddsin*ddsin)*log(zsig)
-           !! Aymeric -- tests particuliers
            if (planet_type=="giant") then
              tetajl(j,l)=teta0+(delt_y*                   &
@@ -207,5 +204,4 @@
         enddo
      enddo
-!     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
 
      ! 3. Initialize fields (if necessary)
@@ -217,8 +213,12 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        if (planet_type=="earth") then
-          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
-          call exner_milieu(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          write(abort_message,*) "Wrong value for disvert_type: ", &
+                              disvert_type
+          call abort_gcm(modname,abort_message,0)
         endif
         CALL massdair(p,masse)
Index: trunk/libf/dyn3d/iniconst.F
===================================================================
--- trunk/libf/dyn3d/iniconst.F	(revision 124)
+++ trunk/libf/dyn3d/iniconst.F	(revision 127)
@@ -1,8 +1,14 @@
 !
-! $Id: iniconst.F 1403 2010-07-01 09:02:53Z fairhead $
+! $Id: iniconst.F 1520 2011-05-23 11:37:09Z emillour $
 !
       SUBROUTINE iniconst
 
       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
@@ -22,4 +28,6 @@
 
 
+      character(len=*),parameter :: modname="iniconst"
+      character(len=80) :: abort_message
 c
 c
@@ -49,14 +57,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-----------------------------------------------------------------------
 
-      if (planet_type.eq."earth") then
-       CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+! vertical discretization: default behavior depends on planet_type flag
+      if (planet_type=="earth") then
+        disvert_type=1
       else
-       CALL disvert_noterre
+        disvert_type=2
       endif
-c
-       RETURN
-       END
+      ! 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: trunk/libf/dyn3d/leapfrog.F
===================================================================
--- trunk/libf/dyn3d/leapfrog.F	(revision 124)
+++ trunk/libf/dyn3d/leapfrog.F	(revision 127)
@@ -178,5 +178,5 @@
 
       character*80 dynhist_file, dynhistave_file
-      character(len=20) :: modname
+      character(len=*),parameter :: modname="leapfrog"
       character*80 abort_message
 
@@ -206,5 +206,4 @@
       endif
       itaufinp1 = itaufin +1
-      modname="leapfrog"
       
 c INITIALISATIONS
@@ -238,7 +237,7 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      if (planet_type.eq."earth") then
+      if (disvert_type==1) then
         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      else
+      else ! we assume that we are in the disvert_type==2 case
         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
       endif
@@ -409,7 +408,7 @@
 
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         if (planet_type.eq."earth") then
-           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-         else
+         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
@@ -531,7 +530,7 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        if (planet_type.eq."earth") then
+        if (disvert_type==1) then
           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
+        else ! we assume that we are in the disvert_type==2 case
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
Index: trunk/libf/dyn3d/logic.h
===================================================================
--- trunk/libf/dyn3d/logic.h	(revision 124)
+++ trunk/libf/dyn3d/logic.h	(revision 127)
@@ -1,4 +1,4 @@
 !
-! $Id: logic.h 1319 2010-02-23 21:29:54Z fairhead $
+! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
 !
 !
@@ -11,5 +11,5 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0,hybrid
+     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid
 
       COMMON/logici/ iflag_phys,iflag_trac
@@ -18,7 +18,8 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
-     &  ,ok_limit,ok_etat0
-      logical hybrid ! vertcal coordinate is hybrid if true (sigma otherwise)
+     &  ,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,iflag_trac
+      integer iflag_phys,iflag_trac
 !-----------------------------------------------------------------------
