Index: trunk/LMDZ.COMMON/libf/dyn3d/calfis.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 776)
@@ -170,5 +170,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
@@ -461,4 +461,5 @@
       if (planet_type=="earth") then
 #ifdef CPP_EARTH
+! PVtheta calls tetalevel, which is in the (Earth) physics
 cIM calcul PV a teta=350, 380, 405K
       CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
@@ -482,6 +483,4 @@
 ! ne pose pas de probleme a priori.
 
-#ifdef CPP_PHYS
-
 !      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
       zdt_split=dtphys/nsplit_phys
@@ -490,4 +489,6 @@
       zdtfic(:,:)=0.
       zdqfic(:,:,:)=0.
+
+#ifdef CPP_PHYS
 
       do isplit=1,nsplit_phys
@@ -563,5 +564,9 @@
          zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
 
-      enddo
+      enddo ! of do isplit=1,nsplit_phys
+
+#endif
+! #endif of #ifdef CPP_PHYS
+
       zdufi(:,:)=zdufic(:,:)/nsplit_phys
       zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
@@ -569,6 +574,4 @@
       zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
 
-#endif
-! #endif of #ifdef CPP_PHYS
 
 500   CONTINUE
Index: trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 776)
@@ -28,5 +28,6 @@
   IMPLICIT NONE
 #ifndef CPP_EARTH
-  WRITE(*,*)'limit_netcdf: Earth-specific program, needs Earth physics'
+#include "iniprint.h"
+  WRITE(lunout,*)'limit_netcdf: Earth-specific program, needs Earth physics'
 #else
 !-------------------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3d/comvert.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/comvert.h	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/comvert.h	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: comvert.h 1520 2011-05-23 11:37:09Z emillour $
+! $Id: comvert.h 1625 2012-05-09 13:14:48Z lguez $
 !
 !-----------------------------------------------------------------------
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/disvert.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/disvert.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/disvert.F90	(revision 776)
@@ -1,7 +1,11 @@
-! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
+! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
 
 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
+
+  use new_unit_m, only: new_unit
+  use ioipsl, only: getin
+  use assert_m, only: assert
 
   IMPLICIT NONE
@@ -18,5 +22,6 @@
 
   real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: ap(llmp1) ! in Pa
+  real, intent(out):: bp(llmp1)
   real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
   real,intent(out) :: presnivs(llm)
@@ -26,20 +31,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
@@ -69,31 +82,12 @@
      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
+     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
+     bp(llmp1) = 0.
 
+     ap = pa * (sig - bp)
+  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 +96,56 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
-  ENDIF
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
+  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
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
+  case("read")
+     ! Read "ap" and "bp". First line is skipped (title line). "ap"
+     ! should be in Pa. First couple of values should correspond to
+     ! the surface, that is : "bp" should be in descending order.
+     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=*) ap(l), bp(l)
+     end do
+     close(unit)
+     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
+          bp(llm + 1) == 0., "disvert: bad ap or bp values")
+  case default
+     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
+  END select
 
   DO l=1, llm
@@ -111,19 +156,4 @@
      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 '
Index: trunk/LMDZ.COMMON/libf/dyn3d/dynetat0.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/dynetat0.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/dynetat0.F	(revision 776)
@@ -6,4 +6,5 @@
 
       USE infotrac
+      use netcdf, only: nf90_get_var
       IMPLICIT NONE
 
@@ -28,5 +29,5 @@
 #include "comconst.h"
 #include "comvert.h"
-#include "comgeom.h"
+#include "comgeom2.h"
 #include "ener.h"
 #include "netcdf.inc"
@@ -40,7 +41,7 @@
 
       CHARACTER*(*) fichnom
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
-      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
-      REAL ps(ip1jmp1),phis(ip1jmp1)
+      REAL vcov(iip1, jjm,llm),ucov(iip1, jjp1,llm),teta(iip1, jjp1,llm)
+      REAL q(iip1,jjp1,llm,nqtot),masse(iip1, jjp1,llm)
+      REAL ps(iip1, jjp1),phis(iip1, jjp1)
 
       REAL time
