Index: trunk/LMDZ.COMMON/libf/bibio/initdynav.F
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/initdynav.F	(revision 617)
+++ 	(revision )
@@ -1,200 +1,0 @@
-!
-! $Id: initdynav.F 1403 2010-07-01 09:02:53Z fairhead $
-!
-      subroutine initdynav(day0,anne0,tstep,t_ops,t_wrt)
-
-#ifdef CPP_IOIPSL
-       USE IOIPSL
-#endif
-       USE infotrac, ONLY : nqtot, ttext
-      use com_io_dyn_mod, only : histaveid,histvaveid,histuaveid,       &
-     &        dynhistave_file,dynhistvave_file,dynhistuave_file
-      implicit none
-
-C
-C   Routine d'initialisation des ecritures des fichiers histoires LMDZ
-C   au format IOIPSL. Initialisation du fichier histoire moyenne.
-C
-C   Appels succesifs des routines: histbeg
-C                                  histhori
-C                                  histver
-C                                  histdef
-C                                  histend
-C
-C   Entree:
-C
-C      infile: nom du fichier histoire a creer
-C      day0,anne0: date de reference
-C      tstep : frequence d'ecriture
-C      t_ops: frequence de l'operation pour IOIPSL
-C      t_wrt: frequence d'ecriture sur le fichier
-C
-C
-C   L. Fairhead, LMD, 03/99
-C
-C =====================================================================
-C
-C   Declarations
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "temps.h"
-#include "ener.h"
-#include "logic.h"
-#include "description.h"
-#include "serre.h"
-#include "iniprint.h"
-
-C   Arguments
-C
-      integer day0, anne0
-      real tstep, t_ops, t_wrt
-
-#ifdef CPP_IOIPSL
-! This routine needs IOIPSL to work
-C   Variables locales
-C
-      integer tau0
-      real zjulian
-      integer iq
-      real rlong(iip1,jjp1), rlat(iip1,jjp1)
-      integer uhoriid, vhoriid, thoriid, zvertiid
-      integer ii,jj
-      integer zan, dayref
-C
-C  Initialisations
-C
-      pi = 4. * atan (1.)
-C
-C  Appel a histbeg: creation du fichier netcdf et initialisations diverses
-C         
-
-      zan = anne0
-      dayref = day0
-      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
-      tau0 = itau_dyn
-      
-      do jj = 1, jjp1
-        do ii = 1, iip1
-          rlong(ii,jj) = rlonv(ii) * 180. / pi
-          rlat(ii,jj)  = rlatu(jj) * 180. / pi
-        enddo
-      enddo
-       
-! Creation de 3 fichiers pour les differentes grilles horizontales
-! Restriction de IOIPSL: seulement 2 coordonnees dans le meme fichier
-! Grille Scalaire       
-      call histbeg(dynhistave_file, iip1, rlong(:,1), jjp1, rlat(1,:),
-     .             1, iip1, 1, jjp1,
-     .             tau0, zjulian, tstep, thoriid,histaveid)
-
-C  Creation du fichier histoire pour les grilles en V et U (oblige pour l'instant,
-C  IOIPSL ne permet pas de grilles avec des nombres de point differents dans
-C  un meme fichier)
-! Grille V
-      do jj = 1, jjm
-        do ii = 1, iip1
-          rlong(ii,jj) = rlonv(ii) * 180. / pi
-          rlat(ii,jj) = rlatv(jj) * 180. / pi
-        enddo
-      enddo
-
-      call histbeg(dynhistvave_file, iip1, rlong(:,1), jjm, rlat(1,:),
-     .             1, iip1, 1, jjm,
-     .             tau0, zjulian, tstep, vhoriid,histvaveid)
-! Grille U
-      do jj = 1, jjp1
-        do ii = 1, iip1
-          rlong(ii,jj) = rlonu(ii) * 180. / pi
-          rlat(ii,jj) = rlatu(jj) * 180. / pi
-        enddo
-      enddo
-
-      call histbeg(dynhistuave_file, iip1, rlong(:,1),jjp1, rlat(1,:),
-     .             1, iip1, 1, jjp1,
-     .             tau0, zjulian, tstep, uhoriid,histuaveid)
-C
-C  Appel a histvert pour la grille verticale
-C
-      call histvert(histaveid,'presnivs','Niveaux Pression
-     &     approximatifs','mb',llm, presnivs/100., zvertiid,'down')
-      call histvert(histuaveid,'presnivs','Niveaux Pression
-     &     approximatifs','mb',llm, presnivs/100., zvertiid,'down')
-      call histvert(histvaveid,'presnivs','Niveaux Pression
-     &     approximatifs','mb',llm, presnivs/100., zvertiid,'down')
-C
-C  Appels a histdef pour la definition des variables a sauvegarder
-C
-C  Vents U
-C
-!      write(6,*)'inithistave',tstep
-      call histdef(histuaveid, 'u', 'vent u moyen ',
-     .             'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-
-C  Vents V
-C
-      call histdef(histvaveid, 'v', 'vent v moyen',
-     .             'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-
-C
-C  Temperature
-C
-      call histdef(histaveid, 'temp', 'temperature moyenne', 'K',
-     .             iip1, jjp1, thoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-C
-C  Temperature potentielle
-C
-      call histdef(histaveid, 'theta', 'temperature potentielle', 'K',
-     .             iip1, jjp1, thoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-C
-C  Geopotentiel
-C
-      call histdef(histaveid, 'phi', 'geopotentiel moyen', '-',
-     .             iip1, jjp1, thoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-C
-C  Traceurs
-C
-!        DO iq=1,nqtot
-!          call histdef(histaveid, ttext(iq), ttext(iq), '-',
-!     .             iip1, jjp1, thoriid, llm, 1, llm, zvertiid,
-!     .             32, 'ave(X)', t_ops, t_wrt)
-!        enddo
-C
-C  Masse
-C
-      call histdef(histaveid, 'masse', 'masse', 'kg',
-     .             iip1, jjp1, thoriid, llm, 1, llm, zvertiid,
-     .             32, 'ave(X)', t_ops, t_wrt)
-C
-C  Pression au sol
-C
-      call histdef(histaveid, 'ps', 'pression naturelle au sol', 'Pa',
-     .             iip1, jjp1, thoriid, 1, 1, 1, -99,
-     .             32, 'ave(X)', t_ops, t_wrt)
-C
-C  Geopotentiel au sol
-C
-!      call histdef(histaveid, 'phis', 'geopotentiel au sol', '-',
-!     .             iip1, jjp1, thoriid, 1, 1, 1, -99,
-!     .             32, 'ave(X)', t_ops, t_wrt)
-!C
-C  Fin
-C
-      call histend(histaveid)
-      call histend(histuaveid)
-      call histend(histvaveid)
-#else
-! tell the user this routine should be run with ioipsl
-      write(lunout,*)"initdynav: Warning this routine should not be",
-     &               " used without ioipsl"
-#endif
-! of #ifdef CPP_IOIPSL
-      return
-      end
Index: trunk/LMDZ.COMMON/libf/bibio/initdynav.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/initdynav.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/bibio/initdynav.F90	(revision 776)
@@ -0,0 +1,194 @@
+! $Id: initdynav.F90 1611 2012-01-25 14:31:54Z lguez $
+
+subroutine initdynav(day0,anne0,tstep,t_ops,t_wrt)
+
+#ifdef CPP_IOIPSL
+  USE IOIPSL
+#endif
+  USE infotrac, ONLY : nqtot, ttext
+  use com_io_dyn_mod, only : histaveid,histvaveid,histuaveid, &
+       dynhistave_file,dynhistvave_file,dynhistuave_file
+  implicit none
+
+
+  !   Routine d'initialisation des ecritures des fichiers histoires LMDZ
+  !   au format IOIPSL. Initialisation du fichier histoire moyenne.
+
+  !   Appels succesifs des routines: histbeg
+  !                                  histhori
+  !                                  histver
+  !                                  histdef
+  !                                  histend
+
+  !   Entree:
+
+  !      infile: nom du fichier histoire a creer
+  !      day0,anne0: date de reference
+  !      tstep : frequence d'ecriture
+  !      t_ops: frequence de l'operation pour IOIPSL
+  !      t_wrt: frequence d'ecriture sur le fichier
+
+
+  !   L. Fairhead, LMD, 03/99
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "comgeom.h"
+  include "temps.h"
+  include "ener.h"
+  include "logic.h"
+  include "description.h"
+  include "serre.h"
+  include "iniprint.h"
+
+  !   Arguments
+
+  integer day0, anne0
+  real tstep, t_ops, t_wrt
+
+#ifdef CPP_IOIPSL
+  ! This routine needs IOIPSL to work
+  !   Variables locales
+
+  integer tau0
+  real zjulian
+  integer iq
+  real rlong(iip1,jjp1), rlat(iip1,jjp1)
+  integer uhoriid, vhoriid, thoriid, zvertiid
+  integer ii,jj
+  integer zan, dayref
+
+  !--------------------------------------------------------------------
+
+  !  Initialisations
+
+  pi = 4. * atan (1.)
+
+  !  Appel a histbeg: creation du fichier netcdf et initialisations diverses
+
+
+  zan = anne0
+  dayref = day0
+  CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
+  tau0 = itau_dyn
+
+  do jj = 1, jjp1
+     do ii = 1, iip1
+        rlong(ii,jj) = rlonv(ii) * 180. / pi
+        rlat(ii,jj)  = rlatu(jj) * 180. / pi
+     enddo
+  enddo
+
+  ! Creation de 3 fichiers pour les differentes grilles horizontales
+  ! Restriction de IOIPSL: seulement 2 coordonnees dans le meme fichier
+  ! Grille Scalaire       
+  call histbeg(dynhistave_file, iip1, rlong(:,1), jjp1, rlat(1,:), &
+       1, iip1, 1, jjp1, &
+       tau0, zjulian, tstep, thoriid,histaveid)
+
+  ! Creation du fichier histoire pour les grilles en V et U (oblige
+  ! pour l'instant, IOIPSL ne permet pas de grilles avec des nombres
+  ! de point differents dans  un meme fichier)
+  ! Grille V
+  do jj = 1, jjm
+     do ii = 1, iip1
+        rlong(ii,jj) = rlonv(ii) * 180. / pi
+        rlat(ii,jj) = rlatv(jj) * 180. / pi
+     enddo
+  enddo
+
+  call histbeg(dynhistvave_file, iip1, rlong(:,1), jjm, rlat(1,:), &
+       1, iip1, 1, jjm, &
+       tau0, zjulian, tstep, vhoriid,histvaveid)
+  ! Grille U
+  do jj = 1, jjp1
+     do ii = 1, iip1
+        rlong(ii,jj) = rlonu(ii) * 180. / pi
+        rlat(ii,jj) = rlatu(jj) * 180. / pi
+     enddo
+  enddo
+
+  call histbeg(dynhistuave_file, iip1, rlong(:,1),jjp1, rlat(1,:), &
+       1, iip1, 1, jjp1, &
+       tau0, zjulian, tstep, uhoriid,histuaveid)
+
+  !  Appel a histvert pour la grille verticale
+
+  call histvert(histaveid,'presnivs','Niveaux Pression approximatifs','mb', &
+       llm, presnivs/100., zvertiid,'down')
+  call histvert(histuaveid,'presnivs','Niveaux Pression approximatifs','mb', &
+       llm, presnivs/100., zvertiid,'down')
+  call histvert(histvaveid,'presnivs','Niveaux Pression approximatifs','mb', &
+       llm, presnivs/100., zvertiid,'down')
+
+  !  Appels a histdef pour la definition des variables a sauvegarder
+
+  !  Vents U
+
+  call histdef(histuaveid, 'u', 'vent u moyen ', &
+       'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Vents V
+
+  call histdef(histvaveid, 'v', 'vent v moyen', &
+       'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+
+  !  Temperature
+
+  call histdef(histaveid, 'temp', 'temperature moyenne', 'K', &
+       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Temperature potentielle
+
+  call histdef(histaveid, 'theta', 'temperature potentielle', 'K', &
+       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Geopotentiel
+
+  call histdef(histaveid, 'phi', 'geopotentiel moyen', '-', &
+       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Traceurs
+
+  !        DO iq=1,nqtot
+  !          call histdef(histaveid, ttext(iq), ttext(iq), '-', &
+  !                  iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+  !                  32, 'ave(X)', t_ops, t_wrt)
+  !        enddo
+
+  !  Masse
+
+  call histdef(histaveid, 'masse', 'masse', 'kg', &
+       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Pression au sol
+
+  call histdef(histaveid, 'ps', 'pression naturelle au sol', 'Pa', &
+       iip1, jjp1, thoriid, 1, 1, 1, -99, &
+       32, 'ave(X)', t_ops, t_wrt)
+
+  !  Geopotentiel au sol
+
+  !      call histdef(histaveid, 'phis', 'geopotentiel au sol', '-', &
+  !                  iip1, jjp1, thoriid, 1, 1, 1, -99, &
+  !                  32, 'ave(X)', t_ops, t_wrt)
+
+  call histend(histaveid)
+  call histend(histuaveid)
+  call histend(histvaveid)
+#else
+  write(lunout,*)"initdynav: Warning this routine should not be", &
+       " used without ioipsl"
+#endif
+  ! of #ifdef CPP_IOIPSL
+
+end subroutine initdynav
Index: trunk/LMDZ.COMMON/libf/bibio/netcdf95.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/netcdf95.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/bibio/netcdf95.F90	(revision 776)
@@ -3,38 +3,10 @@
 
   ! Author: Lionel GUEZ
-
-  ! Three criticisms may be made about the Fortran 90 NetCDF interface:
-
-  ! -- NetCDF procedures are usually functions with side effects.
-  ! First, they have "intent(out)" arguments.
-  ! Furthermore, there is obviously data transfer inside the procedures.
-  ! Any data transfer inside a function is considered as a side effect.
-
-  ! -- The caller of a NetCDF procedure usually has to handle the error
-  ! status. NetCDF procedures would be much friendlier if they behaved
-  ! like the Fortran input/output statements. That is, the error status
-  ! should be an optional output argument.
-  ! If the caller does not request the error status and there is an
-  ! error then the NetCDF procedure should produce an error message
-  ! and stop the program.
-
-  ! -- Some procedures use array arguments with assumed size.
-  ! It would be better to use the pointer attribute.
-
-  ! This module produces a NetCDF interface that answers those three
-  ! criticisms for some (not all) procedures.
-
-  ! "nf95_get_att" is more secure than "nf90_get_att" because it
-  ! checks that the "values" argument is long enough and removes the
-  ! null terminator, if any.
-
-  ! This module replaces some of the official NetCDF procedures.
-  ! This module also provides the procedures "handle_err" and "nf95_gw_var".
-
-  ! This module provides only a partial replacement for some generic
-  ! procedures such as "nf90_def_var".
+  ! See:
+  ! http://www.lmd.jussieu.fr/~lglmd/NetCDF95
 
   use nf95_def_var_m
   use nf95_put_var_m
+  use nf95_get_var_m
   use nf95_gw_var_m
   use nf95_put_att_m
Index: trunk/LMDZ.COMMON/libf/bibio/new_unit.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/new_unit.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/bibio/new_unit.F90	(revision 776)
@@ -0,0 +1,25 @@
+module new_unit_m
+
+  implicit none
+
+contains
+
+  subroutine new_unit(unit)
+
+    integer, intent(out):: unit
+
+    ! Variables local to the procedure:
+    logical opened, exist
+
+    !------------------------------------------------------
+
+    unit = 0
+    do
+       inquire(unit=unit, opened=opened, exist=exist)
+       if (exist .and. .not. opened) exit
+       unit = unit + 1
+    end do
+
+  end subroutine new_unit
+
+end module new_unit_m
Index: trunk/LMDZ.COMMON/libf/bibio/nf95_get_att_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/nf95_get_att_m.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/bibio/nf95_get_att_m.F90	(revision 776)
@@ -1,9 +1,16 @@
 ! $Id$
 module nf95_get_att_m
+
+  use handle_err_m, only: handle_err
+  use netcdf, only: nf90_get_att, nf90_noerr
+  use simple, only: nf95_inquire_attribute
 
   implicit none
 
   interface nf95_get_att
-     module procedure nf95_get_att_text
+     module procedure nf95_get_att_text, nf95_get_att_one_FourByteInt
+
+     ! The difference between the specific procedures is the type of
+     ! argument "values".
   end interface
 
@@ -15,7 +22,4 @@
   subroutine nf95_get_att_text(ncid, varid, name, values, ncerr)
 
-    use netcdf, only: nf90_get_att, nf90_inquire_attribute, nf90_noerr
-    use handle_err_m, only: handle_err
-
     integer,                          intent( in) :: ncid, varid
     character(len = *),               intent( in) :: name
@@ -23,5 +27,5 @@
     integer, intent(out), optional:: ncerr
 
-    ! Variable local to the procedure:
+    ! Variables local to the procedure:
     integer ncerr_not_opt
     integer att_len
@@ -30,15 +34,16 @@
 
     ! Check that the length of "values" is large enough:
-    ncerr_not_opt = nf90_inquire_attribute(ncid, varid, name, len=att_len)
-    call handle_err("nf95_get_att_text nf90_inquire_attribute " &
-         // trim(name), ncerr_not_opt, ncid, varid)
-    if (len(values) < att_len) then
-       print *, "nf95_get_att_text"
-       print *, "varid = ", varid
-       print *, "attribute name: ", name
-       print *, 'length of "values" is not large enough'
-       print *, "len(values) = ", len(values)
-       print *, "number of characters in attribute: ", att_len
-       stop 1
+    call nf95_inquire_attribute(ncid, varid, name, nclen=att_len, &
+         ncerr=ncerr_not_opt)
+    if (ncerr_not_opt == nf90_noerr) then
+       if (len(values) < att_len) then
+          print *, "nf95_get_att_text"
+          print *, "varid = ", varid
+          print *, "attribute name: ", name
+          print *, 'length of "values" is not large enough'
+          print *, "len(values) = ", len(values)
+          print *, "number of characters in attribute: ", att_len
+          stop 1
+       end if
     end if
 
@@ -48,5 +53,6 @@
        ncerr = ncerr_not_opt
     else
-       call handle_err("nf95_get_att_text", ncerr_not_opt, ncid, varid)
+       call handle_err("nf95_get_att_text " // trim(name), ncerr_not_opt, &
+            ncid, varid)
     end if
 
@@ -58,3 +64,42 @@
   end subroutine nf95_get_att_text
 
+  !***********************
+
+  subroutine nf95_get_att_one_FourByteInt(ncid, varid, name, values, ncerr)
+
+    integer,                                    intent( in) :: ncid, varid
+    character(len = *),                         intent( in) :: name
+    integer ,               intent(out) :: values
+    integer, intent(out), optional:: ncerr
+
+    ! Variables local to the procedure:
+    integer ncerr_not_opt
+    integer att_len
+
+    !-------------------
+
+    ! Check that the attribute contains a single value:
+    call nf95_inquire_attribute(ncid, varid, name, nclen=att_len, &
+         ncerr=ncerr_not_opt)
+    if (ncerr_not_opt == nf90_noerr) then
+       if (att_len /= 1) then
+          print *, "nf95_get_att_one_FourByteInt"
+          print *, "varid = ", varid
+          print *, "attribute name: ", name
+          print *, 'the attribute does not contain a single value'
+          print *, "number of values in attribute: ", att_len
+          stop 1
+       end if
+    end if
+
+    ncerr_not_opt = nf90_get_att(ncid, varid, name, values)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_att_one_FourByteInt " // trim(name), &
+            ncerr_not_opt, ncid, varid)
+    end if
+
+  end subroutine nf95_get_att_one_FourByteInt
+
 end module nf95_get_att_m
Index: trunk/LMDZ.COMMON/libf/bibio/nf95_get_var_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/nf95_get_var_m.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/bibio/nf95_get_var_m.F90	(revision 776)
@@ -0,0 +1,391 @@
+module nf95_get_var_m
+
+  use netcdf, only: nf90_get_var
+  use handle_err_m, only: handle_err
+
+  implicit none
+
+  interface nf95_get_var
+     module procedure nf95_get_var_FourByteReal, nf95_get_var_FourByteInt, &
+          nf95_get_var_1D_FourByteReal, nf95_get_var_1D_FourByteInt, &
+          nf95_get_var_2D_FourByteReal, &
+          nf95_get_var_3D_FourByteInt, &
+          nf95_get_var_3D_FourByteReal, &
+          nf95_get_var_4D_FourByteReal, &
+          nf95_get_var_5D_FourByteReal
+  end interface
+
+  private
+  public nf95_get_var
+
+contains
+
+  subroutine nf95_get_var_FourByteReal(ncid, varid, values, start, ncerr)
+
+    integer, intent(in) :: ncid, varid
+    real, intent(out) :: values
+    integer, dimension(:), optional, intent(in) :: start
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_FourByteInt(ncid, varid, values, start, ncerr)
+
+    integer, intent(in) :: ncid, varid
+    integer, intent(out) :: values
+    integer, dimension(:), optional, intent(in) :: start
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_FourByteInt
+
+  !***********************
+
+  subroutine nf95_get_var_1D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer,                         intent(in) :: ncid, varid
+    real, intent(out) :: values(:)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_1D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_1D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_1D_FourByteInt(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer,                         intent(in) :: ncid, varid
+    integer, intent(out) :: values(:)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_1D_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_1D_FourByteInt
+
+  !***********************
+
+  subroutine nf95_get_var_1D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: eightByteReal
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal),     intent(out) :: values(:)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_1D_eightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_1D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_2D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer,                         intent(in) :: ncid, varid
+    real , intent(out) :: values(:, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_2D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_2D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_2D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: EightByteReal
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal), intent(out) :: values(:, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_2D_EightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_2D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_3D_FourByteInt(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer, intent(in):: ncid, varid
+    integer, intent(out):: values(:, :, :)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_3D_FourByteInt", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_3D_FourByteInt
+
+  !***********************
+
+  subroutine nf95_get_var_3D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer,                         intent(in) :: ncid, varid
+    real , intent(out) :: values(:, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_3D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_3D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_3D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: eightByteReal
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal),     intent(out) :: values(:, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_3D_eightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_3D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_4D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer,                         intent(in) :: ncid, varid
+    real , intent(out) :: values(:, :, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_4D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_4D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_4D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: EightByteReal
+
+    integer, intent(in):: ncid, varid
+    real(kind = EightByteReal), intent(out):: values(:, :, :, :)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_4D_EightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_4D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_5D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    integer, intent(in):: ncid, varid
+    real, intent(out):: values(:, :, :, :, :)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_5D_FourByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_5D_FourByteReal
+
+  !***********************
+
+  subroutine nf95_get_var_5D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: EightByteReal
+
+    integer, intent(in):: ncid, varid
+    real(kind = EightByteReal), intent(out):: values(:, :, :, :, :)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_get_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_get_var_5D_EightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_get_var_5D_EightByteReal
+
+end module nf95_get_var_m
Index: trunk/LMDZ.COMMON/libf/bibio/nf95_gw_var_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/nf95_gw_var_m.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/bibio/nf95_gw_var_m.F90	(revision 776)
@@ -1,4 +1,7 @@
 ! $Id$
 module nf95_gw_var_m
+
+  use nf95_get_var_m, only: NF95_GET_VAR
+  use simple, only: nf95_inquire_variable, nf95_inquire_dimension
 
   implicit none
@@ -8,14 +11,14 @@
      ! These procedures read a whole NetCDF variable (coordinate or
      ! primary) into an array.
-     ! The difference between the procedures is the rank of the array
-     ! and the type of Fortran values.
+     ! The difference between the procedures is the rank and type of
+     ! argument "values".
      ! The procedures do not check the type of the NetCDF variable.
 
-!!$     module procedure nf95_gw_var_real_1d, nf95_gw_var_real_2d, &
-!!$          nf95_gw_var_real_3d, nf95_gw_var_real_4d, nf95_gw_var_dble_1d, &
-!!$          nf95_gw_var_dble_3d, nf95_gw_var_int_1d, nf95_gw_var_int_3d
+     ! Not including double precision procedures in the generic
+     ! interface because we use a compilation option that changes default
+     ! real precision.
      module procedure nf95_gw_var_real_1d, nf95_gw_var_real_2d, &
-          nf95_gw_var_real_3d, nf95_gw_var_real_4d, nf95_gw_var_int_1d, &
-          nf95_gw_var_int_3d
+          nf95_gw_var_real_3d, nf95_gw_var_real_4d, nf95_gw_var_real_5d, &
+          nf95_gw_var_int_1d, nf95_gw_var_int_3d
   end interface
 
@@ -29,8 +32,4 @@
     ! Real type, the array has rank 1.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -38,6 +37,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len
-    integer, pointer :: dimids(:)
+    integer nclen
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -46,16 +45,15 @@
 
     if (size(dimids) /= 1) then
-       print *, "nf95_gw_var_real_1d: NetCDF variable is not of rank 1"
-       stop 1
-    end if
-
-    call nf95_inquire_dimension(ncid, dimids(1), len=len)
-    deallocate(dimids) ! pointer
-
-    allocate(values(len))
-    if (len /= 0) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+       print *, "nf95_gw_var_real_1d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 1"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen)
+    deallocate(dimids) ! pointer
+
+    allocate(values(nclen))
+    if (nclen /= 0) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_real_1d
@@ -67,8 +65,4 @@
     ! Real type, the array has rank 2.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -76,6 +70,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len1, len2
-    integer, pointer :: dimids(:)
+    integer nclen1, nclen2
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -84,17 +78,16 @@
 
     if (size(dimids) /= 2) then
-       print *, "nf95_gw_var_real_2d: NetCDF variable is not of rank 2"
-       stop 1
-    end if
-
-    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
-    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
-    deallocate(dimids) ! pointer
-
-    allocate(values(len1, len2))
-    if (len1 /= 0 .and. len2 /= 0) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+       print *, "nf95_gw_var_real_2d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 2"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen1)
+    call nf95_inquire_dimension(ncid, dimids(2), nclen=nclen2)
+    deallocate(dimids) ! pointer
+
+    allocate(values(nclen1, nclen2))
+    if (nclen1 /= 0 .and. nclen2 /= 0) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_real_2d
@@ -106,8 +99,4 @@
     ! Real type, the array has rank 3.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -115,6 +104,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len1, len2, len3
-    integer, pointer :: dimids(:)
+    integer nclen1, nclen2, nclen3
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -123,18 +112,17 @@
 
     if (size(dimids) /= 3) then
-       print *, "nf95_gw_var_real_3d: NetCDF variable is not of rank 3"
-       stop 1
-    end if
-
-    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
-    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
-    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
-    deallocate(dimids) ! pointer
-
-    allocate(values(len1, len2, len3))
-    if (len1 * len2 * len3 /= 0) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+       print *, "nf95_gw_var_real_3d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 3"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen1)
+    call nf95_inquire_dimension(ncid, dimids(2), nclen=nclen2)
+    call nf95_inquire_dimension(ncid, dimids(3), nclen=nclen3)
+    deallocate(dimids) ! pointer
+
+    allocate(values(nclen1, nclen2, nclen3))
+    if (nclen1 * nclen2 * nclen3 /= 0) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_real_3d
@@ -146,8 +134,4 @@
     ! Real type, the array has rank 4.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -155,6 +139,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len_dim(4), i
-    integer, pointer :: dimids(:)
+    integer len_dim(4), i
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -163,18 +147,17 @@
 
     if (size(dimids) /= 4) then
-       print *, "nf95_gw_var_real_4d: NetCDF variable is not of rank 4"
+       print *, "nf95_gw_var_real_4d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 4"
        stop 1
     end if
 
     do i = 1, 4
-       call nf95_inquire_dimension(ncid, dimids(i), len=len_dim(i))
+       call nf95_inquire_dimension(ncid, dimids(i), nclen=len_dim(i))
     end do
     deallocate(dimids) ! pointer
 
     allocate(values(len_dim(1), len_dim(2), len_dim(3), len_dim(4)))
-    if (all(len_dim /= 0)) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+    if (all(len_dim /= 0)) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_real_4d
@@ -182,11 +165,42 @@
   !************************************
 
+  subroutine nf95_gw_var_real_5d(ncid, varid, values)
+
+    ! Real type, the array has rank 5.
+
+    integer, intent(in):: ncid
+    integer, intent(in):: varid
+    real, pointer:: values(:, :, :, :, :)
+
+    ! Variables local to the procedure:
+    integer len_dim(5), i
+    integer, pointer:: dimids(:)
+
+    !---------------------
+
+    call nf95_inquire_variable(ncid, varid, dimids=dimids)
+
+    if (size(dimids) /= 5) then
+       print *, "nf95_gw_var_real_5d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 5"
+       stop 1
+    end if
+
+    do i = 1, 5
+       call nf95_inquire_dimension(ncid, dimids(i), nclen=len_dim(i))
+    end do
+    deallocate(dimids) ! pointer
+
+    allocate(values(len_dim(1), len_dim(2), len_dim(3), len_dim(4), len_dim(5)))
+    if (all(len_dim /= 0)) call NF95_GET_VAR(ncid, varid, values)
+
+  end subroutine nf95_gw_var_real_5d
+
+  !************************************
+
 !!$  subroutine nf95_gw_var_dble_1d(ncid, varid, values)
 !!$
 !!$    ! Double precision, the array has rank 1.
-!!$
-!!$    use netcdf, only: NF90_GET_VAR
-!!$    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-!!$    use handle_err_m, only: handle_err
 !!$
 !!$    integer, intent(in):: ncid
@@ -195,6 +209,6 @@
 !!$
 !!$    ! Variables local to the procedure:
-!!$    integer ierr, len
-!!$    integer, pointer :: dimids(:)
+!!$    integer nclen
+!!$    integer, pointer:: dimids(:)
 !!$
 !!$    !---------------------
@@ -203,16 +217,15 @@
 !!$
 !!$    if (size(dimids) /= 1) then
-!!$       print *, "nf95_gw_var_dble_1d: NetCDF variable is not of rank 1"
-!!$       stop 1
+!!$       print *, "nf95_gw_var_dble_1d:"
+!!$       print *, "varid = ", varid
+!!$       print *, "rank of NetCDF variable is ", size(dimids), ", not 1"
+!!$        stop 1
 !!$    end if
 !!$
-!!$    call nf95_inquire_dimension(ncid, dimids(1), len=len)
+!!$    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen)
 !!$    deallocate(dimids) ! pointer
 !!$
-!!$    allocate(values(len))
-!!$    if (len /= 0) then
-!!$       ierr = NF90_GET_VAR(ncid, varid, values)
-!!$       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-!!$    end if
+!!$    allocate(values(nclen))
+!!$    if (nclen /= 0) call NF95_GET_VAR(ncid, varid, values)
 !!$
 !!$  end subroutine nf95_gw_var_dble_1d
@@ -223,8 +236,4 @@
 !!$
 !!$    ! Double precision, the array has rank 3.
-!!$
-!!$    use netcdf, only: NF90_GET_VAR
-!!$    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-!!$    use handle_err_m, only: handle_err
 !!$
 !!$    integer, intent(in):: ncid
@@ -233,6 +242,6 @@
 !!$
 !!$    ! Variables local to the procedure:
-!!$    integer ierr, len1, len2, len3
-!!$    integer, pointer :: dimids(:)
+!!$    integer nclen1, nclen2, nclen3
+!!$    integer, pointer:: dimids(:)
 !!$
 !!$    !---------------------
@@ -241,21 +250,20 @@
 !!$
 !!$    if (size(dimids) /= 3) then
-!!$       print *, "nf95_gw_var_dble_3d: NetCDF variable is not of rank 3"
+!!$       print *, "nf95_gw_var_dble_3d:"
+!!$       print *, "varid = ", varid
+!!$       print *, "rank of NetCDF variable is ", size(dimids), ", not 3"
 !!$       stop 1
 !!$    end if
 !!$
-!!$    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
-!!$    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
-!!$    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
+!!$    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen1)
+!!$    call nf95_inquire_dimension(ncid, dimids(2), nclen=nclen2)
+!!$    call nf95_inquire_dimension(ncid, dimids(3), nclen=nclen3)
 !!$    deallocate(dimids) ! pointer
 !!$
-!!$    allocate(values(len1, len2, len3))
-!!$    if (len1 * len2 * len3 /= 0) then
-!!$       ierr = NF90_GET_VAR(ncid, varid, values)
-!!$       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-!!$    end if
+!!$    allocate(values(nclen1, nclen2, nclen3))
+!!$    if (nclen1 * nclen2 * nclen3 /= 0) call NF95_GET_VAR(ncid, varid, values)
 !!$
 !!$  end subroutine nf95_gw_var_dble_3d
-
+!!$
   !************************************
 
@@ -264,8 +272,4 @@
     ! Integer type, the array has rank 1.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -273,6 +277,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len
-    integer, pointer :: dimids(:)
+    integer nclen
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -281,16 +285,15 @@
 
     if (size(dimids) /= 1) then
-       print *, "nf95_gw_var_int_1d: NetCDF variable is not of rank 1"
-       stop 1
-    end if
-
-    call nf95_inquire_dimension(ncid, dimids(1), len=len)
-    deallocate(dimids) ! pointer
-
-    allocate(values(len))
-    if (len /= 0) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+       print *, "nf95_gw_var_int_1d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 1"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen)
+    deallocate(dimids) ! pointer
+
+    allocate(values(nclen))
+    if (nclen /= 0) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_int_1d
@@ -302,8 +305,4 @@
     ! Integer type, the array has rank 3.
 
-    use netcdf, only: NF90_GET_VAR
-    use simple, only: nf95_inquire_variable, nf95_inquire_dimension
-    use handle_err_m, only: handle_err
-
     integer, intent(in):: ncid
     integer, intent(in):: varid
@@ -311,6 +310,6 @@
 
     ! Variables local to the procedure:
-    integer ierr, len1, len2, len3
-    integer, pointer :: dimids(:)
+    integer nclen1, nclen2, nclen3
+    integer, pointer:: dimids(:)
 
     !---------------------
@@ -319,18 +318,17 @@
 
     if (size(dimids) /= 3) then
-       print *, "nf95_gw_var_int_3d: NetCDF variable is not of rank 3"
-       stop 1
-    end if
-
-    call nf95_inquire_dimension(ncid, dimids(1), len=len1)
-    call nf95_inquire_dimension(ncid, dimids(2), len=len2)
-    call nf95_inquire_dimension(ncid, dimids(3), len=len3)
-    deallocate(dimids) ! pointer
-
-    allocate(values(len1, len2, len3))
-    if (len1 * len2 * len3 /= 0) then
-       ierr = NF90_GET_VAR(ncid, varid, values)
-       call handle_err("NF90_GET_VAR", ierr, ncid, varid)
-    end if
+       print *, "nf95_gw_var_int_3d:"
+       print *, "varid = ", varid
+       print *, "rank of NetCDF variable is ", size(dimids), ", not 3"
+       stop 1
+    end if
+
+    call nf95_inquire_dimension(ncid, dimids(1), nclen=nclen1)
+    call nf95_inquire_dimension(ncid, dimids(2), nclen=nclen2)
+    call nf95_inquire_dimension(ncid, dimids(3), nclen=nclen3)
+    deallocate(dimids) ! pointer
+
+    allocate(values(nclen1, nclen2, nclen3))
+    if (nclen1 * nclen2 * nclen3 /= 0) call NF95_GET_VAR(ncid, varid, values)
 
   end subroutine nf95_gw_var_int_3d
Index: trunk/LMDZ.COMMON/libf/bibio/nf95_put_var_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/nf95_put_var_m.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/bibio/nf95_put_var_m.F90	(revision 776)
@@ -9,8 +9,4 @@
           nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
           nf95_put_var_4D_FourByteReal
-!!$     module procedure nf95_put_var_1D_FourByteReal, &
-!!$          nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
-!!$          nf95_put_var_4D_FourByteReal, nf90_put_var_1D_EightByteReal, &
-!!$          nf90_put_var_3D_EightByteReal
   end interface
 
@@ -25,7 +21,7 @@
     use handle_err_m, only: handle_err
 
-    integer, intent( in) :: ncid, varid
-    real, intent( in) :: values
-    integer, dimension(:), optional, intent( in) :: start
+    integer, intent(in) :: ncid, varid
+    real, intent(in) :: values
+    integer, dimension(:), optional, intent(in) :: start
     integer, intent(out), optional:: ncerr
 
@@ -52,7 +48,7 @@
     use handle_err_m, only: handle_err
 
-    integer, intent( in) :: ncid, varid
-    integer, intent( in) :: values
-    integer, dimension(:), optional, intent( in) :: start
+    integer, intent(in) :: ncid, varid
+    integer, intent(in) :: values
+    integer, dimension(:), optional, intent(in) :: start
     integer, intent(out), optional:: ncerr
 
@@ -74,6 +70,6 @@
   !***********************
 
-  subroutine nf95_put_var_1D_FourByteReal(ncid, varid, values, start, count, &
-       stride, map, ncerr)
+  subroutine nf95_put_var_1D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
 
     use netcdf, only: nf90_put_var
@@ -82,14 +78,14 @@
     integer,                         intent(in) :: ncid, varid
     real, intent(in) :: values(:)
-    integer, dimension(:), optional, intent(in) :: start, count, stride, map
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-         map)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -103,6 +99,6 @@
   !***********************
 
-  subroutine nf95_put_var_1D_FourByteInt(ncid, varid, values, start, count, &
-       stride, map, ncerr)
+  subroutine nf95_put_var_1D_FourByteInt(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
 
     use netcdf, only: nf90_put_var
@@ -111,14 +107,14 @@
     integer,                         intent(in) :: ncid, varid
     integer, intent(in) :: values(:)
-    integer, dimension(:), optional, intent(in) :: start, count, stride, map
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-         map)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -132,22 +128,52 @@
   !***********************
 
-  subroutine nf95_put_var_2D_FourByteReal(ncid, varid, values, start, count, &
-       stride, map, ncerr)
-
-    use netcdf, only: nf90_put_var
-    use handle_err_m, only: handle_err
-
-    integer,                         intent( in) :: ncid, varid
-    real, intent( in) :: values(:, :)
-    integer, dimension(:), optional, intent( in) :: start, count, stride, map
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-         map)
+  subroutine nf95_put_var_1D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: eightByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal),     intent(in) :: values(:)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_1D_eightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_1D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_2D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real, intent(in) :: values(:, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -161,22 +187,52 @@
   !***********************
 
-  subroutine nf95_put_var_3D_FourByteReal(ncid, varid, values, start, count, &
-       stride, map, ncerr)
-
-    use netcdf, only: nf90_put_var
-    use handle_err_m, only: handle_err
-
-    integer,                         intent( in) :: ncid, varid
-    real, intent( in) :: values(:, :, :)
-    integer, dimension(:), optional, intent( in) :: start, count, stride, map
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-         map)
+  subroutine nf95_put_var_2D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: EightByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal), intent(in) :: values(:, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_2D_EightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_2D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_3D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real, intent(in) :: values(:, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -190,22 +246,52 @@
   !***********************
 
-  subroutine nf95_put_var_4D_FourByteReal(ncid, varid, values, start, count, &
-       stride, map, ncerr)
-
-    use netcdf, only: nf90_put_var
-    use handle_err_m, only: handle_err
-
-    integer,                         intent( in) :: ncid, varid
-    real, intent( in) :: values(:, :, :, :)
-    integer, dimension(:), optional, intent( in) :: start, count, stride, map
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-         map)
+  subroutine nf95_put_var_3D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: eightByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real (kind = EightByteReal),     intent(in) :: values(:, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_3D_eightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_3D_EightByteReal
+
+  !***********************
+
+  subroutine nf95_put_var_4D_FourByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer,                         intent(in) :: ncid, varid
+    real, intent(in) :: values(:, :, :, :)
+    integer, dimension(:), optional, intent(in) :: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -219,61 +305,31 @@
   !***********************
 
-!!$  subroutine nf90_put_var_1D_EightByteReal(ncid, varid, values, start, count, &
-!!$       stride, map, ncerr)
-!!$
-!!$    use typesizes, only: eightByteReal
-!!$    use netcdf, only: nf90_put_var
-!!$    use handle_err_m, only: handle_err
-!!$
-!!$    integer,                         intent( in) :: ncid, varid
-!!$    real (kind = EightByteReal),     intent( in) :: values(:)
-!!$    integer, dimension(:), optional, intent( in) :: start, count, stride, map
-!!$    integer, intent(out), optional:: ncerr
-!!$
-!!$    ! Variable local to the procedure:
-!!$    integer ncerr_not_opt
-!!$
-!!$    !-------------------
-!!$
-!!$    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-!!$         map)
-!!$    if (present(ncerr)) then
-!!$       ncerr = ncerr_not_opt
-!!$    else
-!!$       call handle_err("nf95_put_var_1D_eightByteReal", ncerr_not_opt, ncid, &
-!!$            varid)
-!!$    end if
-!!$
-!!$  end subroutine nf90_put_var_1D_EightByteReal
-!!$
-!!$  !***********************
-!!$
-!!$  subroutine nf90_put_var_3D_EightByteReal(ncid, varid, values, start, count, &
-!!$       stride, map, ncerr)
-!!$
-!!$    use typesizes, only: eightByteReal
-!!$    use netcdf, only: nf90_put_var
-!!$    use handle_err_m, only: handle_err
-!!$
-!!$    integer,                         intent( in) :: ncid, varid
-!!$    real (kind = EightByteReal),     intent( in) :: values(:, :, :)
-!!$    integer, dimension(:), optional, intent( in) :: start, count, stride, map
-!!$    integer, intent(out), optional:: ncerr
-!!$
-!!$    ! Variable local to the procedure:
-!!$    integer ncerr_not_opt
-!!$
-!!$    !-------------------
-!!$
-!!$    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
-!!$         map)
-!!$    if (present(ncerr)) then
-!!$       ncerr = ncerr_not_opt
-!!$    else
-!!$       call handle_err("nf95_put_var_3D_eightByteReal", ncerr_not_opt, ncid, &
-!!$            varid)
-!!$    end if
-!!$
-!!$  end subroutine nf90_put_var_3D_EightByteReal
+  subroutine nf95_put_var_4D_EightByteReal(ncid, varid, values, start, &
+       count_nc, stride, map, ncerr)
+
+    use typesizes, only: EightByteReal
+    use netcdf, only: nf90_put_var
+    use handle_err_m, only: handle_err
+
+    integer, intent(in):: ncid, varid
+    real(kind = EightByteReal), intent(in):: values(:, :, :, :)
+    integer, dimension(:), optional, intent(in):: start, count_nc, stride, map
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count_nc, &
+         stride, map)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_put_var_4D_EightByteReal", ncerr_not_opt, ncid, &
+            varid)
+    end if
+
+  end subroutine nf95_put_var_4D_EightByteReal
 
 end module nf95_put_var_m
Index: trunk/LMDZ.COMMON/libf/bibio/simple.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/simple.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/bibio/simple.F90	(revision 776)
@@ -2,6 +2,10 @@
 module simple
 
+  use handle_err_m, only: handle_err
+  
   implicit none
 
+  private handle_err
+
 contains
 
@@ -9,5 +13,4 @@
 
     use netcdf, only: nf90_open
-    use handle_err_m, only: handle_err
 
     character(len=*), intent(in):: path
@@ -36,8 +39,7 @@
 
     use netcdf, only: nf90_inq_dimid
-    use handle_err_m, only: handle_err
-
-    integer,             intent( in) :: ncid
-    character (len = *), intent( in) :: name
+
+    integer,             intent(in) :: ncid
+    character (len = *), intent(in) :: name
     integer,             intent(out) :: dimid
     integer, intent(out), optional:: ncerr
@@ -52,5 +54,5 @@
        ncerr = ncerr_not_opt
     else
-       call handle_err("nf95_inq_dimid", ncerr_not_opt, ncid)
+       call handle_err("nf95_inq_dimid " // name, ncerr_not_opt, ncid)
     end if
 
@@ -59,20 +61,19 @@
   !************************
 
-  subroutine nf95_inquire_dimension(ncid, dimid, name, len, ncerr)
+  subroutine nf95_inquire_dimension(ncid, dimid, name, nclen, ncerr)
 
     use netcdf, only: nf90_inquire_dimension
-    use handle_err_m, only: handle_err
 
     integer,                       intent( in) :: ncid, dimid
     character (len = *), optional, intent(out) :: name
-    integer,             optional, intent(out) :: len
-    integer, intent(out), optional:: ncerr
-
-    ! Variable local to the procedure:
-    integer ncerr_not_opt
-
-    !-------------------
-
-    ncerr_not_opt = nf90_inquire_dimension(ncid, dimid, name, len)
+    integer,             optional, intent(out) :: nclen
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inquire_dimension(ncid, dimid, name, nclen)
     if (present(ncerr)) then
        ncerr = ncerr_not_opt
@@ -88,8 +89,7 @@
 
     use netcdf, only: nf90_inq_varid
-    use handle_err_m, only: handle_err
 
     integer,             intent(in) :: ncid
-    character (len = *), intent(in) :: name
+    character(len=*), intent(in):: name
     integer,             intent(out) :: varid
     integer, intent(out), optional:: ncerr
@@ -115,5 +115,6 @@
 
     ! In "nf90_inquire_variable", "dimids" is an assumed-size array.
-    ! This is the classical case of an array the size of which is
+    ! This is not optimal.
+    ! We are in the classical case of an array the size of which is
     ! unknown in the calling procedure, before the call.
     ! Here we use a better solution: a pointer argument array.
@@ -121,5 +122,4 @@
 
     use netcdf, only: nf90_inquire_variable, nf90_max_var_dims
-    use handle_err_m, only: handle_err
 
     integer, intent(in):: ncid, varid
@@ -151,5 +151,5 @@
        ncerr = ncerr_not_opt
     else
-       call handle_err("nf95_inquire_variable", ncerr_not_opt, ncid)
+       call handle_err("nf95_inquire_variable", ncerr_not_opt, ncid, varid)
     end if
 
@@ -161,5 +161,4 @@
     
     use netcdf, only: nf90_create
-    use handle_err_m, only: handle_err
 
     character (len = *), intent(in   ) :: path
@@ -186,12 +185,11 @@
   !************************
 
-  subroutine nf95_def_dim(ncid, name, len, dimid, ncerr)
+  subroutine nf95_def_dim(ncid, name, nclen, dimid, ncerr)
 
     use netcdf, only: nf90_def_dim
-    use handle_err_m, only: handle_err
 
     integer,             intent( in) :: ncid
     character (len = *), intent( in) :: name
-    integer,             intent( in) :: len
+    integer,             intent( in) :: nclen
     integer,             intent(out) :: dimid
     integer, intent(out), optional :: ncerr
@@ -202,9 +200,9 @@
     !-------------------
 
-    ncerr_not_opt = nf90_def_dim(ncid, name, len, dimid)
-    if (present(ncerr)) then
-       ncerr = ncerr_not_opt
-    else
-       call handle_err("nf95_def_dim", ncerr_not_opt, ncid)
+    ncerr_not_opt = nf90_def_dim(ncid, name, nclen, dimid)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_def_dim " // name, ncerr_not_opt, ncid)
     end if
 
@@ -216,5 +214,4 @@
 
     use netcdf, only: nf90_redef
-    use handle_err_m, only: handle_err
 
     integer, intent( in) :: ncid
@@ -240,5 +237,4 @@
 
     use netcdf, only: nf90_enddef
-    use handle_err_m, only: handle_err
 
     integer,           intent( in) :: ncid
@@ -265,5 +261,4 @@
 
     use netcdf, only: nf90_close
-    use handle_err_m, only: handle_err
 
     integer, intent( in) :: ncid
@@ -289,5 +284,4 @@
 
     use netcdf, only: nf90_copy_att
-    use handle_err_m, only: handle_err
 
     integer, intent( in):: ncid_in,  varid_in
@@ -305,8 +299,63 @@
        ncerr = ncerr_not_opt
     else
-       call handle_err("nf95_copy_att", ncerr_not_opt, ncid_out)
+       call handle_err("nf95_copy_att " // name, ncerr_not_opt, ncid_out)
     end if
 
   end subroutine nf95_copy_att
 
+  !***********************
+
+  subroutine nf95_inquire_attribute(ncid, varid, name, xtype, nclen, attnum, &
+       ncerr)
+
+    use netcdf, only: nf90_inquire_attribute
+
+    integer,             intent( in)           :: ncid, varid
+    character (len = *), intent( in)           :: name
+    integer,             intent(out), optional :: xtype, nclen, attnum
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inquire_attribute(ncid, varid, name, xtype, nclen, &
+         attnum)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inquire_attribute " // name, ncerr_not_opt, &
+            ncid, varid)
+    end if
+
+  end subroutine nf95_inquire_attribute
+
+  !***********************
+
+  subroutine nf95_inquire(ncid, nDimensions, nVariables, nAttributes, &
+       unlimitedDimId, formatNum, ncerr)
+
+    use netcdf, only: nf90_inquire
+
+    integer,           intent( in) :: ncid
+    integer, optional, intent(out) :: nDimensions, nVariables, nAttributes
+    integer, optional, intent(out) :: unlimitedDimId, formatNum
+    integer, intent(out), optional:: ncerr
+
+    ! Variable local to the procedure:
+    integer ncerr_not_opt
+
+    !-------------------
+
+    ncerr_not_opt = nf90_inquire(ncid, nDimensions, nVariables, nAttributes, &
+         unlimitedDimId, formatNum)
+    if (present(ncerr)) then
+       ncerr = ncerr_not_opt
+    else
+       call handle_err("nf95_inquire", ncerr_not_opt, ncid)
+    end if
+
+  end subroutine nf95_inquire
+
 end module simple
Index: trunk/LMDZ.COMMON/libf/bibio/writedynav.F
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/writedynav.F	(revision 617)
+++ 	(revision )
@@ -1,152 +1,0 @@
-!
-! $Id: writedynav.F 1403 2010-07-01 09:02:53Z fairhead $
-!
-      subroutine writedynav(time, vcov, 
-     ,                ucov,teta,ppk,phi,q,masse,ps,phis)
-
-#ifdef CPP_IOIPSL
-      USE ioipsl
-#endif
-      USE infotrac, ONLY : nqtot, ttext
-      use com_io_dyn_mod, only : histaveid,histvaveid,histuaveid
-      implicit none
-
-C
-C   Ecriture du fichier histoire au format IOIPSL
-C
-C   Appels succesifs des routines: histwrite
-C
-C   Entree:
-C      time: temps de l'ecriture
-C      vcov: vents v covariants
-C      ucov: vents u covariants
-C      teta: temperature potentielle
-C      phi : geopotentiel instantane
-C      q   : traceurs
-C      masse: masse
-C      ps   :pression au sol
-C      phis : geopotentiel au sol
-C      
-C
-C
-C   L. Fairhead, LMD, 03/99
-C
-C =====================================================================
-C
-C   Declarations
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "temps.h"
-#include "ener.h"
-#include "logic.h"
-#include "description.h"
-#include "serre.h"
-#include "iniprint.h"
-
-C
-C   Arguments
-C
-
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
-      REAL teta(ip1jmp1*llm),phi(ip1jmp1,llm),ppk(ip1jmp1*llm)     
-      REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
-      REAL phis(ip1jmp1)                  
-      REAL q(ip1jmp1,llm,nqtot)
-      integer time
-
-
-#ifdef CPP_IOIPSL
-! This routine needs IOIPSL to work
-C   Variables locales
-C
-      integer ndex2d(ip1jmp1),ndexu(ip1jmp1*llm),ndexv(ip1jm*llm)
-      INTEGER iq, ii, ll
-      real tm(ip1jmp1*llm)
-      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm) 
-      logical ok_sync
-      integer itau_w
-C
-C  Initialisations
-C
-      ndexu = 0
-      ndexv = 0
-      ndex2d = 0
-      ok_sync = .TRUE.
-      tm = 999.999
-      vnat = 999.999
-      unat = 999.999
-      itau_w = itau_dyn + time
-
-C Passage aux composantes naturelles du vent
-      call covnat(llm, ucov, vcov, unat, vnat)
-
-C
-C  Appels a histwrite pour l'ecriture des variables a sauvegarder
-C
-C  Vents U 
-C
-      call histwrite(histuaveid, 'u', itau_w, unat, 
-     .               iip1*jjp1*llm, ndexu)
-C
-C  Vents V
-C
-      call histwrite(histvaveid, 'v', itau_w, vnat, 
-     .               iip1*jjm*llm, ndexv)
-C
-C  Temperature potentielle moyennee
-C
-      call histwrite(histaveid, 'theta', itau_w, teta, 
-     .                iip1*jjp1*llm, ndexu)
-C
-C  Temperature moyennee
-C
-      do ii = 1, ijp1llm
-        tm(ii) = teta(ii) * ppk(ii)/cpp
-      enddo
-      call histwrite(histaveid, 'temp', itau_w, tm, 
-     .                iip1*jjp1*llm, ndexu)
-C
-C  Geopotentiel
-C
-      call histwrite(histaveid, 'phi', itau_w, phi, 
-     .                iip1*jjp1*llm, ndexu)
-C
-C  Traceurs
-C
-!        DO iq=1,nqtot
-!          call histwrite(histaveid, ttext(iq), itau_w, q(:,:,iq), 
-!     .                   iip1*jjp1*llm, ndexu)
-!        enddo
-C
-C  Masse
-C
-       call histwrite(histaveid, 'masse', itau_w, masse, 
-     $                   iip1*jjp1*llm, ndexu)
-C
-C  Pression au sol
-C
-       call histwrite(histaveid, 'ps', itau_w, ps, iip1*jjp1, ndex2d)
-C
-C  Geopotentiel au sol
-C
-!       call histwrite(histaveid,'phis',itau_w, phis,iip1*jjp1, ndex2d)
-C
-C  Fin
-C
-      if (ok_sync) then
-          call histsync(histaveid)
-          call histsync(histvaveid)
-          call histsync(histuaveid)
-      ENDIF
-
-#else
-! tell the user this routine should be run with ioipsl
-      write(lunout,*)"writedynav: Warning this routine should not be",
-     &               " used without ioipsl"
-#endif
-! of #ifdef CPP_IOIPSL
-      return
-      end
Index: trunk/LMDZ.COMMON/libf/bibio/writedynav.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/writedynav.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/bibio/writedynav.F90	(revision 776)
@@ -0,0 +1,141 @@
+! $Id: writedynav.F90 1612 2012-01-31 10:11:48Z lguez $
+
+subroutine writedynav(time, vcov, ucov, teta, ppk, phi, q, masse, ps, phis)
+
+#ifdef CPP_IOIPSL
+  USE ioipsl
+#endif
+  USE infotrac, ONLY : nqtot, ttext
+  use com_io_dyn_mod, only : histaveid, histvaveid, histuaveid
+
+  implicit none
+
+  !   Ecriture du fichier histoire au format IOIPSL
+
+  !   Appels succesifs des routines: histwrite
+
+  !   Entree:
+  !      time: temps de l'ecriture
+  !      vcov: vents v covariants
+  !      ucov: vents u covariants
+  !      teta: temperature potentielle
+  !      phi : geopotentiel instantane
+  !      q   : traceurs
+  !      masse: masse
+  !      ps   :pression au sol
+  !      phis : geopotentiel au sol
+
+  !   L. Fairhead, LMD, 03/99
+
+  !   Declarations
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "comgeom.h"
+  include "temps.h"
+  include "ener.h"
+  include "logic.h"
+  include "description.h"
+  include "serre.h"
+  include "iniprint.h"
+
+  !   Arguments
+
+  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) 
+  REAL teta(ip1jmp1*llm), phi(ip1jmp1, llm), ppk(ip1jmp1*llm)     
+  REAL ps(ip1jmp1), masse(ip1jmp1, llm)                   
+  REAL phis(ip1jmp1)                  
+  REAL q(ip1jmp1, llm, nqtot)
+  integer time
+
+#ifdef CPP_IOIPSL
+  ! This routine needs IOIPSL to work
+  !   Variables locales
+
+  integer ndex2d(ip1jmp1), ndexu(ip1jmp1*llm), ndexv(ip1jm*llm)
+  INTEGER iq, ii, ll
+  real tm(ip1jmp1*llm)
+  REAL vnat(ip1jm, llm), unat(ip1jmp1, llm) 
+  logical ok_sync
+  integer itau_w
+
+  !-----------------------------------------------------------------
+
+  !  Initialisations
+
+  ndexu = 0
+  ndexv = 0
+  ndex2d = 0
+  ok_sync = .TRUE.
+  tm = 999.999
+  vnat = 999.999
+  unat = 999.999
+  itau_w = itau_dyn + time
+
+  ! Passage aux composantes naturelles du vent
+  call covnat(llm, ucov, vcov, unat, vnat)
+
+  !  Appels a histwrite pour l'ecriture des variables a sauvegarder
+
+  !  Vents U 
+
+  call histwrite(histuaveid, 'u', itau_w, unat,  &
+       iip1*jjp1*llm, ndexu)
+
+  !  Vents V
+
+  call histwrite(histvaveid, 'v', itau_w, vnat,  &
+       iip1*jjm*llm, ndexv)
+
+  !  Temperature potentielle moyennee
+
+  call histwrite(histaveid, 'theta', itau_w, teta,  &
+       iip1*jjp1*llm, ndexu)
+
+  !  Temperature moyennee
+
+  do ii = 1, ijp1llm
+     tm(ii) = teta(ii) * ppk(ii)/cpp
+  enddo
+  call histwrite(histaveid, 'temp', itau_w, tm,  &
+       iip1*jjp1*llm, ndexu)
+
+  !  Geopotentiel
+
+  call histwrite(histaveid, 'phi', itau_w, phi,  &
+       iip1*jjp1*llm, ndexu)
+
+  !  Traceurs
+
+  !  DO iq=1, nqtot
+  !       call histwrite(histaveid, ttext(iq), itau_w, q(:, :, iq), &
+  !                   iip1*jjp1*llm, ndexu)
+  ! enddo
+
+  !  Masse
+
+  call histwrite(histaveid, 'masse', itau_w, masse,  &
+       iip1*jjp1*llm, ndexu)
+
+  !  Pression au sol
+
+  call histwrite(histaveid, 'ps', itau_w, ps, iip1*jjp1, ndex2d)
+
+  ! Geopotentiel au sol
+
+  ! call histwrite(histaveid, 'phis', itau_w, phis, iip1*jjp1, ndex2d)
+
+  if (ok_sync) then
+     call histsync(histaveid)
+     call histsync(histvaveid)
+     call histsync(histuaveid)
+  ENDIF
+
+#else
+  write(lunout, *) "writedynav: Warning this routine should not be", &
+       " used without ioipsl"
+#endif
+  ! of #ifdef CPP_IOIPSL
+
+end subroutine writedynav
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'
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/bands.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/bands.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/bands.F90	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: bands.F90 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: bands.F90 1615 2012-02-10 15:42:26Z emillour $
 !
   module Bands
@@ -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: trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 776)
@@ -27,8 +27,5 @@
      $                  pdpsfi)
 #ifdef CPP_PHYS
-! Ehouarn: For now, calfis_p needs Earth physics
-c
-c    Auteur :  P. Le Van, F. Hourdin 
-c   .........
+! Ehouarn: if using (parallelized) physics
       USE dimphy
       USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
@@ -225,5 +222,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)
       
@@ -512,6 +509,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,
@@ -519,4 +517,5 @@
      $           ntetaSTD,rtetaSTD,PVteta)
 c
+#endif
       ENDIF
 
@@ -666,4 +665,5 @@
       zdqfic_omp(:,:,:)=0.
 
+#ifdef CPP_PHYS
       do isplit=1,nsplit_phys
 
@@ -742,4 +742,7 @@
       enddo
 
+#endif
+! of #ifdef CPP_PHYS
+
       zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
       zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1600 2011-12-06 13:16:30Z jghattas $
+! $Id: ce0l.F90 1615 2012-02-10 15:42:26Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -31,5 +31,6 @@
   IMPLICIT NONE
 #ifndef CPP_EARTH
-  WRITE(*,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
+#include "iniprint.h"
+  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
 #else
 !-------------------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/comvert.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/comvert.h	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/disvert.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/disvert.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/dynetat0.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/dynetat0.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/dynredem.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/dynredem_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F	(revision 776)
@@ -1,4 +1,4 @@
 !
-! $Id: dynredem_p.F 1577 2011-10-20 15:06:47Z fairhead $
+! $Id: dynredem_p.F 1635 2012-07-12 11:37:16Z lguez $
 !
 c
@@ -9,4 +9,6 @@
       USE parallel
       USE infotrac
+      use netcdf95, only: NF95_PUT_VAR
+ 
       IMPLICIT NONE
 c=======================================================================
@@ -19,5 +21,5 @@
 #include "comconst.h"
 #include "comvert.h"
-#include "comgeom.h"
+#include "comgeom2.h"
 #include "temps.h"
 #include "ener.h"
@@ -30,5 +32,5 @@
 c   ----------
       INTEGER iday_end
-      REAL phis(ip1jmp1)
+      REAL phis(iip1, jjp1)
       CHARACTER*(*) fichnom
 
@@ -56,4 +58,5 @@
       character*30 unites
 
+
 c-----------------------------------------------------------------------
       if (mpi_rank==0) then
@@ -69,14 +72,14 @@
       mmois0=1
       jjour0=1
-#endif                
+#endif        
 
       DO l=1,length
        tab_cntrl(l) = 0.
       ENDDO
-       tab_cntrl(1)  =  REAL(iim)
-       tab_cntrl(2)  =  REAL(jjm)
-       tab_cntrl(3)  =  REAL(llm)
-       tab_cntrl(4)  =  REAL(day_ref)
-       tab_cntrl(5)  =  REAL(annee_ref)
+       tab_cntrl(1)  = REAL(iim)
+       tab_cntrl(2)  = REAL(jjm)
+       tab_cntrl(3)  = REAL(llm)
+       tab_cntrl(4)  = REAL(day_ref)
+       tab_cntrl(5)  = REAL(annee_ref)
        tab_cntrl(6)  = rad
        tab_cntrl(7)  = omeg
@@ -118,6 +121,6 @@
       ENDIF
 
-       tab_cntrl(30) =  REAL(iday_end)
-       tab_cntrl(31) =  REAL(itau_dyn + itaufin)
+       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
@@ -165,9 +168,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)
@@ -182,9 +181,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)
@@ -199,9 +194,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)
@@ -216,9 +207,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)
@@ -233,9 +220,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)
@@ -250,9 +233,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)
@@ -267,9 +246,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)
@@ -284,9 +259,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)
@@ -301,9 +272,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)
@@ -316,9 +283,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
@@ -337,9 +300,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)
@@ -356,9 +315,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:
@@ -377,9 +332,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:
@@ -398,9 +349,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:
@@ -510,5 +457,4 @@
       ierr = NF_ENDDEF(nid) ! sortir du mode de definition
       ierr = NF_CLOSE(nid) ! fermer le fichier
