Index: LMDZ5/branches/testing/libf/dyn3dpar/bands.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/bands.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/bands.F90	(revision 1665)
@@ -93,6 +93,6 @@
    SUBROUTINE  Set_Bands 
      USE parallel
-#ifdef CPP_EARTH
-! Ehouarn: what follows is only related to // physics; for now only for Earth 
+#ifdef CPP_PHYS
+! Ehouarn: what follows is only related to // physics
      USE mod_phys_lmdz_para, ONLY : jj_para_begin,jj_para_end
 #endif
@@ -106,6 +106,5 @@
       enddo
           
-#ifdef CPP_EARTH
-! Ehouarn: what follows is only related to // physics; for now only for Earth          
+#ifdef CPP_PHYS
       do i=0,MPI_Size-1
         jj_Nb_physic(i)=jj_para_end(i)-jj_para_begin(i)+1
@@ -332,6 +331,6 @@
     subroutine AdjustBands_physic
       use times
-#ifdef CPP_EARTH
-! Ehouarn: what follows is only related to // physics; for now only for Earth 
+#ifdef CPP_PHYS
+! Ehouarn: what follows is only related to // physics
       USE mod_phys_lmdz_para, only : klon_mpi_para_nb
 #endif
@@ -359,6 +358,5 @@
       medium=medium/mpi_size      
       NbTot=0
-#ifdef CPP_EARTH
-! Ehouarn: what follows is only related to // physics; for now only for Earth 
+#ifdef CPP_PHYS
       do i=0,mpi_size-1
         Inc(i)=nint(klon_mpi_para_nb(i)*(medium-value(i))/value(i))
Index: LMDZ5/branches/testing/libf/dyn3dpar/calfis_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/calfis_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/calfis_p.F	(revision 1665)
@@ -27,9 +27,6 @@
      $                  pdqfi,
      $                  pdpsfi)
-#ifdef CPP_EARTH
-! Ehouarn: For now, calfis_p needs Earth physics
-c
-c    Auteur :  P. Le Van, F. Hourdin 
-c   .........
+#ifdef CPP_PHYS
+! If using physics
       USE dimphy
       USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
@@ -146,5 +143,6 @@
       REAL clesphy0( longcles )
 
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
+! Ehouarn: for now calfis_p needs some informations from physics to compile
 c    Local variables :
 c    -----------------
@@ -222,5 +220,5 @@
       PARAMETER(ntetaSTD=3)
       REAL rtetaSTD(ntetaSTD)
-      DATA rtetaSTD/350., 380., 405./
+      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
       REAL PVteta(klon,ntetaSTD)
       
@@ -489,6 +487,7 @@
 
 
-      IF (is_sequential) THEN
-c
+      IF (is_sequential.and.(planet_type=="earth")) THEN
+#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,
@@ -496,4 +495,5 @@
      $           ntetaSTD,rtetaSTD,PVteta)
 c
+#endif
       ENDIF
 
@@ -627,7 +627,4 @@
 c$OMP BARRIER
       
-      if (planet_type=="earth") then
-#ifdef CPP_EARTH
-
 !$OMP MASTER
 !      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
@@ -639,4 +636,6 @@
       zdqfic_omp(:,:,:)=0.
 
+      if (planet_type=="earth") then
+#ifdef CPP_PHYS
       do isplit=1,nsplit_phys
 
@@ -687,11 +686,12 @@
       enddo
 
+#endif
+! of #ifdef CPP_PHYS
+      endif !of if (planet_type=="earth")
+
       zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
       zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
       zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
       zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
-
-#endif
-      endif !of if (planet_type=="earth")
 c$OMP BARRIER
 
@@ -1110,10 +1110,10 @@
       firstcal = .FALSE.
 
-#else
-      write(lunout,*)
-     & "calfis_p: for now can only work with parallel physics"
-      stop
-#endif
-! of #ifdef CPP_EARTH
+#else 
+      write(lunout,*) 
+     & "calfis_p: for now can only work with parallel physics" 
+      stop 
+#endif 
+! of #ifdef CPP_PHYS
       RETURN
       END