@@ -70,9 +71,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
-#endif
+      ierr = nf90_get_var(nid, nvarid, tab_cntrl)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echoue pour <controle>"
@@ -142,9 +139,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonu)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonu)
-#endif
+      ierr = nf90_get_var(nid, nvarid, rlonu)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>"
@@ -157,9 +150,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatu)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatu)
-#endif
+      ierr = nf90_get_var(nid, nvarid, rlatu)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>"
@@ -172,9 +161,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonv)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, rlonv)
-#endif
+      ierr = nf90_get_var(nid, nvarid, rlonv)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>"
@@ -187,9 +172,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatv)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, rlatv)
-#endif
+      ierr = nf90_get_var(nid, nvarid, rlatv)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour rlatv"
@@ -202,9 +183,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cu)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, cu)
-#endif
+      ierr = nf90_get_var(nid, nvarid, cu)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <cu>"
@@ -217,9 +194,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cv)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, cv)
-#endif
+      ierr = nf90_get_var(nid, nvarid, cv)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <cv>"
@@ -232,9 +205,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aire)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, aire)
-#endif
+      ierr = nf90_get_var(nid, nvarid, aire)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <aire>"
@@ -247,9 +216,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phis)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, phis)
-#endif
+      ierr = nf90_get_var(nid, nvarid, phis)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>"
@@ -262,9 +227,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, time)
-#endif
+      ierr = nf90_get_var(nid, nvarid, time)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee <temps>"
@@ -277,9 +238,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ucov)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, ucov)
-#endif
+      ierr = nf90_get_var(nid, nvarid, ucov)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <ucov>"
@@ -292,9 +249,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vcov)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, vcov)
-#endif
+      ierr = nf90_get_var(nid, nvarid, vcov)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <vcov>"
@@ -307,9 +260,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, teta)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, teta)
-#endif
+      ierr = nf90_get_var(nid, nvarid, teta)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <teta>"
@@ -325,11 +274,7 @@
      &                    "> est absent"
            write(lunout,*)"          Il est donc initialise a zero"
-           q(:,:,iq)=0.
+           q(:,:,:,iq)=0.
         ELSE
-#ifdef NC_DOUBLE
-          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q(1,1,iq))
-#else
-          ierr = NF_GET_VAR_REAL(nid, nvarid, q(1,1,iq))
-#endif
+           ierr = NF90_GET_VAR(nid, nvarid, q(:,:,:,iq))
           IF (ierr .NE. NF_NOERR) THEN
             write(lunout,*)"dynetat0: Lecture echouee pour "//tname(iq)
@@ -345,9 +290,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, masse)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, masse)
-#endif
+      ierr = nf90_get_var(nid, nvarid, masse)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <masse>"
@@ -360,9 +301,5 @@
          CALL abort
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ps)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, ps)
-#endif
+      ierr = nf90_get_var(nid, nvarid, ps)
       IF (ierr .NE. NF_NOERR) THEN
          write(lunout,*)"dynetat0: Lecture echouee pour <ps>"
Index: trunk/LMDZ.COMMON/libf/dyn3d/dynredem.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/dynredem.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/dynredem.F	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: dynredem.F 1403 2010-07-01 09:02:53Z fairhead $
+! $Id: dynredem.F 1635 2012-07-12 11:37:16Z lguez $
 !
 c
@@ -8,4 +8,5 @@
 #endif
       USE infotrac
+      use netcdf95, only: NF95_PUT_VAR
  
       IMPLICIT NONE
@@ -19,5 +20,5 @@
 #include "comconst.h"
 #include "comvert.h"
-#include "comgeom.h"
+#include "comgeom2.h"
 #include "temps.h"
 #include "ener.h"
@@ -31,5 +32,5 @@
 c   ----------
       INTEGER iday_end
-      REAL phis(ip1jmp1)
+      REAL phis(iip1, jjp1)
       CHARACTER*(*) fichnom
 
@@ -138,5 +139,5 @@
 c
       ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 27,
-     .                       "Fichier demarrage dynamique")
+     .                       "Fichier demmarage dynamique")
 c
 c Definir les dimensions du fichiers:
@@ -166,9 +167,5 @@
      .                       "Parametres de controle")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,tab_cntrl)
 c
       ierr = NF_REDEF (nid)
@@ -183,9 +180,5 @@
      .                       "Longitudes des points U")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,rlonu)
 c
       ierr = NF_REDEF (nid)
@@ -200,9 +193,5 @@
      .                       "Latitudes des points U")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu)
-#endif
+      call NF95_PUT_VAR (nid,nvarid,rlatu)
 c
       ierr = NF_REDEF (nid)
@@ -217,9 +206,5 @@
      .                       "Longitudes des points V")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,rlonv)
 c
       ierr = NF_REDEF (nid)
@@ -234,9 +219,5 @@
      .                       "Latitudes des points V")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,rlatv)
 c
       ierr = NF_REDEF (nid)
@@ -251,9 +232,5 @@
      .                       "Numero naturel des couches s")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsigs)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsigs)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,nivsigs)
 c
       ierr = NF_REDEF (nid)