-
 
       PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
@@ -524,4 +470,7 @@
       USE infotrac
       USE control_mod
+      use netcdf, only: NF90_get_VAR
+      use netcdf95, only: NF95_PUT_VAR
+ 
       IMPLICIT NONE
 c=================================================================
@@ -536,9 +485,10 @@
 #include "temps.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
      
@@ -546,5 +496,5 @@
       INTEGER nid, nvarid, nid_trac, nvarid_trac
       REAL trac_tmp(ip1jmp1,llm)      
-      INTEGER ierr, ierr_file
+      INTEGER ierr, ierr_file 
       INTEGER iq
       INTEGER length
@@ -567,5 +517,5 @@
       
       do iq=1,nqtot
-        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
+        call Gather_Field(q(:,:,:,iq),ip1jmp1,llm,0)
       enddo
       
@@ -589,9 +539,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/))
       PRINT*, "Enregistrement pour ", nb, time
 
@@ -605,15 +551,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
-       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
+      ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl)
+       tab_cntrl(31) = REAL(itau_dyn + itaufin)
+      call NF95_PUT_VAR(nid,nvarid,tab_cntrl)
 
 c  Ecriture des champs
@@ -624,9 +562,5 @@
          CALL abort
       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)
@@ -635,9 +569,5 @@
          CALL abort
       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)