Index: LMDZ5/branches/testing/libf/dyn3dpar/ce0l.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/ce0l.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/ce0l.F90	(revision 1665)
@@ -19,5 +19,5 @@
   USE dimphy
   USE comgeomphy
-  USE mod_phys_lmdz_para
+  USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
   USE mod_const_mpi
   USE infotrac
@@ -31,4 +31,5 @@
   IMPLICIT NONE
 #ifndef CPP_EARTH
+#include "iniprint.h"
   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
 #else
@@ -42,5 +43,10 @@
 #include "temps.h"
 #include "logic.h"
+#ifdef CPP_MPI
+      include 'mpif.h'
+#endif
+
   INTEGER, PARAMETER            :: longcles=20
+  INTEGER                       :: ierr
   REAL,    DIMENSION(longcles)  :: clesphy0
   REAL,    DIMENSION(iip1,jjp1) :: masque
@@ -50,5 +56,7 @@
   CALL conf_gcm( 99, .TRUE. , clesphy0 )
 
+#ifdef CPP_MPI
   CALL init_mpi
+#endif
 
   CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
@@ -115,7 +123,10 @@
   END IF
   
+#ifdef CPP_MPI
 !$OMP MASTER
-  CALL finalize_parallel
+  CALL MPI_FINALIZE(ierr)
+  IF (ierr /= 0) CALL abort_gcm('ce0l','Error in MPI_FINALIZE',1)
 !$OMP END MASTER
+#endif
 
 #endif
Index: LMDZ5/branches/testing/libf/dyn3dpar/comvert.h
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/comvert.h	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/conf_gcm.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/conf_gcm.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/conf_gcm.F	(revision 1665)
@@ -167,4 +167,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/dyn3dpar/control_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/control_mod.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/disvert.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/disvert.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/dynetat0.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/dynetat0.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/dynredem.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/dynredem.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/dynredem_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/dynredem_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/dynredem_p.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/dyn3dpar/etat0_netcdf.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/etat0_netcdf.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/exner_hyb.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/exner_hyb.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/exner_hyb_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/exner_hyb_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/exner_hyb_p.F	(revision 1665)
@@ -60,11 +60,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/dyn3dpar/exner_milieu.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/exner_milieu.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/exner_milieu_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/exner_milieu_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/exner_milieu_p.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.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/dyn3dpar/filtreg_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/filtreg_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/filtreg_p.F	(revision 1665)
@@ -208,24 +208,39 @@
                IF( ifiltre.EQ.-2 )   THEN
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matrinvn(1,1,j), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
                ELSE IF ( griscal )     THEN
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matriceun(1,1,j), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0,
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
                ELSE 
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matricevn(1,1,j), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
@@ -236,8 +251,14 @@
                IF( ifiltre.EQ.-2 )   THEN
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matrinvs(1,1,j-jfiltsu+1), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
+     &                            champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
@@ -245,8 +266,14 @@
                   
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matriceus(1,1,j-jfiltsu+1), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
+     &                            champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
@@ -254,8 +281,14 @@
                   
                   DO j = jdfil,jffil
+#ifdef BLAS
                      CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
      &                    matricevs(1,1,j-jfiltsv+1), iim, 
      &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
      &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
+#else
+                     champ_fft(:,j-jdfil+1,:)
+     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
+     &                            champ_loc(:iim,j,:))
+#endif
                   ENDDO
                   
Index: LMDZ5/branches/testing/libf/dyn3dpar/gcm.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/gcm.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/gcm.F	(revision 1665)
@@ -20,6 +20,5 @@
       USE control_mod
 
-! Ehouarn: for now these only apply to Earth:
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
       USE mod_grid_phy_lmdz
       USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
@@ -87,6 +86,4 @@
 
       REAL zdtvr
-c      INTEGER nbetatmoy, nbetatdem,nbetat
-      INTEGER nbetatmoy, nbetatdem
 
 c   variables dynamiques