@@ -268,9 +245,5 @@
      .                       "Numero naturel des couches sigma")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsig)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsig)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,nivsig)
 c
       ierr = NF_REDEF (nid)
@@ -285,9 +258,5 @@
      .                       "Coefficient A pour hybride")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,ap)
 c
       ierr = NF_REDEF (nid)
@@ -302,9 +271,5 @@
      .                       "Coefficient B pour hybride")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,bp)
 c
       ierr = NF_REDEF (nid)
@@ -317,9 +282,5 @@
 cIM 220306 END
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,presnivs)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,presnivs)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,presnivs)
 c
 c Coefficients de passage cov. <-> contra. <--> naturel
@@ -338,9 +299,5 @@
      .                       "Coefficient de passage pour U")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,cu)
 c
       ierr = NF_REDEF (nid)
@@ -357,9 +314,5 @@
      .                       "Coefficient de passage pour V")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,cv)
 c
 c Aire de chaque maille:
@@ -378,9 +331,5 @@
      .                       "Aires de chaque maille")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,aire)
 c
 c Geopentiel au sol:
@@ -399,9 +348,5 @@
      .                       "Geopotentiel au sol")
       ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,phis)
 c
 c Definir les variables pour pouvoir les enregistrer plus tard:
@@ -524,4 +469,6 @@
       USE infotrac
       USE control_mod
+      use netcdf, only: NF90_get_VAR
+      use netcdf95, only: NF95_PUT_VAR
  
       IMPLICIT NONE
@@ -538,9 +485,10 @@
 #include "iniprint.h"
 
+
       INTEGER l
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
-      REAL teta(ip1jmp1,llm)                   
-      REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
-      REAL q(ip1jmp1,llm,nqtot)
+      REAL vcov(iip1,jjm,llm),ucov(iip1, jjp1,llm) 
+      REAL teta(iip1, jjp1,llm)                   
+      REAL ps(iip1, jjp1),masse(iip1, jjp1,llm)                   
+      REAL q(iip1, jjp1, llm, nqtot)
       CHARACTER*(*) fichnom
      
@@ -576,9 +524,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR1_DOUBLE (nid,nvarid,nb,time)
-#else
-      ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,time,start=(/nb/))
       write(lunout,*) "dynredem1: Enregistrement pour ", nb, time
 
@@ -592,15 +536,7 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
-#endif
+      ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl)
        tab_cntrl(31) = REAL(itau_dyn + itaufin)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,tab_cntrl)
 
 c  Ecriture des champs
@@ -612,9 +548,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ucov)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,ucov)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,ucov)
 
       ierr = NF_INQ_VARID(nid, "vcov", nvarid)
@@ -624,9 +556,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,vcov)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,vcov)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,vcov)
 
       ierr = NF_INQ_VARID(nid, "teta", nvarid)
@@ -636,9 +564,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,teta)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,teta)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,teta)
 
       IF (type_trac == 'inca') THEN
@@ -662,9 +586,5 @@
                CALL abort_gcm(modname,abort_message,ierr)
             ENDIF
-#ifdef NC_DOUBLE
-            ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
-#else
-            ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
-#endif
+            call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
         ELSE ! type_trac = inca
 ! lecture de la valeur du traceur dans start_trac.nc
@@ -681,18 +601,10 @@
                    CALL abort_gcm(modname,abort_message,ierr)
                 ENDIF
-#ifdef NC_DOUBLE
-                ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
-#else
-                ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
-#endif
+                call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
                 
              ELSE
                 write(lunout,*) "dynredem1: ",trim(tname(iq)),
      &              " est present dans start_trac.nc"
-#ifdef NC_DOUBLE
-               ierr = NF_GET_VAR_DOUBLE(nid_trac, nvarid_trac, trac_tmp)
-#else
-               ierr = NF_GET_VAR_REAL(nid_trac, nvarid_trac, trac_tmp)
-#endif
+               ierr = NF90_GET_VAR(nid_trac, nvarid_trac, trac_tmp)
                 IF (ierr .NE. NF_NOERR) THEN
                    abort_message="dynredem1: Lecture echouee pour"//
@@ -708,9 +620,5 @@
                    CALL abort_gcm(modname,abort_message,ierr)
                 ENDIF
-#ifdef NC_DOUBLE
-                ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,trac_tmp)
-#else
-                ierr = NF_PUT_VAR_REAL (nid,nvarid,trac_tmp)
-#endif
+                call NF95_PUT_VAR(nid, nvarid, trac_tmp)
                
              ENDIF ! IF (ierr .NE. NF_NOERR)
