Index: LMDZ5/branches/testing/libf/dyn3d/calfis.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/calfis.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/calfis.F	(revision 1665)
@@ -165,5 +165,5 @@
       PARAMETER(ntetaSTD=3)
       REAL rtetaSTD(ntetaSTD)
-      DATA rtetaSTD/350., 380., 405./
+      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
       REAL PVteta(ngridmx,ntetaSTD)
 c
@@ -434,5 +434,6 @@
 c
       if (planet_type=="earth") then
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
+! PVtheta calls tetalevel, which is in the physics
 cIM calcul PV a teta=350, 380, 405K
       CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
@@ -450,6 +451,4 @@
 
 
-      if (planet_type=="earth") then
-#ifdef CPP_EARTH
 
 !      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
@@ -460,5 +459,8 @@
       zdqfic(:,:,:)=0.
 
-      do isplit=1,nsplit_phys
+      if (planet_type=="earth") then
+#ifdef CPP_PHYS
+
+       do isplit=1,nsplit_phys
 
          jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
@@ -503,5 +505,10 @@
          zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
 
-      enddo
+       enddo ! of do isplit=1,nsplit_phys
+
+#endif
+! of #ifdef CPP_PHYS
+      endif ! of if (planet_type=="earth")
+
       zdufi(:,:)=zdufic(:,:)/nsplit_phys
       zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
@@ -509,6 +516,4 @@
       zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
 
-#endif
-      endif !of if (planet_type=="earth")
 
 500   CONTINUE
Index: LMDZ5/branches/testing/libf/dyn3d/ce0l.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/ce0l.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/ce0l.F90	(revision 1665)
@@ -28,4 +28,5 @@
   IMPLICIT NONE
 #ifndef CPP_EARTH
+#include "iniprint.h"
   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
 #else
Index: LMDZ5/branches/testing/libf/dyn3d/comvert.h
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/comvert.h	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/comvert.h	(revision 1665)
@@ -9,5 +9,5 @@
      &               aps(llm),bps(llm),scaleheight
 
-      common/comverti/disvert_type
+      common/comverti/disvert_type, pressure_exner
 
       real ap     ! hybrid pressure contribution at interlayers