@@ -646,9 +576,5 @@
          CALL abort
       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
@@ -675,9 +601,5 @@
                CALL abort
             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
@@ -691,17 +613,9 @@
                    CALL abort
                 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
                 PRINT*, 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
                    PRINT*, "Lecture echouee pour", tname(iq)
@@ -713,9 +627,5 @@
                    CALL abort
                 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)
@@ -728,11 +638,7 @@
                 CALL abort
              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
+       END IF   !type_trac
       
       ENDDO
@@ -746,9 +652,5 @@
          CALL abort
       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)
@@ -757,9 +659,5 @@
          CALL abort
       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/dyn3dpar/etat0_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/etat0_netcdf.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/exner_hyb.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/exner_hyb_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/exner_milieu_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.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.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/dyn3dpar/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 776)
@@ -20,6 +20,6 @@
       USE control_mod
 
-! Ehouarn: for now these only apply to Earth:
-#ifdef CPP_EARTH
+! Ehouarn: the following are needed with (parallel) physics:
+#ifdef CPP_PHYS
       USE mod_grid_phy_lmdz
       USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
@@ -182,14 +182,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_PHYS
-! Ehouarn: NB: For now only Earth physics is parallel
       CALL Init_interface_dyn_phys
 #endif
@@ -203,7 +199,4 @@
 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