@@ -725,9 +633,5 @@
                    CALL abort_gcm(modname,abort_message,ierr)
              ENDIF
-#ifdef NC_DOUBLE
-             ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
-#else
-             ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
-#endif
+             call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
           ENDIF ! (ierr_file .ne. 2)
        END IF   !type_trac
@@ -742,9 +646,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,masse)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,masse)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,masse)
 c
       ierr = NF_INQ_VARID(nid, "ps", nvarid)
@@ -754,9 +654,5 @@
          CALL abort_gcm(modname,abort_message,ierr)
       ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ps)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,ps)
-#endif
+      call NF95_PUT_VAR(nid,nvarid,ps)
 
       ierr = NF_CLOSE(nid)
Index: trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: etat0_netcdf.F90 1520 2011-05-23 11:37:09Z emillour $
+! $Id: etat0_netcdf.F90 1625 2012-05-09 13:14:48Z lguez $
 !
 !-------------------------------------------------------------------------------
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/exner_hyb.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/exner_hyb.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/exner_hyb.F	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/exner_milieu.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/exner_milieu.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/exner_milieu.F	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 776)
@@ -21,10 +21,8 @@
 ! 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:
+! Ehouarn: the following are needed with (parallel) physics:
 #ifdef CPP_PHYS
       USE dimphy
       USE comgeomphy
-#endif
-#ifdef CPP_EARTH
       USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
 #endif
@@ -177,12 +175,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
@@ -465,8 +459,7 @@
          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
      ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
-#endif ! CPP_PHYS
+#endif
          call_iniphys=.false.
       ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
-!#endif
 
 c  numero de stockage pour les fichiers de redemarrage:
Index: trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
+! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $
 !
 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F	(revision 617)
+++ 	(revision )
@@ -1,85 +1,0 @@
-!
-! $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
-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: trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F90	(revision 776)
@@ -0,0 +1,84 @@
+!
+! $Id: iniconst.F90 1625 2012-05-09 13:14:48Z lguez $
+!
+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: trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/inigrads.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/inigrads.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/inigrads.F	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3d/integrd.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 776)
@@ -1,8 +1,9 @@
 !
-! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $
+! $Id: integrd.F 1616 2012-02-17 11:59:00Z emillour $
 !
       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: trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 776)
@@ -124,5 +124,6 @@
 
       REAL  SSUM
-      REAL time_0 , finvmaold(ip1jmp1,llm)
+      REAL time_0
+!     REAL finvmaold(ip1jmp1,llm)
 
 cym      LOGICAL  lafin
@@ -243,7 +244,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
@@ -271,10 +272,10 @@
 c   ----------------------------------
 
-   1  CONTINUE
+   1  CONTINUE ! Matsuno Forward step begins here
 
       jD_cur = jD_ref + day_ini - day_ref +                             &
-     &          int (itau * dtvr / daysec) 
+     &          itau/day_step 
       jH_cur = jH_ref + start_time +                                    &
-     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     &          mod(itau,day_step)/float(day_step) 
       jD_cur = jD_cur + int(jH_cur)
       jH_cur = jH_cur - int(jH_cur)
@@ -307,9 +308,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-----------------------------------------------------------------------
@@ -357,4 +358,5 @@
       call tpot2t(ijp1llm,teta,temp,pk)
       tsurpk = cpp*temp/pk
+      ! compute geopotential phi()
       CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
 
@@ -372,4 +374,6 @@
 
 !      IF( forward. OR . leapf )  THEN
+! Ehouarn: NB: at this point p with ps are not synchronized
+!              (whereas mass and ps are...)
       IF((.not.forward).OR. leapf )  THEN
         ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
@@ -398,6 +402,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                                    )
 
        IF ((planet_type.eq."titan").and.(tidal)) then
@@ -431,14 +435,14 @@
 
          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
 
            jD_cur = jD_ref + day_ini - day_ref +                        &
-     &          int (itau * dtvr / daysec) 
+     &          itau/day_step 
            jH_cur = jH_ref + start_time +                               &
-     &          (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     &          mod(itau,day_step)/float(day_step) 
            jD_cur = jD_cur + int(jH_cur)
            jH_cur = jH_cur - int(jH_cur)
@@ -545,7 +549,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
@@ -613,16 +617,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)
@@ -749,4 +757,5 @@
 
               CLOSE(99)
+              !!! Ehouarn: Why not stop here and now?
             ENDIF ! of IF (itau.EQ.itaufin)
 
Index: trunk/LMDZ.COMMON/libf/dyn3d/wrgrads.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/wrgrads.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3d/wrgrads.F	(revision 776)
@@ -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'