@@ -30,3 +30,7 @@
                            !     using 'z2sig.def' (or 'esasig.def) file
 
+      logical pressure_exner
+!     compute pressure inside layers using Exner function, else use mean
+!     of pressure values at interfaces
+
  !-----------------------------------------------------------------------
Index: LMDZ5/branches/testing/libf/dyn3d/conf_gcm.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/conf_gcm.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/conf_gcm.F	(revision 1665)
@@ -155,4 +155,12 @@
       nday = 10
       CALL getin('nday',nday)
+
+!Config  Key  = starttime
+!Config  Desc = Heure de depart de la simulation
+!Config  Def  = 0
+!Config  Help = Heure de depart de la simulation
+!Config         en jour
+      starttime = 0
+      CALL getin('starttime',starttime)
 
 !Config  Key  = day_step
Index: LMDZ5/branches/testing/libf/dyn3d/control_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/control_mod.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/control_mod.F90	(revision 1665)
@@ -10,5 +10,5 @@
   IMPLICIT NONE
 
-  REAL    :: periodav
+  REAL    :: periodav, starttime
   INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
   INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy
Index: LMDZ5/branches/testing/libf/dyn3d/disvert.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/disvert.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/disvert.F90	(revision 1665)
@@ -4,4 +4,7 @@
 
   ! Auteur : P. Le Van
+
+  use new_unit_m, only: new_unit
+  use ioipsl, only: getin
 
   IMPLICIT NONE
@@ -26,20 +29,28 @@
   real zk, zkm1, dzk1, dzk2, k0, k1
 
-  INTEGER l
+  INTEGER l, unit
   REAL dsigmin
   REAL alpha, beta, deltaz
-  INTEGER iostat
   REAL x
   character(len=*),parameter :: modname="disvert"
 
+  character(len=6):: vert_sampling
+  ! (allowed values are "param", "tropo", "strato" and "read")
+
   !-----------------------------------------------------------------------
+
+  print *, "Call sequence information: disvert"
 
   ! default scaleheight is 8km for earth
   scaleheight=8.
 
-  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
+  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
+  call getin('vert_sampling', vert_sampling)
+  print *, 'vert_sampling = ' // vert_sampling
 
-  IF (iostat == 0) THEN
-     ! cas 1 on lit les options dans sigma.def:
+  select case (vert_sampling)
+  case ("param")
+     ! On lit les options dans sigma.def:
+     OPEN(99, file='sigma.def', status='old', form='formatted')
      READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
@@ -68,32 +79,8 @@
 
      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
-
+  case("tropo")
      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
+        dsig(l) = 1.0 + 7.0 * SIN(x)**2
      ENDDO
      dsig = dsig / sum(dsig)
@@ -102,5 +89,37 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
-  ENDIF
+  case("strato")
+     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
+
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig(l) = sig(l+1) + dsig(l)
+     ENDDO
+  case("read")
+     call new_unit(unit)
+     open(unit, file="hybrid.txt", status="old", action="read", &
+          position="rewind")
+     read(unit, fmt=*) ! skip title line
+     do l = 1, llm + 1
+        read(unit, fmt=*) sig(l)
+     end do
+     close(unit)
+  case default
+     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
+  END select
 
   DO l=1, llm
Index: LMDZ5/branches/testing/libf/dyn3d/dynetat0.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/dynetat0.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/dynetat0.F	(revision 1665)
@@ -119,4 +119,5 @@
       day_ini = tab_cntrl(30)
       itau_dyn = tab_cntrl(31)
+      start_time = tab_cntrl(32)
 c   .................................................................
 c
Index: LMDZ5/branches/testing/libf/dyn3d/dynredem.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/dynredem.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/dynredem.F	(revision 1665)
@@ -120,4 +120,6 @@
        tab_cntrl(30) = REAL(iday_end)
        tab_cntrl(31) = REAL(itau_dyn + itaufin)
+c start_time: start_time of simulation (not necessarily 0.)
+       tab_cntrl(32) = start_time
 c
 c    .........................................................
Index: LMDZ5/branches/testing/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/etat0_netcdf.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/etat0_netcdf.F90	(revision 1665)
@@ -251,7 +251,7 @@
 !*******************************************************************************
   CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  if (disvert_type.eq.1) then
+  if (pressure_exner) then
     CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
-  else ! we assume that we are in the disvert_type==2 case
+  else
     CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
   endif
Index: LMDZ5/branches/testing/libf/dyn3d/exner_hyb.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/exner_hyb.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/exner_hyb.F	(revision 1665)
@@ -56,11 +56,4 @@
       ! 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
Index: LMDZ5/branches/testing/libf/dyn3d/exner_milieu.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/exner_milieu.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/exner_milieu.F	(revision 1665)
@@ -53,11 +53,4 @@
       ! 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
Index: LMDZ5/branches/testing/libf/dyn3d/gcm.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/gcm.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/gcm.F	(revision 1665)
@@ -21,6 +21,5 @@
 ! A nettoyer. On ne veut qu'une ou deux routines d'interface 
 ! dynamique -> physique pour l'initialisation
-! Ehouarn: for now these only apply to Earth:
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
       USE dimphy
       USE comgeomphy
@@ -88,5 +87,4 @@
 
       REAL zdtvr
-      INTEGER nbetatmoy, nbetatdem,nbetat
 
 c   variables dynamiques
@@ -181,11 +179,8 @@
 ! A nettoyer. On ne veut qu'une ou deux routines d'interface 
 ! dynamique -> physique pour l'initialisation
-! Ehouarn : temporarily (?) keep this only for Earth
-      if (planet_type.eq."earth") then
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
       CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
       call InitComgeomphy
 #endif
-      endif
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 c-----------------------------------------------------------------------
@@ -305,4 +300,15 @@
 C on remet le calendrier à zero si demande
 c
+      IF (start_time /= starttime) then
+        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
+     &,' fichier restart ne correspond pas à celle lue dans le run.def'
+        IF (raz_date == 1) then
+          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
+          start_time = starttime
+        ELSE
+          WRITE(lunout,*)'Je m''arrete'
+          CALL abort
+        ENDIF
+      ENDIF
       IF (raz_date == 1) THEN
         annee_ref = anneeref
@@ -373,7 +379,4 @@
 #endif
 
-c  nombre d'etats dans les fichiers demarrage et histoire
-      nbetatdem = nday / iecri
-      nbetatmoy = nday / periodav + 1
 
       if (iflag_phys.eq.1) then
@@ -428,14 +431,11 @@
          WRITE(lunout,*)
      .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
-! Earth:
-         if (planet_type.eq."earth") then
-#ifdef CPP_EARTH
+! Physics:
+#ifdef CPP_PHYS
          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
 #endif
-         endif ! of if (planet_type.eq."earth")
          call_iniphys=.false.
       ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
-!#endif
 
 c  numero de stockage pour les fichiers de redemarrage:
@@ -459,5 +459,5 @@
 #endif
 
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
 ! Create start file (startphy.nc) and boundary conditions (limit.nc) 
 ! for the Earth verstion
Index: LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90	(revision 1665)
@@ -644,7 +644,7 @@
 ! -----------------------------------------------------------------
     CALL pression( ip1jmp1, ap, bp, psi, p )
-    if (disvert_type==1) then
+    if (pressure_exner) then
       CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
-    else ! we assume that we are in the disvert_type==2 case
+    else
       CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
     endif
Index: LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90	(revision 1665)
@@ -222,12 +222,8 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        if (disvert_type.eq.1) then
+        if (pressure_exner) then
           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        elseif (disvert_type.eq.2) then
+        else
           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/branches/testing/libf/dyn3d/iniconst.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/iniconst.F	(revision 1664)
+++ 	(revision )
@@ -1,84 +1,0 @@
-!
-! $Id$
-!
-      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
-c
-c      P. Le Van
-c
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "temps.h"
-#include "comvert.h"
-#include "iniprint.h"
-
-      character(len=*),parameter :: modname="iniconst"
-      character(len=80) :: abort_message
-c
-c
-c
-c-----------------------------------------------------------------------
-c   dimension des boucles:
-c   ----------------------
-
-      im      = iim
-      jm      = jjm
-      lllm    = llm
-      imp1    = iim 
-      jmp1    = jjm + 1
-      lllmm1  = llm - 1
-      lllmp1  = llm + 1
-
-c-----------------------------------------------------------------------
-
-      dtphys  = iphysiq * dtvr
-      unsim   = 1./iim
-      pi      = 2.*ASIN( 1. )
-
-c-----------------------------------------------------------------------
-c
-
-      r       = cpp * kappa
-
-      write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
-c
-c-----------------------------------------------------------------------
-
-! 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/branches/testing/libf/dyn3d/iniconst.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/iniconst.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/dyn3d/iniconst.F90	(revision 1665)
@@ -0,0 +1,84 @@
+!
+! $Id$
+!
+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
+  !
+  !      P. Le Van
+  !
+  !   Declarations:
+  !   -------------
+  !
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "temps.h"
+  include "comvert.h"
+  include "iniprint.h"
+
+  character(len=*),parameter :: modname="iniconst"
+  character(len=80) :: abort_message
+  !
+  !
+  !
+  !-----------------------------------------------------------------------
+  !   dimension des boucles:
+  !   ----------------------
+
+  im      = iim
+  jm      = jjm
+  lllm    = llm
+  imp1    = iim 
+  jmp1    = jjm + 1
+  lllmm1  = llm - 1
+  lllmp1  = llm + 1
+
+  !-----------------------------------------------------------------------
+
+  dtphys  = iphysiq * dtvr
+  unsim   = 1./iim
+  pi      = 2.*ASIN( 1. )
+
+  !-----------------------------------------------------------------------
+  !
+
+  r       = cpp * kappa
+
+  write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
+  !
+  !-----------------------------------------------------------------------
+
+  ! 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
+
+  pressure_exner = disvert_type == 1 ! default value
+  call getin('pressure_exner', pressure_exner)
+
+  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 SUBROUTINE iniconst
Index: LMDZ5/branches/testing/libf/dyn3d/inidissip.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/inidissip.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/inidissip.F90	(revision 1665)
@@ -28,7 +28,8 @@
 ! Local variables:
   REAL fact,zvert(llm),zz
-  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
+  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
+  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
   REAL ullm,vllm,umin,vmin,zhmin,zhmax
-  REAL zllm,z1llm
+  REAL zllm
 
   INTEGER l,ij,idum,ii
@@ -78,16 +79,11 @@
   DO l = 1,50
      IF(lstardis) THEN
-        CALL divgrad2(1,zh,deltap,niterh,zh)
+        CALL divgrad2(1,zh,deltap,niterh,divgra)
      ELSE
-        CALL divgrad (1,zh,niterh,zh)
+        CALL divgrad (1,zh,niterh,divgra)
      ENDIF
 
-     CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-     zllm  = ABS( zhmax )
-     z1llm = 1./zllm
-     DO ij = 1,ip1jmp1
-        zh(ij) = zh(ij)* z1llm
-     ENDDO
+     zllm  = ABS(maxval(divgra))
+     zh = divgra / zllm
   ENDDO
 
@@ -123,30 +119,19 @@
            !cccc             CALL covcont( 1,zu,zv,zu,zv )
            IF(lstardis) THEN
-              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
+              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
            ELSE
-              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
+              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
            ENDIF
         ELSE
            IF(lstardis) THEN
-              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
+              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
            ELSE
-              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
+              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
            ENDIF
         ENDIF
 
-        CALL minmax(iip1*jjp1,zu,umin,ullm )
-        CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-        ullm = ABS  ( ullm )
-        vllm = ABS  ( vllm )
-
-        zllm  = MAX( ullm,vllm )
-        z1llm = 1./ zllm
-        DO ij = 1, ip1jmp1
-           zu(ij) = zu(ij)* z1llm
-        ENDDO
-        DO ij = 1, ip1jm
-           zv(ij) = zv(ij)* z1llm
-        ENDDO
+        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
+        zu = gx / zllm
+        zv = gy / zllm
      end DO
 
Index: LMDZ5/branches/testing/libf/dyn3d/inigrads.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/inigrads.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/inigrads.F	(revision 1665)
@@ -9,5 +9,5 @@
       implicit none
 
-      integer if,im,jm,lm,i,j,l,lnblnk
+      integer if,im,jm,lm,i,j,l
       real x(im),y(jm),z(lm),fx,fy,fz,dt
       real xmin,xmax,ymin,ymax
@@ -40,5 +40,5 @@
       ivar(if)=0
 
-      fichier(if)=file(1:lnblnk(file))
+      fichier(if)=trim(file)
 
       firsttime(if)=.true.
@@ -70,7 +70,7 @@
 
       print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if))