@@ -211,5 +204,5 @@
 c$OMP END PARALLEL 
 #endif
-!      endif ! of if (planet_type.eq."earth")
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 c
@@ -476,8 +469,5 @@
          WRITE(lunout,*)
      .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
-
-! Initialisation de la physique: pose probleme quand on tourne 
-! SANS physique, car iniphysiq.F est dans le repertoire phy[]... 
-! Il faut une cle CPP_PHYS
+! Physics
 #ifdef CPP_PHYS
          CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
@@ -486,5 +476,5 @@
          call_iniphys=.false.
       ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
-!#endif
+
 
 c-----------------------------------------------------------------------
@@ -525,4 +515,12 @@
 #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
+
       if (planet_type.eq."mars") then
 ! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/gr_dyn_fi_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/gr_dyn_fi_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/gr_dyn_fi_p.F	(revision 776)
@@ -1,9 +1,8 @@
 !
-! $Id: gr_dyn_fi_p.F 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: gr_dyn_fi_p.F 1615 2012-02-10 15:42:26Z emillour $
 !
       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: trunk/LMDZ.COMMON/libf/dyn3dpar/gr_fi_dyn_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/gr_fi_dyn_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/gr_fi_dyn_p.F	(revision 776)
@@ -1,9 +1,8 @@
 !