@@ -189,13 +186,10 @@
       call ini_getparam("out.def")
       call Read_Distrib
-! 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,mpi_size,distrib_phys)
 #endif
-      endif ! of if (planet_type.eq."earth")
       CALL set_bands
-#ifdef CPP_EARTH
-! Ehouarn: For now only Earth physics is parallel
+#ifdef CPP_PHYS
       CALL Init_interface_dyn_phys
 #endif
@@ -209,12 +203,9 @@
 c$OMP END PARALLEL
 
-! Ehouarn : temporarily (?) keep this only for Earth
-      if (planet_type.eq."earth") then
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
 c$OMP PARALLEL
       call InitComgeomphy
 c$OMP END PARALLEL 
 #endif
-      endif ! of if (planet_type.eq."earth")
 
 c-----------------------------------------------------------------------
@@ -323,4 +314,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
@@ -390,7 +392,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
@@ -445,11 +444,9 @@
          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))
@@ -484,4 +481,12 @@
  301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
  302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
+#endif
+
+#ifdef CPP_PHYS
+! Create start file (startphy.nc) and boundary conditions (limit.nc) 
+! for the Earth verstion
+       if (iflag_phys>=100) then
+          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
+       endif
 #endif
 
Index: LMDZ5/branches/testing/libf/dyn3dpar/gr_dyn_fi_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/gr_dyn_fi_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/gr_dyn_fi_p.F	(revision 1665)
@@ -3,7 +3,6 @@
 !
       SUBROUTINE gr_dyn_fi_p(nfield,im,jm,ngrid,pdyn,pfi)
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
 ! Interface with parallel physics,
-! for now this routine only works with Earth physics
       USE mod_interface_dyn_phys
       USE dimphy
@@ -40,10 +39,6 @@
       ENDDO
 c$OMP END DO NOWAIT
-#else
-      write(lunout,*) "gr_fi_dyn_p : This routine should not be called",
-     &   "without parallelized physics"
-      stop
 #endif
-! of #ifdef CPP_EARTH
+! of #ifdef CPP_PHYS
       RETURN
       END
Index: LMDZ5/branches/testing/libf/dyn3dpar/gr_fi_dyn_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/gr_fi_dyn_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/gr_fi_dyn_p.F	(revision 1665)
@@ -3,7 +3,6 @@
 !
       SUBROUTINE gr_fi_dyn_p(nfield,ngrid,im,jm,pfi,pdyn)
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
 ! Interface with parallel physics,
-! for now this routine only works with Earth physics
       USE mod_interface_dyn_phys
       USE dimphy
@@ -52,10 +51,6 @@
       ENDDO
 c$OMP END DO NOWAIT
-#else
-      write(lunout,*) "gr_fi_dyn_p : This routine should not be called",
-     &   "without parallelized physics"
-      stop
 #endif
-! of #ifdef CPP_EARTH
+! of #ifdef CPP_PHYS
       RETURN
       END
Index: LMDZ5/branches/testing/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/guide_p_mod.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/guide_p_mod.F90	(revision 1665)
@@ -455,5 +455,5 @@
 !       Calcul niveaux pression milieu de couches 
 	CALL pression_p( ip1jmp1, ap, bp, ps, p )
-	if (disvert_type==1) then
+	if (pressure_exner) then
           CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
 	else
@@ -755,7 +755,7 @@
     ELSE
 	CALL pression_p( ip1jmp1, ap, bp, psi, p )
-	if (disvert_type==1) then
+	if (pressure_exner) 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
+        else
           CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
         endif
Index: LMDZ5/branches/testing/libf/dyn3dpar/iniacademic.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/iniacademic.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/iniconst.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/iniconst.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/iniconst.F90	(revision 1665)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/inidissip.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/inidissip.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/inigrads.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/inigrads.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/integrd_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/integrd_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/integrd_p.F	(revision 1665)
@@ -4,5 +4,5 @@
       SUBROUTINE integrd_p
      $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
-     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis,finvmaold)
+     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
       USE parallel
       USE control_mod, only : planet_type