-      print*,file(1:lnblnk(file))//'.dat'
+      print*,trim(file)//'.dat'
 
-      OPEN (unit(if)+1,FILE=file(1:lnblnk(file))//'.dat'
+      OPEN (unit(if)+1,FILE=trim(file)//'.dat'
      s   ,FORM='unformatted',
      s   ACCESS='direct'
Index: LMDZ5/branches/testing/libf/dyn3d/integrd.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/integrd.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/integrd.F	(revision 1665)
@@ -4,5 +4,6 @@
       SUBROUTINE integrd
      $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
-     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold )
+     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis !,finvmaold
+     &  )
 
       use control_mod, only : planet_type
@@ -34,20 +35,30 @@
 #include "temps.h"
 #include "serre.h"
+#include "iniprint.h"
 
 c   Arguments:
 c   ----------
 
-      INTEGER nq
-
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
-      REAL q(ip1jmp1,llm,nq)
-      REAL ps(ip1jmp1),masse(ip1jmp1,llm),phis(ip1jmp1)
-
-      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
-      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1),massem1(ip1jmp1,llm)
-
-      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
-      REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm)
+      integer,intent(in) :: nq ! number of tracers to handle in this routine
+      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
+      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
+      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
+      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
+      real,intent(inout) :: ps(ip1jmp1) ! surface pressure
+      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
+      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
+      ! values at previous time step
+      real,intent(inout) :: vcovm1(ip1jm,llm)
+      real,intent(inout) :: ucovm1(ip1jmp1,llm)
+      real,intent(inout) :: tetam1(ip1jmp1,llm)
+      real,intent(inout) :: psm1(ip1jmp1)
+      real,intent(inout) :: massem1(ip1jmp1,llm)
+      ! the tendencies to add
+      real,intent(in) :: dv(ip1jm,llm)
+      real,intent(in) :: du(ip1jmp1,llm)
+      real,intent(in) :: dteta(ip1jmp1,llm)
+      real,intent(in) :: dp(ip1jmp1)
+      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
+!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
 
 c   Local:
@@ -55,5 +66,6 @@
 
       REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
-      REAL massescr( ip1jmp1,llm ), finvmasse(ip1jmp1,llm)
+      REAL massescr( ip1jmp1,llm )
+!      REAL finvmasse(ip1jmp1,llm)
       REAL p(ip1jmp1,llmp1)
       REAL tpn,tps,tppn(iim),tpps(iim)
@@ -61,5 +73,5 @@
       REAL deltap( ip1jmp1,llm )
 
-      INTEGER  l,ij,iq
+      INTEGER  l,ij,iq,i,j
 
       REAL SSUM
@@ -88,7 +100,12 @@
       DO ij = 1,ip1jmp1
         IF( ps(ij).LT.0. ) THEN
-         PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij)
-         print *, ' dans integrd'
-         stop 1
+         write(lunout,*) "integrd: negative surface pressure ",ps(ij)
+         write(lunout,*) " at node ij =", ij
+         ! since ij=j+(i-1)*jjp1 , we have
+         j=modulo(ij,jjp1)
+         i=1+(ij-j)/jjp1
+         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
+     &                   " lat = ",rlatu(j)*180./pi, " deg"
+         stop
         ENDIF
       ENDDO
@@ -110,6 +127,8 @@
       CALL massdair (     p  , masse         )
 