-! $Id: gr_fi_dyn_p.F 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: gr_fi_dyn_p.F 1615 2012-02-10 15:42:26Z emillour $
 !
       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: trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 776)
@@ -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: trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/iniconst.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/iniconst.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F90	(revision 776)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/inigrads.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/inigrads.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/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/dyn3dpar/integrd_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 776)
@@ -1,8 +1,8 @@
 !
-! $Id: integrd_p.F 1550 2011-07-05 09:44:55Z lguez $
+! $Id: integrd_p.F 1616 2012-02-17 11:59:00Z emillour $
 !
       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: trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 776)
@@ -132,5 +132,5 @@
       REAL  SSUM
       REAL time_0 
-      REAL,SAVE :: finvmaold(ip1jmp1,llm)
+!      REAL,SAVE :: finvmaold(ip1jmp1,llm)
 
 cym      LOGICAL  lafin
@@ -272,7 +272,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
@@ -283,10 +283,10 @@
 c et du parallelisme !!
 
-   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) 
       if (jH_cur > 1.0 ) then
         jD_cur = jD_cur +1.
@@ -324,6 +324,7 @@
          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
@@ -343,5 +344,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
@@ -353,6 +354,7 @@
 
 
-          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)
@@ -370,5 +372,5 @@
 cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