@@ -33,20 +33,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 ps0(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) :: ps0(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:
@@ -54,5 +64,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,SAVE :: p(ip1jmp1,llmp1)
       REAL tpn,tps,tppn(iim),tpps(iim)
@@ -60,5 +71,5 @@
       REAL,SAVE :: deltap( ip1jmp1,llm )
 
-      INTEGER  l,ij,iq
+      INTEGER  l,ij,iq,i,j
 
       REAL SSUM
@@ -126,8 +137,12 @@
        
         IF( .NOT. checksum ) THEN
-         PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. '
-     &         , ps(stop_it)
-         print *, ' dans integrd'
-         stop 1
+         write(lunout,*) "integrd: negative surface pressure ",
+     &                                                ps(stop_it)
+         write(lunout,*) " at node ij =", stop_it
+         ! since ij=j+(i-1)*jjp1 , we have
+         j=modulo(stop_it,jjp1)
+         i=1+(stop_it-j)/jjp1
+         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
+     &                   " lat = ",rlatu(j)*180./pi, " deg"
         ENDIF
 
@@ -167,17 +182,19 @@
       CALL massdair_p (     p  , masse         )
 
-c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
-      ijb=ij_begin
-      ije=ij_end
-      
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
-      DO  l = 1,llm
-        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
-      ENDDO
-c$OMP END DO NOWAIT
-
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
+! Ehouarn : we don't use/need finvmaold and finvmasse,
+!           so might as well not compute them
+!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
+!      ijb=ij_begin
+!      ije=ij_end
+!      
+!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
+!      DO  l = 1,llm
+!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
+!      ENDDO
+!c$OMP END DO NOWAIT
+!
+!      jjb=jj_begin
+!      jje=jj_end
+!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
 c
 
@@ -330,11 +347,12 @@
       ENDIF
       
-c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
-
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      DO l = 1, llm      
-        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)        
-      ENDDO
-c$OMP END DO NOWAIT
+! Ehouarn: forget about finvmaold
+!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
+!
+!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+!      DO l = 1, llm      
+!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)        
+!      ENDDO
+!c$OMP END DO NOWAIT
 
       endif ! of if (planet_type.eq."earth")
Index: LMDZ5/branches/testing/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/leapfrog_p.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/leapfrog_p.F	(revision 1665)
@@ -75,5 +75,4 @@
 
       real zqmin,zqmax
-      INTEGER nbetatmoy, nbetatdem,nbetat
 
 c   variables dynamiques
@@ -125,5 +124,5 @@
       REAL  SSUM
       REAL time_0 
-      REAL,SAVE :: finvmaold(ip1jmp1,llm)
+!      REAL,SAVE :: finvmaold(ip1jmp1,llm)
 
 cym      LOGICAL  lafin
@@ -234,7 +233,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
@@ -245,9 +244,14 @@
 c et du parallelisme !!
 
-   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) 
+      if (jH_cur > 1.0 ) then
+        jD_cur = jD_cur +1.
+        jH_cur = jH_cur -1.
+      endif
 
 
@@ -280,7 +284,8 @@
          massem1= masse
          psm1= ps
-         
-         finvmaold = masse
-         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
+
+! Ehouarn: finvmaold is actually not used       
+!         finvmaold = masse
+!         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
 c$OMP END MASTER
 c$OMP BARRIER
@@ -300,5 +305,5 @@
            tetam1   (ijb:ije,l) = teta  (ijb:ije,l)
            massem1  (ijb:ije,l) = masse (ijb:ije,l)
-           finvmaold(ijb:ije,l)=masse(ijb:ije,l)
+!           finvmaold(ijb:ije,l)=masse(ijb:ije,l)
                  
            if (pole_sud) ije=ij_end-iip1
@@ -309,7 +314,7 @@
 c$OMP ENDDO  
 
-
-          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1, 
-     .                    llm, -2,2, .TRUE., 1 )
+! Ehouarn: finvmaold not used
+!          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1, 
+!     .                    llm, -2,2, .TRUE., 1 )
 
        endif ! of if (FirstCaldyn)