-      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
-      CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
+! Ehouarn : we don't use/need finvmaold and finvmasse,
+!           so might as well not compute them
+!      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
+!      CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
 c
 
@@ -218,6 +237,6 @@
        ENDDO
 
-
-      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
+! Ehouarn: forget about finvmaold
+!      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
 
       endif ! of if (planet_type.eq."earth")
Index: LMDZ5/branches/testing/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/leapfrog.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/leapfrog.F	(revision 1665)
@@ -72,5 +72,4 @@
 
       real zqmin,zqmax
-      INTEGER nbetatmoy, nbetatdem,nbetat
 
 c   variables dynamiques
@@ -118,5 +117,6 @@
 
       REAL  SSUM
-      REAL time_0 , finvmaold(ip1jmp1,llm)
+      REAL time_0
+!     REAL finvmaold(ip1jmp1,llm)
 
 cym      LOGICAL  lafin
@@ -212,7 +212,7 @@
       dq(:,:,:)=0.
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      if (disvert_type==1) then
+      if (pressure_exner) then
         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-      else ! we assume that we are in the disvert_type==2 case
+      else
         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
       endif
@@ -222,9 +222,12 @@
 c   ----------------------------------
 
-   1  CONTINUE
-
-      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec) 
-      jH_cur = jH_ref +                                                 &
-     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+   1  CONTINUE ! Matsuno Forward step begins here
+
+      jD_cur = jD_ref + day_ini - day_ref +                             &
+     &          itau/day_step
+      jH_cur = jH_ref + start_time +                                    &
+     &          mod(itau,day_step)/float(day_step) 
+      jD_cur = jD_cur + int(jH_cur)
+      jH_cur = jH_cur - int(jH_cur)
 
 
@@ -255,9 +258,9 @@
 
 c   ...    P.Le Van .26/04/94  ....