-   2  CONTINUE
+   2  CONTINUE ! Matsuno backward or leapfrog step begins here
 
 c$OMP MASTER
@@ -515,6 +517,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
@@ -616,4 +618,5 @@
       call start_timer(timer_caldyn)
 
+      ! compute geopotential phi()
 ! ADAPTATION GCM POUR CP(T)
 !      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
@@ -699,6 +702,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                                    )
 
        IF ((planet_type.eq."titan").and.(tidal)) then
@@ -773,14 +776,14 @@
 
 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) 
+     $        + itau/day_step
            jH_cur = jH_ref + start_time +                                &
-     &              (itau * dtvr / daysec - int(itau * dtvr / daysec)) 
+     &              mod(itau,day_step)/float(day_step) 
 !         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
            if (jH_cur > 1.0 ) then
@@ -803,6 +806,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 
@@ -1131,7 +1136,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
@@ -1295,4 +1300,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
@@ -1305,5 +1314,6 @@
           ENDDO
 c$OMP END MASTER
-        endif
+         endif ! of if (1 == 0)
+        endif ! of of (pole_nord)
         
         if (pole_sud) then
@@ -1321,4 +1331,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
@@ -1331,5 +1345,6 @@
           ENDDO
 c$OMP END MASTER
-        endif
+         endif ! of if (1 == 0)
+        endif ! of if (pole_sud)
 
 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/mod_interface_dyn_phys.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/mod_interface_dyn_phys.F90	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/mod_interface_dyn_phys.F90	(revision 776)