@@ -327,5 +332,5 @@
 cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
-   2  CONTINUE
+   2  CONTINUE ! Matsuno backward or leapfrog step begins here
 
 c$OMP MASTER
@@ -472,6 +477,6 @@
          call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
      &                                jj_Nb_caldyn,0,0,TestRequest)
-         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
-     &                                jj_Nb_caldyn,0,0,TestRequest)
+!         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
+!     &                                jj_Nb_caldyn,0,0,TestRequest)
  
         do j=1,nqtot
@@ -555,4 +560,5 @@
       call start_timer(timer_caldyn)
 
+      ! compute geopotential phi()
       CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 
@@ -629,6 +635,6 @@
 
        CALL integrd_p ( 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                                    )
 
 !       CALL FTRACE_REGION_END("integrd")
@@ -693,15 +699,19 @@
 
 c$OMP BARRIER
-         if (disvert_type==1) then
+         if (pressure_exner) 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
+         else
            CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
          endif
 c$OMP BARRIER
            jD_cur = jD_ref + day_ini - day_ref
-     $        + int (itau * dtvr / daysec) 
-           jH_cur = jH_ref +                                            &
-     &              (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     $        + itau/day_step
+           jH_cur = jH_ref + start_time +                                &
+     &              mod(itau,day_step)/float(day_step) 
 !         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
+           if (jH_cur > 1.0 ) then
+             jD_cur = jD_cur +1.
+             jH_cur = jH_cur -1.
+           endif
 
 c rajout debug
@@ -719,6 +729,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 
@@ -1036,7 +1048,7 @@
         CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
 c$OMP BARRIER
-        if (disvert_type==1) then
+        if (pressure_exner) 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
+        else
           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
@@ -1185,4 +1197,8 @@
 c$OMP END DO NOWAIT
 
+         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
 c$OMP MASTER               
           DO ij =  1,iim
@@ -1195,5 +1211,6 @@
           ENDDO
 c$OMP END MASTER
-        endif
+         endif ! of if (1 == 0)
+        endif ! of of (pole_nord)
         
         if (pole_sud) then
@@ -1211,4 +1228,8 @@
 c$OMP END DO NOWAIT
 
+         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
 c$OMP MASTER               
           DO ij =  1,iim
@@ -1221,5 +1242,6 @@
           ENDDO
 c$OMP END MASTER
-        endif
+         endif ! of if (1 == 0)
+        endif ! of if (pole_sud)
 
 
@@ -1431,5 +1453,4 @@
 c$OMP BARRIER
 c$OMP MASTER
-              nbetat = nbetatdem
               CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
        
@@ -1634,5 +1655,4 @@
 c$OMP BARRIER
 c$OMP MASTER
-                nbetat = nbetatdem
                 CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
 
Index: LMDZ5/branches/testing/libf/dyn3dpar/mod_interface_dyn_phys.F90
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/mod_interface_dyn_phys.F90	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/mod_interface_dyn_phys.F90	(revision 1665)
@@ -7,7 +7,6 @@
   
   
-#ifdef CPP_EARTH
+#ifdef CPP_PHYS
 ! Interface with parallel physics,
-! for now this routine only works with Earth physics
 CONTAINS
   
@@ -56,4 +55,4 @@
   END SUBROUTINE Init_interface_dyn_phys 
 #endif
-! of #ifdef CPP_EARTH
+! of #ifdef CPP_PHYS
 END MODULE mod_interface_dyn_phys
Index: LMDZ5/branches/testing/libf/dyn3dpar/temps.h
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/temps.h	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/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/dyn3dpar/wrgrads.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dpar/wrgrads.F	(revision 1664)
+++ LMDZ5/branches/testing/libf/dyn3dpar/wrgrads.F	(revision 1665)
@@ -22,5 +22,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
@@ -55,5 +55,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
             print*,'initialisation ecriture de ',var(ivar(if),if)
@@ -96,8 +96,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'