-
-      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
-      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
-
-   2  CONTINUE
+! Ehouarn: finvmaold is actually not used
+!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
+!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
+
+   2  CONTINUE ! Matsuno backward or leapfrog step begins here
 
 c-----------------------------------------------------------------------
@@ -302,4 +305,5 @@
 c   --------------------------------
 
+      ! compute geopotential phi()
       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 
@@ -315,5 +319,6 @@
 
       IF( forward. OR . leapf )  THEN
-
+! Ehouarn: NB: at this point p with ps are not synchronized
+!              (whereas mass and ps are...)
          CALL caladvtrac(q,pbaru,pbarv,
      *        p, masse, dq,  teta,
@@ -340,6 +345,6 @@
 
        CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
-     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
-     $              finvmaold                                    )
+     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
+!     $              finvmaold                                    )
 
 
@@ -364,7 +369,7 @@
 
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         if (disvert_type==1) then
+         if (pressure_exner) then
            CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
-         else ! we assume that we are in the disvert_type==2 case
+         else
            CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
          endif
@@ -372,8 +377,14 @@
 !           rdaym_ini  = itau * dtvr / daysec
 !           rdayvrai   = rdaym_ini  + day_ini
-           jD_cur = jD_ref + day_ini - day_ref
-     $        + int (itau * dtvr / daysec) 
-           jH_cur = jH_ref +                                            &
-     &              (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+!           jD_cur = jD_ref + day_ini - day_ref
+!     $        + int (itau * dtvr / daysec) 
+!           jH_cur = jH_ref +                                            &
+!     &              (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+           jD_cur = jD_ref + day_ini - day_ref +                        &
+     &          itau/day_step
+           jH_cur = jH_ref + start_time +                               &
+     &              mod(itau,day_step)/float(day_step) 
+           jD_cur = jD_cur + int(jH_cur)
+           jH_cur = jH_cur - int(jH_cur)
 !         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
 !         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
@@ -394,6 +405,8 @@
 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
            IF (planet_type.eq."earth") THEN
+#ifdef CPP_EARTH
             CALL diagedyn(ztit,2,1,1,dtphys
      &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
+#endif
            ENDIF
          ENDIF ! of IF (ip_ebil_dyn.ge.1 )
@@ -472,7 +485,7 @@
 
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        if (disvert_type==1) then
+        if (pressure_exner) then
           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else ! we assume that we are in the disvert_type==2 case
+        else
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
@@ -529,16 +542,20 @@
         ENDDO
 
-        DO ij =  1,iim
-          tppn(ij)  = aire(  ij    ) * ps (  ij    )
-          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
-        ENDDO
-          tpn  = SSUM(iim,tppn,1)/apoln
-          tps  = SSUM(iim,tpps,1)/apols
-
-        DO ij = 1, iip1
-          ps(  ij    ) = tpn
-          ps(ij+ip1jm) = tps
-        ENDDO
-
+        if (1 == 0) then
+!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
+!!!                     2) should probably not be here anyway
+!!! but are kept for those who would want to revert to previous behaviour
+           DO ij =  1,iim
+             tppn(ij)  = aire(  ij    ) * ps (  ij    )
+             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
+           ENDDO
+             tpn  = SSUM(iim,tppn,1)/apoln
+             tps  = SSUM(iim,tpps,1)/apols
+
+           DO ij = 1, iip1
+             ps(  ij    ) = tpn
+             ps(ij+ip1jm) = tps
+           ENDDO
+        endif ! of if (1 == 0)
 
       END IF ! of IF(apdiss)
@@ -622,5 +639,4 @@
              ! Ehouarn: output only during LF or Backward Matsuno
 	     if (leapf.or.(.not.leapf.and.(.not.forward))) then
-              nbetat = nbetatdem
               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
               unat=0.
@@ -652,9 +668,10 @@
 !              if (planet_type.eq."earth") then
 ! Write an Earth-format restart file
-                CALL dynredem1("restart.nc",0.0,
+                CALL dynredem1("restart.nc",start_time,
      &                         vcov,ucov,teta,q,masse,ps)
 !              endif ! of if (planet_type.eq."earth")
 
               CLOSE(99)
+              !!! Ehouarn: Why not stop here and now?
             ENDIF ! of IF (itau.EQ.itaufin)
 
@@ -740,5 +757,4 @@
               IF(MOD(itau,iecri         ).EQ.0) THEN
 c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
-                nbetat = nbetatdem
                 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
                 unat=0.
@@ -763,5 +779,5 @@
               IF(itau.EQ.itaufin) THEN
 !                if (planet_type.eq."earth") then
-                  CALL dynredem1("restart.nc",0.0,
+                  CALL dynredem1("restart.nc",start_time,
      &                           vcov,ucov,teta,q,masse,ps)
 !                endif ! of if (planet_type.eq."earth")
Index: LMDZ5/branches/testing/libf/dyn3d/temps.h
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/temps.h	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/temps.h	(revision 1665)
@@ -14,10 +14,12 @@
 
       COMMON/temps/itaufin, dt, day_ini, day_end, annee_ref, day_ref,   &
-     &             itau_dyn, itau_phy, jD_ref, jH_ref, calend
+     &             itau_dyn, itau_phy, jD_ref, jH_ref, calend,          &
+     &             start_time
+
 
       INTEGER   itaufin
       INTEGER itau_dyn, itau_phy
       INTEGER day_ini, day_end, annee_ref, day_ref
-      REAL      dt, jD_ref, jH_ref
+      REAL      dt, jD_ref, jH_ref, start_time
       CHARACTER (len=10) :: calend
 
Index: LMDZ5/branches/testing/libf/dyn3d/wrgrads.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3d/wrgrads.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3d/wrgrads.F	(revision 1665)
@@ -26,5 +26,5 @@
 c   local
 
-      integer im,jm,lm,i,j,l,lnblnk,iv,iii,iji,iif,ijf
+      integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
 
       logical writectl
@@ -59,5 +59,5 @@
             nvar(if)=ivar(if)
             var(ivar(if),if)=name
-            tvar(ivar(if),if)=titlevar(1:lnblnk(titlevar))
+            tvar(ivar(if),if)=trim(titlevar)
             nld(ivar(if),if)=nl
 c           print*,'initialisation ecriture de ',var(ivar(if),if)
@@ -101,8 +101,8 @@
       file=fichier(if)
 c   WARNING! on reecrase le fichier .ctl a chaque ecriture
-      open(unit(if),file=file(1:lnblnk(file))//'.ctl'
+      open(unit(if),file=trim(file)//'.ctl'
      &         ,form='formatted',status='unknown')
       write(unit(if),'(a5,1x,a40)')
-     &       'DSET ','^'//file(1:lnblnk(file))//'.dat'
+     &       'DSET ','^'//trim(file)//'.dat'
 
       write(unit(if),'(a12)') 'UNDEF 1.0E30'