@@ -1,4 +1,4 @@
 ! 
-! $Id: mod_interface_dyn_phys.F90 1279 2009-12-10 09:02:56Z fairhead $
+! $Id: mod_interface_dyn_phys.F90 1615 2012-02-10 15:42:26Z emillour $
 !
 MODULE mod_interface_dyn_phys
@@ -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: trunk/LMDZ.COMMON/libf/dyn3dpar/wrgrads.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/wrgrads.F	(revision 617)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/wrgrads.F	(revision 776)
@@ -17,4 +17,8 @@
       integer if,nl
       real field(imx*jmx*lmx)
+
+      integer, parameter:: wp = selected_real_kind(p=6, r=36)
+      real(wp) field4(imx*jmx*lmx)
+
       character*10 name,file
       character*10 titlevar
@@ -22,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
@@ -29,5 +33,5 @@
       writectl=.false.
 
-      print*,if,iid(if),jid(if),ifd(if),jfd(if)
+c     print*,if,iid(if),jid(if),ifd(if),jfd(if)
       iii=iid(if)
       iji=jid(if)
@@ -38,6 +42,6 @@
       lm=lmd(if)
 
-      print*,'im,jm,lm,name,firsttime(if)'
-      print*,im,jm,lm,name,firsttime(if)
+c     print*,'im,jm,lm,name,firsttime(if)'
+c     print*,im,jm,lm,name,firsttime(if)
 
       if(firsttime(if)) then
@@ -55,8 +59,8 @@
             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)
-            print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
+c           print*,'initialisation ecriture de ',var(ivar(if),if)
+c           print*,'if ivar(if) nld ',if,ivar(if),nld(ivar(if),if)
          endif
          writectl=.true.
@@ -81,6 +85,7 @@
       endif
 
-      print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
-      print*,ivar(if),nvar(if),var(ivar(if),if),writectl
+c     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
+c     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
+      field4(1:imd(if)*jmd(if)*nl)=field(1:imd(if)*jmd(if)*nl)
       do l=1,nl
          irec(if)=irec(if)+1
@@ -89,5 +94,5 @@
 c    s ,(l-1)*imd(if)*jmd(if)+(ijf-1)*imd(if)+iif
          write(unit(if)+1,rec=irec(if))
-     s   ((field((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i)
+     s   ((field4((l-1)*imd(if)*jmd(if)+(j-1)*imd(if)+i)
      s   ,i=iii,iif),j=iji,ijf)
       enddo
@@ -96,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'
