Index: /trunk/DOC/chantiers/commit_importants.log
===================================================================
--- /trunk/DOC/chantiers/commit_importants.log	(revision 1301)
+++ /trunk/DOC/chantiers/commit_importants.log	(revision 1302)
@@ -1389,2 +1389,63 @@
 - guide_mod.F : upate to keep up with Earth model
 
+**********************
+**** commit_v1302 ****
+**********************
+Ehouarn: Updates in the dynamics (seq and //) to keep up with updates
+in LMDZ5 (up to LMDZ5 trunk, rev 2070):
+* compilation:
+- bld.cfg : adaptations for compilation with RRTM and PRISM.
+- makelmdz_fcm : adaptations for compilation with RRTM and OASIS-MCT ocean coupled mode
+  and change handling of RRTM and SISVAT so that bld.cfg is never overwritten (as is the case
+  for COSP), since bld.cfg is sometimes a link to the one located in LMDZ.COMMON.
+- makelmdz : adaptations for compilation with RRTM and in 1D, and some cleanup.
+- create_make_gcm : cleanup and improved scaning of dependencies 
+
+* dyn3d_common:
+- iniconst.F90 : change argument list in call to disvert
+- disvert.F90 : change argument list; add "sima", "starto_**" cases
+- grilles_gcm_netcdf.F : removed because obsolete
+- q_sat.F : made free-form compatible
+- exner_hyb_m.F90 : now a module; some cleanup around beta and alpha;
+  made pkf an optional argument and removed the averaging of pks at the poles (because was
+  is not necessary; note that it changes model results wrt previous versions of code)
+- exner_milieu_m.F90 : now a module; some cleanup around beta and alpha;
+  made pkf an optional argument and removed the averaging of pks at the poles (because was
+  is not necessary; note that it changes model results wrt previous versions of code)
+- iniacademic.F90 : adaptations for changes with exner_* routines
+- removed PVtheta.F and tourabs.F
+
+* dyn3dpar :
+- parallel_lmdz.F90 : updates for XIOS and Oasis-MCT coupler
+- mod_const_mpi.F90 : updates for XIOS and Oasis-MCT coupler
+- conf_gcm.F : remove obsolete "grilles_gcm_netcdf" flag
+- logic.h : remove obsolete "grilles_gcm_netcdf" flag
+- ce0l.F : follow up on updates from Earth model
+- leapfrog_p.F : enable runing a number of dyn steps given by a negative "nday" value. Adaptations
+  for exner_* routines. Add computation of geopotential before calling physics.
+- gcm.F : enable runing a number of dyn steps given by a negative "nday" value
+- guide_p_mod.F90 : follow up on updates from Earth model
+- exner_hyb_p_m.F90 : now a module; some cleanup around beta and alpha;
+  made pkf an optional argument and removed the averaging of pks at the poles (because was
+  is not necessary; note that it changes model results wrt previous versions of code)
+- exner_milieu_p_m.F90 : now a module; some cleanup around beta and alpha;
+  made pkf an optional argument and removed the averaging of pks at the poles (because was
+  is not necessary; note that it changes model results wrt previous versions of code)
+- calfis_p.F : removed call to PVtheta
+
+* dyn3d :
+- conf_gcm.F : remove obsolete "grilles_gcm_netcdf" flag and made free-form compatible
+  (because it is included in the 1D Earth model, which is free form)
+- logic.h : remove obsolete "grilles_gcm_netcdf" flag
+- ce0l.F : follow up on updates from Earth model
+- leapfrog.F : enable runing a number of dyn steps given by a negative "nday" value. Adaptations
+  for exner_* routines. Add computation of geopotential before calling physics.
+- gcm.F : enable runing a number of dyn steps given by a negative "nday" value
+- guide_mod.F90 : follow up on updates from Earth model
+- calfis.F : removed call to PVtheta
+
+* bibio :
+- pres2lev_mod.F90 : move from dyn3d_common to bibio
+- wxios.F90 : follow up on updates on XIOS usage
+
+
Index: /trunk/LMDZ.COMMON/arch/arch-X64_ADA.path
===================================================================
--- /trunk/LMDZ.COMMON/arch/arch-X64_ADA.path	(revision 1301)
+++ /trunk/LMDZ.COMMON/arch/arch-X64_ADA.path	(revision 1302)
@@ -3,3 +3,11 @@
 IOIPSL_INCDIR=/smphome/rech/psl/rpsl035/IOIPSL_PLUS/modipsl/lib
 IOIPSL_LIBDIR=/smphome/rech/psl/rpsl035/IOIPSL_PLUS/modipsl/lib
+XIOS_INCDIR=$LMDGCM/../XIOS/inc
+XIOS_LIBDIR=$LMDGCM/../XIOS/lib
+ORCH_INCDIR=$LMDGCM/../../lib
+ORCH_LIBDIR=$LMDGCM/../../lib
+OASIS_INCDIR=$LMDGCM/../../oasis3-mct/BLD/build/lib/psmile.MPI1
+OASIS_LIBDIR=$LMDGCM/../../oasis3-mct/BLD/lib
+INCA_LIBDIR=$LMDGCM/../INCA/config/lib
+INCA_INCDIR=$LMDGCM/../INCA/config/lib
 
Index: /trunk/LMDZ.COMMON/bld.cfg
===================================================================
--- /trunk/LMDZ.COMMON/bld.cfg	(revision 1301)
+++ /trunk/LMDZ.COMMON/bld.cfg	(revision 1302)
@@ -26,5 +26,6 @@
 src::dyn     %SRC_PATH/%DYN
 src::phys    %SRC_PATH/%PHYS
-#src::sisvat  %SRC_PATH/%PHYS/sisvat
+src::sisvat  %SISVAT
+src::rrtm    %RRTM
 src::grid    %SRC_PATH/grid
 src::filtrez %SRC_PATH/filtrez
@@ -88,4 +89,5 @@
 bld::excl_dep        use::mod_prism_get_proto
 bld::excl_dep        use::mod_prism_put_proto
+bld::excl_dep        use::mod_prism
 bld::excl_dep        use::xios
 bld::excl_dep        use::iaxis
Index: /trunk/LMDZ.COMMON/create_make_gcm
===================================================================
--- /trunk/LMDZ.COMMON/create_make_gcm	(revision 1301)
+++ /trunk/LMDZ.COMMON/create_make_gcm	(revision 1302)
@@ -5,9 +5,10 @@
 #set -xv
 
-rm -f create_make_gcm.log
 
 machine=`hostname`
 os=`uname`
 gcm=`pwd`
+logfile=$gcm/create_make_gcm.log
+rm -f $logfile ; touch -f $logfile
 libf=$gcm/libf
 libo=$gcm/libo
@@ -51,5 +52,4 @@
 echo "FILTRE   = filtre"
 echo "PHYS  = "
-echo "DYN  = dyn "
 echo 'LIBPHY = $(LIBO)/libphy$(PHYS).a'
 echo 'DIRMAIN=dyn$(DIM)d$(FLAG_PARA)'
@@ -59,14 +59,8 @@
 echo
 echo '# Les differentes librairies pour l"edition des liens:'
-echo
-echo 'dyn3d            = $(LIBO)/libdyn3d.a $(LIBO)/lib$(FILTRE).a'
-echo 'dyn3dpar      = $(LIBO)/libdyn3dpar.a $(LIBO)/lib$(FILTRE).a'
-echo 'dyn2d            = $(LIBO)/libdyn2d.a'
-echo 'dyn1d            = $(LIBO)/libdyn1d.a'
-echo 'L_DYN      = -ldyn$(DIM)d$(FLAG_PARA)'
+echo 'L_DYN      = -ldyn$(DIM)d$(FLAG_PARA) -ldyn3d_common'
 echo 'L_FILTRE   = -l$(FILTRE)'
 echo 'L_PHY = -lphy$(PHYS) '
 echo 'L_BIBIO    = -lbibio'
-echo 'L_DYN3D_COMMON    = -ldyn3d_common'
 echo 'L_ADJNT    ='
 echo 'L_COSP     = -lcosp'
@@ -100,7 +94,8 @@
 echo
 echo "PROG = code"
+echo 'DYN  = dyn$(DIM)d'
 echo
 #echo 'main : chimie $(DYN) bibio phys $(OPTION_DEP) '
-echo 'main : $(DYN) bibio dyn3d_common phys $(OPTION_DEP) '
+echo 'main : $(DYN) bibio phys $(OPTION_DEP) '
 echo '	cd $(LIBO) ; $(RANLIB) lib*.a ; cd $(GCM) ;\'
 echo '	cd $(LOCAL_DIR); \'
@@ -108,13 +103,13 @@
 echo '	$(LINK) $(PROG).o -L$(LIBO) $(L_DYN) $(L_ADJNT) $(L_COSP) $(L_AERONOMARS) $(L_PHY) $(L_DYN) $(L_DYN3D_COMMON) $(L_BIBIO) $(L_DYN3D_COMMON) $(L_PHY) $(L_DYN) $(L_FILTRE) $(OPLINK) $(OPTION_LINK) -o $(LOCAL_DIR)/$(PROG).e ; $(RM) $(PROG).o '
 echo
-echo 'dyn : $(LIBO)/libdyn$(DIM)d$(FLAG_PARA).a $(FILTRE)$(DIM)d'
-echo
 echo 'phys : $(LIBPHY)'
+echo 'dyn3d : $(LIBO)/libdyn$(DIM)d$(FLAG_PARA).a $(LIBO)/lib$(FILTRE).a $(LIBO)/libdyn3d_common.a'
+echo 'dyn2d : $(LIBO)/libdyn2d.a'
+echo 'dyn1d :'
+echo
 echo
 #echo 'chimie : $(LIBO)/libchimie.a'
 echo
 echo 'bibio : $(LIBO)/libbibio.a'
-echo
-echo 'dyn3d_common : $(LIBO)/libdyn3d_common.a'
 echo
 echo 'adjnt : $(LIBO)/libadjnt.a'
@@ -138,6 +133,6 @@
 echo
 cd $libf >/dev/null 2>&1
-for diri in ` ls `
-do
+for diri in ` ls ` ; do
+
    if [ -d $diri ] ; then
    if [ "`ls $diri/*.F $diri/*/*.F`" != "" ] || [ "`ls $diri/*.F90 $diri/*/*.F90`" != "" ]  ; then 
@@ -146,6 +141,5 @@
       listlib=""
 # Liste des fichiers .F et .F90 n'etant pas des programmes principaux
-      for fili in `ls $diri/*.F $diri/*/*.F`
-      do
+      for fili in `ls $diri/*.F $diri/*/*.F` ; do
          # Check if file is a routine or main program
          # i.e: look for the "program" keword preceeded by leading spaces
@@ -156,6 +150,5 @@
          fi
       done
-      for fili in `ls $diri/*.F90 $diri/*/*.F90`
-      do
+      for fili in `ls $diri/*.F90 $diri/*/*.F90` ; do
          # Check if file is a routine or main program
          # i.e. look for "program" keyword (with possibly some leading spaces)
@@ -168,5 +161,5 @@
 #
       echo
-      echo $listlib >> create_make_gcm.log
+      echo $listlib >> $logfile
       echo
       echo "#======================================================================="
@@ -185,11 +178,5 @@
       echo "# Compilation des membres de la bibliotheque lib"$diri".a"
       echo
-      for trufile in $listlib
-      do
-#         if [ -f $fili.F90 ] ; then
-#           trufile=$fili.F90
-#         else
-#           trufile=$fili.F
-#         fi
+      for trufile in $listlib ; do
          fili=`echo $trufile | awk -F/ ' { print $NF } ' | cut -d. -f1`
          F90=0 ; egrep -i '^ *use ' $trufile > /dev/null 2>&1 && F90=1
@@ -208,10 +195,9 @@
 ##########################################################################
 
-         for stri in ` ( sed -e '/^[cC\!]/d' $trufile | sed -n "/\#include/s/\#include//p" | sed 's/\"//g' ; sed -e '/^[cC\!]/d' $trufile | egrep -i '^ *use ' | sed -e 's/,/ /' | awk ' { print $2 } ' ) `
-         do
-            echo $trufile $fili $stri >> create_make_gcm.log
+         for stri in ` ( sed -e '/^[cC\!]/d' $trufile | sed -e "/\#include /s/\#include /XXnclude /" -e "/include /s/include /XXnclude /" -e "/INCLUDE /s/INCLUDE /XXnclude /" | sed -n "/XXnclude /s/XXnclude//p" | sed 's/\"//g' | sort | uniq ; sed -e '/^[cC\!]/d' $trufile | egrep -i '^ *use ' | sed -e 's/,/ /' | awk ' { print $2 } ' | sort | uniq ) ` ; do
+            echo $trufile $fili $stri >> $logfile
             strj=`echo $stri | tr [A-Z] [a-z]`
             str2=""
-            for dirinc in filtrez bibio dyn3d_common grid dyn3d phydev $diri $diri/*/ ; do
+            for dirinc in filtrez bibio dyn3d_common grid dyn3d phy* $diri $diri/*/ ; do
 # Recherche dans l'ordre hierarchique inverse car seule la derniere
 # ligne est conservee
@@ -223,4 +209,5 @@
                    libstr=$dirinc
                fi
+               echo dirinc $dirinc >> $logfile
                if [ -f $dirinc/$stri ] ; then
                   str2='$(LIBF)/'$dirstr/$stri
@@ -234,9 +221,12 @@
                   str2='$(LIBO)/lib'$strlib'.a('$stri'.o)'
                fi
+               echo stri $stri  >> $logfile
+               echo str2 $str2  >> $logfile
             done
+               echo dernier str2 $str2  >> $logfile
             if [ "$str2" = "" ] ; then
-               echo Warn: Include $diri $trufile $fili $stri pas trouve >> $gcm/create_make_gcm.log
+               echo Warn: Include $diri $trufile $fili $stri pas trouve >> $logfile
             elif [ "$str2" = "$str1" ] ; then
-               echo Warn: Include $diri $trufile $fili $stri $str1 trouve 2x >> $gcm/create_make_gcm.log
+               echo Warn: Include $diri $trufile $fili $stri $str1 trouve 2x >> $logfile
             else
                echo $str1 \\
Index: /trunk/LMDZ.COMMON/libf/bibio/pres2lev_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/bibio/pres2lev_mod.F90	(revision 1302)
+++ /trunk/LMDZ.COMMON/libf/bibio/pres2lev_mod.F90	(revision 1302)
@@ -0,0 +1,79 @@
+! $Id: pres2lev.F 1179 2009-06-11 14:18:47Z jghattas $
+!
+MODULE pres2lev_mod
+
+CONTAINS 
+
+!******************************************************
+SUBROUTINE pres2lev(varo,varn,lmo,lmn,po,pn,ni,nj,ok_invertp)
+!
+! interpolation lineaire pour passer
+! a une nouvelle discretisation verticale pour
+! les variables de GCM
+! Francois Forget (01/1995)
+! MOdif remy roca 12/97 pour passer de pres2sig
+! Modif F.Codron 07/08 po en 3D
+!**********************************************************
+
+  IMPLICIT NONE
+
+!   Declarations:
+! ==============
+!
+!  ARGUMENTS
+!  """""""""
+  LOGICAL, INTENT(IN) :: ok_invertp
+  INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
+  INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
+  
+  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
+  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
+  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
+
+  REAL, INTENT(IN)  :: varo(ni*nj,lmo) ! var dans l'ancienne grille
+  REAL, INTENT(OUT) :: varn(ni*nj,lmn) ! var dans la nouvelle grille
+
+  REAL :: zvaro(ni*nj,lmo),zpo(ni*nj,lmo)
+
+! Autres variables
+! """"""""""""""""
+  INTEGER ::  ln ,lo, k
+  REAL    :: coef
+
+
+! Inversion de l'ordre des niveaux verticaux
+  IF (ok_invertp) THEN
+    DO lo=1,lmo
+      DO k=1,ni*nj
+        zpo(k,lo)=po(k,lmo+1-lo)
+        zvaro(k,lo)=varo(k,lmo+1-lo)
+      ENDDO
+    ENDDO
+  ELSE
+    DO lo=1,lmo
+      DO k=1,ni*nj
+        zpo(k,lo)=po(k,lo)
+        zvaro(k,lo)=varo(k,lo)
+      ENDDO
+    ENDDO
+  ENDIF 
+
+  DO ln=1,lmn
+    DO lo=1,lmo-1
+      DO k=1,ni*nj
+        IF (pn(k,ln) >= zpo(k,1) ) THEN
+          varn(k,ln) = zvaro(k,1)
+        ELSE IF (pn(k,ln) <= zpo(k,lmo)) THEN
+          varn(k,ln) = zvaro(k,lmo)
+        ELSE IF ( pn(k,ln) <= zpo(k,lo) .AND. pn(k,ln) > zpo(k,lo+1) ) THEN
+          coef = (pn(k,ln)-zpo(k,lo)) / (zpo(k,lo+1)-zpo(k,lo))
+          varn(k,ln) = zvaro(k,lo) + coef*(zvaro(k,lo+1)-zvaro(k,lo))
+        ENDIF
+         
+      ENDDO  
+    ENDDO
+  ENDDO                
+
+END SUBROUTINE pres2lev    
+
+END MODULE pres2lev_mod
Index: /trunk/LMDZ.COMMON/libf/bibio/wxios.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/bibio/wxios.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/bibio/wxios.F90	(revision 1302)
@@ -26,17 +26,4 @@
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    !   str + i   =>   str_i   !!!!!!!!!!!!!!!!!!!!
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-    SUBROUTINE concat(str, str2, str_str2)
-        CHARACTER(len=*), INTENT(IN) :: str, str2
-        CHARACTER(len=20), INTENT(OUT) :: str_str2
-        
-        
-        str_str2 = TRIM(ADJUSTL(str//"_"//TRIM(ADJUSTL(str2))))
-        !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ",str,"+",str2,"=",str_str2
-    END SUBROUTINE concat
-    
-    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !   36day => 36d etc     !!!!!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -109,5 +96,5 @@
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-    SUBROUTINE wxios_init(xios_ctx_name, locom, outcom)
+    SUBROUTINE wxios_init(xios_ctx_name, locom, outcom, type_ocean)
         IMPLICIT NONE
         INCLUDE 'iniprint.h'
@@ -116,4 +103,5 @@
       INTEGER, INTENT(IN), OPTIONAL :: locom
       INTEGER, INTENT(OUT), OPTIONAL :: outcom
+      CHARACTER(len=6), INTENT(IN), OPTIONAL :: type_ocean
 
     
@@ -142,6 +130,9 @@
         g_ctx_name = xios_ctx_name
         
-        CALL wxios_context_init()
-        
+        ! Si couple alors init fait dans cpl_init
+        IF (.not. PRESENT(type_ocean)) THEN
+            CALL wxios_context_init()
+        ENDIF
+
     END SUBROUTINE wxios_init
 
@@ -158,6 +149,8 @@
         g_ctx = xios_ctx
 
-        IF (prt_level >= 10) WRITE(lunout,*) "wxios_context_init: Current context is ",trim(g_ctx_name)
-
+        IF (prt_level >= 10) THEN
+          WRITE(lunout,*) "wxios_context_init: Current context is ",trim(g_ctx_name)
+          WRITE(lunout,*) "     now call xios_solve_inheritance()"
+        ENDIF
         !Une premiÃ¨re analyse des hÃ©ritages:
         CALL xios_solve_inheritance()
@@ -303,28 +296,34 @@
     ! Pour dÃ©clarer un axe vertical !!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    SUBROUTINE wxios_add_vaxis(axisgroup_id, axis_file, axis_size, axis_value)
-        IMPLICIT NONE
-        INCLUDE 'iniprint.h'
-
-        CHARACTER (len=*), INTENT(IN) :: axisgroup_id, axis_file
+    SUBROUTINE wxios_add_vaxis(axis_id, axis_size, axis_value)
+        IMPLICIT NONE
+        INCLUDE 'iniprint.h'
+
+        CHARACTER (len=*), INTENT(IN) :: axis_id
         INTEGER, INTENT(IN) :: axis_size
         REAL, DIMENSION(axis_size), INTENT(IN) :: axis_value
         
-        TYPE(xios_axisgroup) :: axgroup
-        TYPE(xios_axis) :: ax
-        CHARACTER(len=20) :: axis_id
-        
-        
-        !PrÃ©paration du nom de l'axe:
-        CALL concat(axisgroup_id, axis_file, axis_id)
+!        TYPE(xios_axisgroup) :: axgroup
+!        TYPE(xios_axis) :: ax
+!        CHARACTER(len=50) :: axis_id 
+        
+!        IF (len_trim(axisgroup_id).gt.len(axis_id)) THEN
+!          WRITE(lunout,*) "wxios_add_vaxis: error, size of axis_id too small!!"
+!          WRITE(lunout,*) "     increase it to at least ",len_trim(axisgroup_id)
+!          CALL abort_gcm("wxios_add_vaxis","len(axis_id) too small",1)
+!        ENDIF
+!        axis_id=trim(axisgroup_id)
         
         !On rÃ©cupÃ¨re le groupe d'axes qui va bien:
-        CALL xios_get_axisgroup_handle(axisgroup_id, axgroup)
+        !CALL xios_get_axisgroup_handle(axisgroup_id, axgroup)
         
         !On ajoute l'axe correspondant Ã  ce fichier:
-        CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id)))
+        !CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id)))
         
         !Et on le parametrise:
-        CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value)
+        !CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value)
+        
+        ! Ehouarn: New way to declare axis, without axis_group:
+        CALL xios_set_axis_attr(trim(axis_id),size=axis_size,value=axis_value)
         
         !VÃ©rification:
@@ -332,5 +331,5 @@
             IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_vaxis: Axis created: ", TRIM(ADJUSTL(axis_id))
         ELSE
-            WRITE(*,*) "wxios_add_vaxis: Invalid axis: ", TRIM(ADJUSTL(axis_id))
+            WRITE(lunout,*) "wxios_add_vaxis: Invalid axis: ", TRIM(ADJUSTL(axis_id))
         END IF
 
@@ -367,18 +366,23 @@
         
             IF (xios_is_valid_file("X"//fname)) THEN
-                IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: New file: ", "X"//fname
-                IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                IF (prt_level >= 10) THEN
+                  WRITE(lunout,*) "wxios_add_file: New file: ", "X"//fname
+                  WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                ENDIF
             ELSE
-                WRITE(*,*) "wxios_add_file: Error, invalid file: ", "X"//trim(fname)
-                WRITE(*,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
+                WRITE(lunout,*) "wxios_add_file: Error, invalid file: ", "X"//trim(fname)
+                WRITE(lunout,*) "wxios_add_file: output_freq=",TRIM(ADJUSTL(nffreq)),"; output_lvl=",flvl
             END IF
         ELSE
-            IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " dÃ©fined using XML."
-                CALL xios_set_file_attr(fname, enabled=.TRUE.)
+            IF (prt_level >= 10) THEN
+              WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " dÃ©fined using XML."
+            ENDIF
+            ! Ehouarn: add an enable=.true. on top of xml definitions... why???
+            CALL xios_set_file_attr(fname, enabled=.TRUE.)
         END IF
     END SUBROUTINE wxios_add_file
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    ! Pour crÃ©er un champ      !!!!!!!!!!!!!!!
+    ! Pour crÃ©er un champ      !!!!!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUBROUTINE wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit)
@@ -418,5 +422,5 @@
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    ! Pour dÃ©clarer un champ      !!!!!!!!!!!!!!!
+    ! Pour dÃ©clarer un champ      !!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op)
@@ -432,5 +436,5 @@
         CHARACTER(len=*), INTENT(IN) :: op
         
-        CHARACTER(len=20) :: axis_id
+        CHARACTER(len=20) :: axis_id ! Ehouarn: dangerous...
         CHARACTER(len=100) :: operation
         TYPE(xios_file) :: f
@@ -441,6 +445,10 @@
         
         
-        !PrÃ©paration du nom de l'axe:
-        CALL concat("presnivs", fname, axis_id)
+        ! Ajout Abd pour NMC:
+        IF (fid.LE.6) THEN
+          axis_id="presnivs"
+        ELSE
+          axis_id="plev"
+        ENDIF
         
         !on prÃ©pare le nom de l'opÃ©ration:
@@ -448,8 +456,7 @@
         
         
-        
         !On selectionne le bon groupe de champs:
         IF (fdim.EQ.2) THEN
-            CALL xios_get_fieldgroup_handle("fields_2D", fieldgroup)
+          CALL xios_get_fieldgroup_handle("fields_2D", fieldgroup)
         ELSE
           CALL xios_get_fieldgroup_handle("fields_3D", fieldgroup)
@@ -515,26 +522,29 @@
             !Sinon on se contente de l'activer:
             CALL xios_set_field_attr(fieldname, enabled=.TRUE.)
+            !NB: This will override an enable=.false. set by a user in the xml file;
+            !   then the only way to not output the field is by changing its
+            !   output level
         ENDIF       
         
     END SUBROUTINE wxios_add_field_to_file
     
-    SUBROUTINE wxios_update_calendar(ito)
-        INTEGER, INTENT(IN) :: ito
-        CALL xios_update_calendar(ito)
-    END SUBROUTINE wxios_update_calendar
-    
-    SUBROUTINE wxios_write_2D(fieldname, fdata)
-        CHARACTER(len=*), INTENT(IN) :: fieldname
-        REAL, DIMENSION(:,:), INTENT(IN) :: fdata
-        
-        CALL xios_send_field(fieldname, fdata)
-    END SUBROUTINE wxios_write_2D
-    
-    SUBROUTINE wxios_write_3D(fieldname, fdata)
-        CHARACTER(len=*), INTENT(IN) :: fieldname
-        REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata
-        
-        CALL xios_send_field(fieldname, fdata)
-    END SUBROUTINE wxios_write_3D
+!    SUBROUTINE wxios_update_calendar(ito)
+!        INTEGER, INTENT(IN) :: ito
+!        CALL xios_update_calendar(ito)
+!    END SUBROUTINE wxios_update_calendar
+!    
+!    SUBROUTINE wxios_write_2D(fieldname, fdata)
+!        CHARACTER(len=*), INTENT(IN) :: fieldname
+!        REAL, DIMENSION(:,:), INTENT(IN) :: fdata
+!
+!        CALL xios_send_field(fieldname, fdata)
+!    END SUBROUTINE wxios_write_2D
+    
+!    SUBROUTINE wxios_write_3D(fieldname, fdata)
+!        CHARACTER(len=*), INTENT(IN) :: fieldname
+!        REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata
+!        
+!        CALL xios_send_field(fieldname, fdata)
+!    END SUBROUTINE wxios_write_3D
     
     SUBROUTINE wxios_closedef()
Index: /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1302)
@@ -172,9 +172,4 @@
       REAL unskap, pksurcp
       save unskap
-
-cIM diagnostique PVteta, Amip2
-      INTEGER,PARAMETER :: ntetaSTD=3
-      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
-      REAL PVteta(ngridmx,ntetaSTD)
 
       REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
@@ -651,14 +646,4 @@
       ENDDO
 c
-      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,
-     $           ztfi,zplay,zplev,
-     $           ntetaSTD,rtetaSTD,PVteta)
-#endif
-      endif
-c
 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
       CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
@@ -713,7 +698,5 @@
      .             zdqfi,
      .             zdpsrf,
-cIM diagnostique PVteta, Amip2          
-     .             pducov,
-     .             PVteta)
+     .             pducov)
 
       else if ( planet_type=="generic" ) then
Index: /trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90	(revision 1302)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
+! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -30,5 +30,5 @@
 #ifndef CPP_EARTH
 #include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific program, needs Earth physics'
+  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
 #else
 !-------------------------------------------------------------------------------
@@ -94,12 +94,12 @@
   END IF
 
-  IF (grilles_gcm_netcdf) THEN
-     WRITE(lunout,'(//)')
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,'(//)')
-     CALL grilles_gcm_netcdf_sub(masque,phis)
-  END IF
+  
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
+
 #endif
 ! of #ifndef CPP_EARTH #else
@@ -108,2 +108,3 @@
 !
 !-------------------------------------------------------------------------------
+
Index: /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1302)
@@ -2,8 +2,8 @@
 ! $Id: conf_gcm.F 1418 2010-07-19 15:11:24Z jghattas $
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit )
-c
+!
       USE control_mod
 #ifdef CPP_IOIPSL
@@ -18,18 +18,18 @@
 
       IMPLICIT NONE
-c-----------------------------------------------------------------------
-c     Auteurs :   L. Fairhead , P. Le Van  .
-c
-c     Arguments :
-c
-c     tapedef   :
-c     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
-c     -metres  du zoom  avec  celles lues sur le fichier start .
-c
+!-----------------------------------------------------------------------
+!     Auteurs :   L. Fairhead , P. Le Van  .
+!
+!     Arguments :
+!
+!     tapedef   :
+!     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
+!     -metres  du zoom  avec  celles lues sur le fichier start .
+!
        LOGICAL etatinit
        INTEGER tapedef
 
-c   Declarations :
-c   --------------
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -43,8 +43,8 @@
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -54,32 +54,32 @@
       INTEGER i
       LOGICAL use_filtre_fft
-c
-c  -------------------------------------------------------------------
-c
-c       .........     Version  du 29/04/97       ..........
-c
-c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
-c      tetatemp   ajoutes  pour la dissipation   .
-c
-c   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
-c
-c  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
-c    Sinon , choix de fxynew  , a derivee sinusoidale  ..
-c
-c   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
-c         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
-c                de limit.dat ( dic)                        ...........
-c           Sinon  etatinit = . FALSE .
-c
-c   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
-c    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
-c   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
-c    lectba .  
-c   Ces parmetres definissant entre autres la grille et doivent etre
-c   pareils et coherents , sinon il y aura  divergence du gcm .
-c
-c-----------------------------------------------------------------------
-c   initialisations:
-c   ----------------
+!
+!  -------------------------------------------------------------------
+!
+!       .........     Version  du 29/04/97       ..........
+!
+!   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
+!      tetatemp   ajoutes  pour la dissipation   .
+!
+!   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
+!
+!  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
+!    Sinon , choix de fxynew  , a derivee sinusoidale  ..
+!
+!   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
+!         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
+!                de limit.dat ( dic)                        ...........
+!           Sinon  etatinit = . FALSE .
+!
+!   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
+!    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
+!   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
+!    lectba .  
+!   Ces parmetres definissant entre autres la grille et doivent etre
+!   pareils et coherents , sinon il y aura  divergence du gcm .
+!
+!-----------------------------------------------------------------------
+!   initialisations:
+!   ----------------
 
 !Config  Key  = lunout
@@ -91,5 +91,5 @@
       CALL getin('lunout', lunout)
       IF (lunout /= 5 .and. lunout /= 6) THEN
-        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write', 
+        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write',                     &
      &          STATUS='unknown',FORM='formatted')
       ENDIF
@@ -103,7 +103,7 @@
       CALL getin('prt_level',prt_level)
 
-c-----------------------------------------------------------------------
-c  Parametres de controle du run:
-c-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+!  Parametres de controle du run:
+!-----------------------------------------------------------------------
 !Config  Key  = planet_type
 !Config  Desc = planet type ("earth", "mars", "venus", ...)
@@ -264,6 +264,6 @@
        CALL getin('dissip_period',dissip_period)
 
-ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
-ccc
+!cc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
+!cc
 
 !Config  Key  = lstardis
@@ -430,5 +430,5 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!    ...............................................................
 
 !Config  Key  =  read_start
@@ -587,12 +587,6 @@
       CALL getin('ok_etat0',ok_etat0)
 
-!Config  Key  = grilles_gcm_netcdf
-!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
-!Config  Def  = n
-      grilles_gcm_netcdf = .FALSE.
-      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
-
-c----------------------------------------
-c Parameters for zonal averages in the case of Titan
+!----------------------------------------
+! Parameters for zonal averages in the case of Titan
       moyzon_mu = .false.
       moyzon_ch = .false.
@@ -601,13 +595,13 @@
        CALL getin('moyzon_ch', moyzon_ch)
       endif
-c----------------------------------------
-
-c----------------------------------------
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
-C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
-c
-c----------------------------------------
+!----------------------------------------
+
+!----------------------------------------
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
+! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
+!
+!----------------------------------------
       IF( etatinit ) then
 
@@ -645,6 +639,6 @@
 
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
-     &   'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
+        write(lunout,*)                                                        &
+     &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -654,6 +648,6 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
-     &  'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
+        write(lunout,*)                                                        &
+     &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -662,8 +656,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -737,5 +731,5 @@
 c
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
+        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',     &
      &    ' est differente de celle lue sur le fichier  start '
         STOP
@@ -752,5 +746,5 @@
 
       IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de grossismx passee par ',
+        write(lunout,*)'conf_gcm: La valeur de grossismx passee par ',       &
      &  'run.def est differente de celle lue sur le fichier  start '
         STOP
@@ -766,5 +760,5 @@
 
       IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de grossismy passee par ',
+        write(lunout,*)'conf_gcm: La valeur de grossismy passee par ',        &
      & 'run.def est differente de celle lue sur le fichier  start '
         STOP
@@ -772,5 +766,5 @@
       
       IF( grossismx.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
@@ -781,5 +775,5 @@
 
       IF( grossismy.LT.1. )  THEN
-        write(lunout,*)
+        write(lunout,*)                                                        &
      &       'conf_gcm: ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
@@ -789,8 +783,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -805,6 +799,6 @@
          IF( fxyhypbb )     THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** fxyhypb lu sur le fichier start est ',
-     *       'F alors  qu il est  T  sur  run.def  ***'
+            write(lunout,*)' *** fxyhypb lu sur le fichier start est ',     &
+     &       'F alors  qu il est  T  sur  run.def  ***'
               STOP
          ENDIF
@@ -812,10 +806,10 @@
          IF( .NOT.fxyhypbb )   THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' ***  fxyhypb lu sur le fichier start est ',
-     *        'T alors  qu il est  F  sur  run.def  ****  '
+            write(lunout,*)' ***  fxyhypb lu sur le fichier start est ',     &
+     &        'T alors  qu il est  F  sur  run.def  ****  '
               STOP
          ENDIF
       ENDIF
-c
+!
 !Config  Key  = dzoomx
 !Config  Desc = extension en longitude
@@ -828,6 +822,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de dzoomx passee par ',
-     *  'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomx passee par ',         &
+     &  'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -844,6 +838,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de dzoomy passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de dzoomy passee par ',          &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -859,6 +853,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(taux - tauxx).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de taux passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de taux passee par ',           &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -874,6 +868,6 @@
       IF( fxyhypb )  THEN
        IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
-        write(lunout,*)'conf_gcm: La valeur de tauy passee par ',
-     * 'run.def est differente de celle lue sur le fichier  start '
+        write(lunout,*)'conf_gcm: La valeur de tauy passee par ',           &
+     & 'run.def est differente de celle lue sur le fichier  start '
         STOP
        ENDIF
@@ -895,6 +889,6 @@
           IF( ysinuss )     THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** ysinus lu sur le fichier start est F',
-     *       ' alors  qu il est  T  sur  run.def  ***'
+            write(lunout,*)' *** ysinus lu sur le fichier start est F',     &
+     &       ' alors  qu il est  T  sur  run.def  ***'
             STOP
           ENDIF
@@ -902,6 +896,6 @@
           IF( .NOT.ysinuss )   THEN
             write(lunout,*)' ********  PBS DANS  CONF_GCM  ******** '
-            write(lunout,*)' *** ysinus lu sur le fichier start est T',
-     *        ' alors  qu il est  F  sur  run.def  ****  '
+            write(lunout,*)' *** ysinus lu sur le fichier start est T',     &
+     &        ' alors  qu il est  F  sur  run.def  ****  '
               STOP
           ENDIF
@@ -910,5 +904,5 @@
 
       endif ! etatinit
-c----------------------------------------
+!----------------------------------------
 
 
@@ -962,5 +956,4 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
       if (planet_type=="titan") then
        write(lunout,*)' moyzon_mu = ', moyzon_mu
Index: /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 1302)
@@ -107,7 +107,4 @@
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL masse(ip1jmp1,llm)                ! masse d'air
       REAL phis(ip1jmp1)                     ! geopotentiel au sol
@@ -133,5 +130,4 @@
       data call_iniphys/.true./
 
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 c+jld variables test conservation energie
 c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
@@ -500,5 +496,9 @@
 c   ------------------------
 
-      day_end=day_ini+nday
+      if (nday>=0) then ! standard case
+        day_end=day_ini+nday
+      else ! special case when nday <0, run -nday dynamical steps
+        day_end=day_ini-nday/day_step
+      endif
       if (less1day) then
         day_end=day_ini+floor(time_0+fractday)
Index: /trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 1302)
@@ -437,5 +437,5 @@
 ! Sauvegarde du guidage?
     f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)  
-    IF (f_out) CALL guide_out("S",jjp1,1,ps)
+    IF (f_out) CALL guide_out("SP",jjp1,1,ps)
     
     if (guide_u) then
@@ -447,5 +447,7 @@
         if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
-        IF (f_out) CALL guide_out("U",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
+        IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
+        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
         ucov=ucov+f_add
     endif
@@ -459,5 +461,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
-        IF (f_out) CALL guide_out("T",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
         teta=teta+f_add
     endif
@@ -471,5 +473,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
-        IF (f_out) CALL guide_out("P",jjp1,1,f_add(1:ip1jmp1,1)/factt)
+        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
         ps=ps+f_add(1:ip1jmp1,1)
         CALL pression(ip1jmp1,ap,bp,ps,p)
@@ -485,5 +487,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
-        IF (f_out) CALL guide_out("Q",jjp1,llm,f_add/factt)
+        IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
         q=q+f_add
     endif
@@ -497,5 +499,7 @@
         if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
-        IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:)/factt)
+        IF (f_out) CALL guide_out("v",jjm,llm,vcov)
+        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
+        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
         vcov=vcov+f_add(1:ip1jm,:)
     endif
@@ -589,4 +593,6 @@
   SUBROUTINE guide_interp(psi,teta)
   
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
   IMPLICIT NONE
 
@@ -610,6 +616,5 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-  REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: prefkap,unskap
@@ -676,7 +681,7 @@
     CALL pression( ip1jmp1, ap, bp, psi, p )
     if (pressure_exner) then
-      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
     else
-      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
+      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
     endif
 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
@@ -1507,5 +1512,5 @@
     
     ! Variables entree
-    CHARACTER, INTENT(IN)                          :: varname
+    CHARACTER*(*), INTENT(IN)                          :: varname
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
@@ -1516,7 +1521,9 @@
     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
+    INTEGER       :: vid_au,vid_av
     INTEGER, DIMENSION (3) :: dim3
     INTEGER, DIMENSION (4) :: dim4,count,start
-    INTEGER                :: ierr, varid
+    INTEGER                :: ierr, varid,l
+    REAL, DIMENSION (iip1,hsize,vsize) :: field2
 
     print *,'Guide: output timestep',timestep,'var ',varname
@@ -1542,5 +1549,7 @@
         ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
         ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
+        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
         ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
         
         ierr=NF_ENDDEF(nid)
@@ -1555,4 +1564,6 @@
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
 #else
         ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
@@ -1563,4 +1574,6 @@
         ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
 #endif
 ! --------------------------------------------------------------------
@@ -1579,4 +1592,6 @@
         IF (guide_u) THEN
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1584,4 +1599,6 @@
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1606,62 +1623,36 @@
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
+    IF (varname=="SP") timestep=timestep+1
+
+    ierr = NF_INQ_VARID(nid,varname,varid)
     SELECT CASE (varname)
-    CASE ("S")
-        timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"SP",varid)
+    CASE ("SP","ps")
         start=(/1,1,timestep,0/)
         count=(/iip1,jjp1,1,0/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("P")
-        ierr = NF_INQ_VARID(nid,"ps",varid)
-        start=(/1,1,timestep,0/)
-        count=(/iip1,jjp1,1,0/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("U")
-        ierr = NF_INQ_VARID(nid,"ucov",varid)
+    CASE ("v","va","vcov")
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjm,llm,1/)
+    CASE DEFAULT
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
+    END SELECT
+
+    SELECT CASE (varname)
+    CASE("u","ua")
+        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
+        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
+    CASE("v","va")
+        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE DEFAULT
+        field2=field
+    END SELECT
+
+
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
+    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
 #else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
+    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
 #endif
-    CASE ("V")
-        ierr = NF_INQ_VARID(nid,"vcov",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjm,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("T")
-        ierr = NF_INQ_VARID(nid,"teta",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("Q")
-        ierr = NF_INQ_VARID(nid,"q",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    END SELECT
- 
+
     ierr = NF_CLOSE(nid)
 
Index: /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1302)
@@ -20,4 +20,6 @@
      &                       ok_dynzon,periodav,ok_dyn_ave,iecri,
      &                       ok_dyn_ins,output_grads_dyn
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
       use cpdet_mod, only: cpdet,tpot2t,t2tpot
       use sponge_mod, only: callsponge,mode_sponge,sponge
@@ -217,5 +219,10 @@
       endif
 
-      itaufin   = nday*day_step
+      if (nday>=0) then
+         itaufin   = nday*day_step
+      else
+         ! to run a given (-nday) number of dynamical steps
+         itaufin   = -nday
+      endif
       if (less1day) then
 c MODIF VENUS: to run less than one day:
@@ -262,7 +269,7 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 
@@ -476,8 +483,11 @@
          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
          if (pressure_exner) then
-           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
          endif
+
+! Compute geopotential (physics might need it)
+         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 
            jD_cur = jD_ref + day_ini - day_ref +                        &
@@ -551,7 +561,7 @@
           CALL massdair(p,masse)
           if (pressure_exner) then
-            CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
           else
-            CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
           endif
           
@@ -604,7 +614,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
         endif
         CALL massdair(p,masse)
Index: /trunk/LMDZ.COMMON/libf/dyn3d/logic.h
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/logic.h	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/logic.h	(revision 1302)
@@ -11,5 +11,5 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
+     &  ,ok_limit,ok_etat0,hybrid                                       &
      &  ,moyzon_mu,moyzon_ch
 
@@ -19,5 +19,6 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0
+
       logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                      ! (only used if disvert_type==2)
Index: unk/LMDZ.COMMON/libf/dyn3d_common/PVtheta.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/PVtheta.F	(revision 1301)
+++ 	(revision )
@@ -1,196 +1,0 @@
-      SUBROUTINE PVtheta(ilon,ilev,pucov,pvcov,pteta,
-     $           ztfi,zplay,zplev,
-     $           nbteta,theta,PVteta)
-      IMPLICIT none
-
-c=======================================================================
-c
-c   Auteur:  I. Musat
-c   -------
-c
-c   Objet:
-c   ------
-c
-c    *******************************************************************
-c    Calcul de la vorticite potentielle PVteta sur des iso-theta selon
-c    la methodologie du NCEP/NCAR :
-c    1) on calcule la stabilite statique N**2=g/T*(dT/dz+g/cp) sur les
-c       niveaux du modele => N2
-c    2) on interpole les vents, la temperature et le N**2 sur des isentropes
-c       (en fait sur des iso-theta) lineairement en log(theta) =>
-c       ucovteta, vcovteta, N2teta
-c    3) on calcule la vorticite absolue sur des iso-theta => vorateta
-c    4) on calcule la densite rho sur des iso-theta => rhoteta 
-c
-c       rhoteta = (T/theta)**(cp/R)*p0/(R*T)
-c
-c    5) on calcule la vorticite potentielle sur des iso-theta => PVteta
-c
-c       PVteta = (vorateta * N2 * theta)/(g * rhoteta) ! en PVU
-c
-c       NB: 1PVU=10**(-6) K*m**2/(s * kg)
-c
-c       PVteta =  vorateta * N2/(g**2 * rhoteta) ! en 1/(Pa*s)
-c
-c
-c    *******************************************************************
-c
-c
-c     Variables d'entree : ilon,ilev,pucov,pvcov,pteta,ztfi,zplay,zplev,nbteta,theta
-c                       -> sur la grille dynamique
-c     Variable de sortie : PVteta
-c                       -> sur la grille physique 
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-c
-c variables Input
-c
-      INTEGER ilon, ilev
-      REAL pvcov(iip1,jjm,ilev)
-      REAL pucov(iip1,jjp1,ilev)
-      REAL pteta(iip1,jjp1,ilev)
-      REAL ztfi(ilon,ilev)
-      REAL zplay(ilon,ilev), zplev(ilon,ilev+1)
-      INTEGER nbteta
-      REAL theta(nbteta)
-c
-c variable Output
-c
-      REAL PVteta(ilon,nbteta)
-c
-c variables locales
-c
-      INTEGER i, j, l, ig0
-      REAL SSUM
-      REAL teta(ilon, ilev)
-      REAL ptetau(ip1jmp1, ilev), ptetav(ip1jm, ilev)
-      REAL ucovteta(ip1jmp1,ilev), vcovteta(ip1jm,ilev)
-      REAL N2(ilon,ilev-1), N2teta(ilon,nbteta)
-      REAL ztfiteta(ilon,nbteta)
-      REAL rhoteta(ilon,nbteta)
-      REAL vorateta(iip1,jjm,nbteta)
-      REAL voratetafi(ilon,nbteta), vorpol(iim)
-c
-#include "comgeom2.h"
-#include "comconst.h"
-#include "comvert.h"
-c
-c projection teta sur la grille physique
-c
-      DO l=1,llm
-       teta(1,l)   =  pteta(1,1,l)
-       ig0         = 2
-       DO j = 2, jjm
-        DO i = 1, iim
-         teta(ig0,l)    = pteta(i,j,l)
-         ig0            = ig0 + 1
-        ENDDO
-       ENDDO
-       teta(ig0,l)    = pteta(1,jjp1,l)
-      ENDDO
-c
-c calcul pteta sur les grilles U et V
-c
-      DO l=1, llm
-       DO j=1, jjp1
-        DO i=1, iip1
-         ig0=i+(j-1)*iip1
-         ptetau(ig0,l)=pteta(i,j,l)
-        ENDDO !i
-       ENDDO !j
-       DO j=1, jjm
-        DO i=1, iip1
-         ig0=i+(j-1)*iip1
-         ptetav(ig0,l)=0.5*(pteta(i,j,l)+pteta(i,j+1,l))
-        ENDDO !i
-       ENDDO !j
-      ENDDO !l
-c
-c projection pucov, pvcov sur une surface de theta constante
-c
-      DO l=1, nbteta
-cIM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l),
-       CALL tetalevel(ip1jmp1,llm,.true.,ptetau,theta(l),
-     .                pucov,ucovteta(:,l))
-cIM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l),
-       CALL tetalevel(ip1jm,llm,.true.,ptetav,theta(l),
-     .                pvcov,vcovteta(:,l))
-      ENDDO !l
-c
-c calcul vorticite absolue sur une iso-theta : vorateta
-c
-      CALL tourabs(nbteta,vcovteta,ucovteta,vorateta)
-c
-c projection vorateta sur la grille physique => voratetafi
-c
-      DO l=1,nbteta
-       DO j=2,jjm
-        ig0=1+(j-2)*iim
-        DO i=1,iim
-         voratetafi(ig0+i+1,l) = vorateta( i ,j-1,l) * alpha4(i+1,j) +
-     $                           vorateta(i+1,j-1,l) * alpha1(i+1,j) +
-     $                           vorateta(i  ,j  ,l) * alpha3(i+1,j) +
-     $                           vorateta(i+1,j  ,l) * alpha2(i+1,j)
-        ENDDO
-        voratetafi(ig0 +1,l) = voratetafi(ig0 +1+ iim,l)
-       ENDDO
-      ENDDO
-c
-      DO l=1,nbteta
-       DO i=1,iim
-        vorpol(i)  = vorateta(i,1,l)*aire(i,1)
-       ENDDO
-       voratetafi(1,l)= SSUM(iim,vorpol,1)/apoln
-      ENDDO
-c
-      DO l=1,nbteta
-       DO i=1,iim
-        vorpol(i)  = vorateta(i,jjm,l)* aire(i,jjm +1)
-       ENDDO
-       voratetafi(ilon,l)= SSUM(iim,vorpol,1)/apols
-      ENDDO
-c 
-c calcul N**2 sur la grille physique => N2
-c
-      DO l=1, llm-1 
-       DO i=1, ilon
-        N2(i,l) = (g**2 * zplay(i,l) * 
-     $            (ztfi(i,l+1)-ztfi(i,l)) )/
-     $            (R*ztfi(i,l)*ztfi(i,l)*
-     $            (zplev(i,l)-zplev(i,l+1)) )+
-     $            (g**2)/(ztfi(i,l)*CPP)
-       ENDDO !i
-      ENDDO !l
-c
-c calcul N2 sur une iso-theta => N2teta 
-c
-      DO l=1, nbteta
-       CALL tetalevel(ilon,llm-1,.true.,teta,theta(l),
-     $                N2,N2teta(:,l))
-       CALL tetalevel(ilon,llm,.true.,teta,theta(l),
-     $                ztfi,ztfiteta(:,l))
-      ENDDO !l=1, nbteta
-c
-c calcul rho et PV sur une iso-theta : rhoteta, PVteta
-c
-      DO l=1, nbteta
-       DO i=1, ilon
-        rhoteta(i,l)=(ztfiteta(i,l)/theta(l))**(CPP/R)*
-     $  (preff/(R*ztfiteta(i,l)))
-c
-c PVteta en PVU
-c
-        PVteta(i,l)=(theta(l)*g*voratetafi(i,l)*N2teta(i,l))/
-     $              (g**2*rhoteta(i,l))
-c
-c PVteta en 1/(Pa*s)
-c
-        PVteta(i,l)=(voratetafi(i,l)*N2teta(i,l))/
-     $              (g**2*rhoteta(i,l))
-       ENDDO !i
-      ENDDO !l
-c
-      RETURN
-      END 
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/disvert.F90	(revision 1302)
@@ -1,10 +1,14 @@
 ! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
 
-SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
+SUBROUTINE disvert()
 
   ! Auteur : P. Le Van
 
+#ifdef CPP_IOIPSL
+  use ioipsl, only: getin
+#else
+  USE ioipsl_getincom, only: getin
+#endif
   use new_unit_m, only: new_unit
-  use ioipsl, only: getin
   use assert_m, only: assert
 
@@ -13,34 +17,44 @@
   include "dimensions.h"
   include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
   include "iniprint.h"
   include "logic.h"
 
-  ! s = sigma ** kappa : coordonnee verticale
-  ! dsig(l) : epaisseur de la couche l ds la coord. s
-  ! sig(l) : sigma a l'interface des couches l et l-1
-  ! ds(l) : distance entre les couches l et l-1 en coord.s
-
-  real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1) ! in Pa
-  real, intent(out):: bp(llmp1)
-  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
-  real,intent(out) :: presnivs(llm)
-  real,intent(out) :: scaleheight
-
+!-------------------------------------------------------------------------------
+! Purpose: Vertical distribution functions for LMDZ.
+!          Triggered by the levels number llm.
+!-------------------------------------------------------------------------------
+! Read    in "comvert.h":
+! pa                         !--- PURE PRESSURE COORDINATE FOR P<pa (in Pascals)
+! preff                      !--- REFERENCE PRESSURE                 (101325 Pa)
+! Written in "comvert.h":
+! ap(llm+1), bp(llm+1)       !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
+! aps(llm),  bps(llm)        !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
+! dpres(llm)                 !--- PRESSURE DIFFERENCE FOR EACH LAYER
+! presnivs(llm)              !--- PRESSURE AT EACH MID-LAYER
+! scaleheight                !--- VERTICAL SCALE HEIGHT            (Earth: 8kms)
+! nivsig(llm+1)              !--- SIGMA INDEX OF EACH LAYER INTERFACE
+! nivsigs(llm)               !--- SIGMA INDEX OF EACH MID-LAYER
+!-------------------------------------------------------------------------------
+! Local variables:
   REAL sig(llm+1), dsig(llm)
-  real zk, zkm1, dzk1, dzk2, k0, k1
+  REAL sig0(llm+1), zz(llm+1)
+  REAL zk, zkm1, dzk1, dzk2, z, k0, k1
 
   INTEGER l, unit
   REAL dsigmin
+  REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
+
   REAL alpha, beta, deltaz
   REAL x
   character(len=*),parameter :: modname="disvert"
 
-  character(len=6):: vert_sampling
+  character(len=24):: vert_sampling
   ! (allowed values are "param", "tropo", "strato" and "read")
 
   !-----------------------------------------------------------------------
 
-  print *, "Call sequence information: disvert"
+  WRITE(lunout,*) TRIM(modname)//" starts"
 
   ! default scaleheight is 8km for earth
@@ -49,5 +63,13 @@
   vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
   call getin('vert_sampling', vert_sampling)
-  print *, 'vert_sampling = ' // vert_sampling
+  WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
+  if (llm==39 .and. vert_sampling=="strato") then
+     dsigmin=0.3 ! Vieille option par dÃ©faut pour CMIP5
+  else
+     dsigmin=1.
+  endif
+  call getin('dsigmin', dsigmin)
+  WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
+
 
   select case (vert_sampling)
@@ -86,8 +108,8 @@
 
      ap = pa * (sig - bp)
-  case("tropo")
+  case("sigma")
      DO l = 1, llm
         x = 2*asin(1.) * (l - 0.5) / (llm + 1)
-        dsig(l) = 1.0 + 7.0 * SIN(x)**2
+        dsig(l) = dsigmin + 7.0 * SIN(x)**2
      ENDDO
      dsig = dsig / sum(dsig)
@@ -98,4 +120,19 @@
 
      bp(1)=1.
+     bp(2: llm) = sig(2:llm)
+     bp(llmp1) = 0.
+     ap(:)=0.
+  case("tropo")
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+        dsig(l) = dsigmin + 7.0 * SIN(x)**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.
@@ -104,14 +141,4 @@
      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)
@@ -131,4 +158,154 @@
      ap(1)=0.
      ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
+  case("strato_correct")
+!==================================================================
+! Fredho 2014/05/18, Saint-Louis du Senegal
+! Cette version de la discretisation strato est corrige au niveau
+! du passage des sig aux ap, bp
+! la version precedente donne un coude dans l'epaisseur des couches
+! vers la tropopause
+!==================================================================
+
+
+     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)
+     sig0(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig0(l) = sig0(l+1) + dsig(l)
+     ENDDO
+     sig=racinesig(sig0)
+
+     bp(1)=1.
+     bp(2:llm)=EXP(1.-1./sig(2: llm)**2)
+     bp(llmp1)=0.
+
+     ap(1)=0.
+     ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
+     ap(llm+1)=0.
+
+  CASE("strato_custom0") 
+!=======================================================
+! Version Transitoire
+    ! custumize strato distribution with specific alpha & beta values and function
+    ! depending on llm (experimental and temporary)!
+    SELECT CASE (llm)
+      CASE(55)
+        alpha=0.45
+        beta=4.0
+      CASE(63)
+        alpha=0.45
+        beta=5.0
+      CASE(71)
+        alpha=3.05
+        beta=65.
+      CASE(79)
+        alpha=3.20
+        ! alpha=2.05 ! FLOTT 79 (PLANTE)
+        beta=70.
+    END SELECT
+    ! Or used values provided by user in def file:
+    CALL getin("strato_alpha",alpha)
+    CALL getin("strato_beta",beta)
+    
+    ! Build geometrical distribution
+    scaleheight=7.
+    zz(1)=0.
+    IF (llm==55.OR.llm==63) THEN
+      DO l=1,llm
+        z=zz(l)/scaleheight
+        zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5))     &
+                            +5.0*EXP((l-llm)/beta)
+      ENDDO
+    ELSEIF (llm==71) THEN !.OR.llm==79) THEN      ! FLOTT 79 (PLANTE)
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
+      ENDDO
+    ELSEIF (llm==79) THEN
+      DO l=1,llm
+        z=zz(l)
+        zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.))     &
+                            +0.03*TANH(z/.25)
+      ENDDO
+    ENDIF ! of IF (llm==55.OR.llm==63) ...
+    
+
+    ! Build sigma distribution
+    sig0=EXP(-zz(:)/scaleheight)
+    sig0(llm+1)=0.
+!    sig=ridders(sig0)
+    sig=racinesig(sig0)
+    
+    ! Compute ap() and bp()
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
+  CASE("strato_custom") 
+!===================================================================
+! David Cugnet, FranÃ§ois Lott, Lionel Guez, Ehouoarn Millour, Fredho
+! 2014/05
+! custumize strato distribution
+! Al the parameter are given in km assuming a given scalehigh
+    vert_scale=7.     ! scale hight
+    vert_dzmin=0.02   ! width of first layer
+    vert_dzlow=1.     ! dz in the low atmosphere
+    vert_z0low=8.     ! height at which resolution recches dzlow
+    vert_dzmid=3.     ! dz in the mid atmsophere
+    vert_z0mid=70.    ! height at which resolution recches dzmid
+    vert_h_mid=20.    ! width of the transition
+    vert_dzhig=11.    ! dz in the high atmsophere
+    vert_z0hig=80.    ! height at which resolution recches dz
+    vert_h_hig=20.    ! width of the transition
+!===================================================================
+
+    call getin('vert_scale',vert_scale)
+    call getin('vert_dzmin',vert_dzmin)
+    call getin('vert_dzlow',vert_dzlow)
+    call getin('vert_z0low',vert_z0low)
+    CALL getin('vert_dzmid',vert_dzmid)
+    CALL getin('vert_z0mid',vert_z0mid)
+    call getin('vert_h_mid',vert_h_mid)
+    call getin('vert_dzhig',vert_dzhig)
+    call getin('vert_z0hig',vert_z0hig)
+    call getin('vert_h_hig',vert_h_hig)
+
+    scaleheight=vert_scale ! for consistency with further computations
+    ! Build geometrical distribution
+    zz(1)=0.
+    DO l=1,llm
+       z=zz(l)
+       zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+                &
+&      (vert_dzmid-vert_dzlow)* &
+&           (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
+&     +(vert_dzhig-vert_dzmid-vert_dzlow)*                                  &
+&           (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
+    ENDDO
+
+
+!===================================================================
+! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
+! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
+! sig0 is p/p0
+! sig is an intermediate distribution introduce to estimate bp
+! 1.   sig0=exp(-z/H)
+! 2.   inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
+! 3.   bp=exp(1-1/sig**2)
+! 4.   ap deduced from  the combination of 2 and 3 so that sig0=ap/p0+bp
+!===================================================================
+
+    sig0=EXP(-zz(:)/vert_scale)
+    sig0(llm+1)=0.
+    sig=racinesig(sig0)
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
   case("read")
      ! Read "ap" and "bp". First line is skipped (title line). "ap"
@@ -178,3 +355,101 @@
   write(lunout, *) presnivs
 
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION ridders(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
+! Notes:   Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
+! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
+!       Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=9999
+  c1=Pa/Preff; c2=1.-c1
+  DO l=1,ns
+    xx=HUGE(1.)
+    x1=0.0; f1=distrib(x1,c1,c2,sig(l))
+    x2=1.0; f2=distrib(x2,c1,c2,sig(l))
+    DO it=1,maxit
+      x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
+      s=SQRT(f3**2-f1*f2);                 IF(s==0.) EXIT
+      x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
+      xx=x4; f4=distrib(x4,c1,c2,sig(l));  IF(f4==0.) EXIT
+      IF(SIGN(f3,f4)/=f3) THEN;      x1=x3; f1=f3; x2=xx; f2=f4
+      ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
+      ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
+      ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
+      END IF
+      IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT       !--- ERROR ON SIG <= 0.01m            
+    END DO
+    IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
+     &e for level ',l
+    sg(l)=xx
+  END DO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION ridders
+
+FUNCTION racinesig(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Fredho 2014/05/18
+! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
+! Notes:   Uses Newton Raphson search
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=100
+  c1=Pa/Preff; c2=1.-c1
+  DO l=2,ns-1
+    sg(l)=sig(l)
+    DO it=1,maxit
+       f1=exp(1-1./sg(l)**2)*(1.-c1)
+       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
+    ENDDO
+!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
+  ENDDO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION racinesig
+
+
+
+
 END SUBROUTINE disvert
+
+!-------------------------------------------------------------------------------
+
+FUNCTION distrib(x,c1,c2,x0) RESULT(res)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN) :: x, c1, c2, x0
+  REAL             :: res
+!-------------------------------------------------------------------------------
+  res=c1*x+c2*EXP(1-1/(x**2))-x0
+
+END FUNCTION distrib
+
+
+
Index: unk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb.F	(revision 1301)
+++ 	(revision )
@@ -1,152 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
-c
-c     Auteurs :  P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c                                 -------- z                                   
-c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
-c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
-c    ( voir note de Fr.Hourdin )  ,
-c
-c    on determine successivement , du haut vers le bas des couches, les 
-c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
-c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
-c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
-c
-c
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL unpl2k,dellta
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-c
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_hyb"
-      
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-      if (llm.eq.1) then
-        
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
-          pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
-
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-     
-      unpl2k    = 1.+ 2.* kappa
-c
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
-
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-      DO     ij      = 1, ngrid
-       alpha(ij,llm) = 0.
-       beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-        DO ij = 1, ngrid
-        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
-        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
-        beta (ij,l)  =   p(ij,l  ) / dellta   
-        ENDDO
-c
-      ENDDO
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-      DO   ij   = 1, ngrid
-       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
-     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
-      ENDDO
-c
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-      ENDDO
-c
-c
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-
-      RETURN
-      END
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90	(revision 1302)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_hyb_m.F90	(revision 1302)
@@ -0,0 +1,146 @@
+module exner_hyb_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
+
+    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_hyb"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+       ENDDO
+
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    DO     ij      = 1, ngrid
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       DO ij = 1, ngrid
+          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
+          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
+          beta (ij,l)  =   p(ij,l  ) / dellta   
+       ENDDO
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    DO   ij   = 1, ngrid
+       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
+            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
+    ENDDO
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       DO   ij   = 1, ngrid
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb
+
+end module exner_hyb_m
+
Index: unk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu.F	(revision 1301)
+++ 	(revision )
@@ -1,144 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL dum1
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      logical,save :: firstcall=.true.
-      character(len=*),parameter :: modname="exner_milieu"
-
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-!!!! Specific behaviour for Shallow Water (1 vertical layer) case:
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        
-        DO   ij  = 1, ngrid
-          pks(ij) = (cpp/preff) * ps(ij) 
-          pk(ij,1) = .5*pks(ij)
-        ENDDO
-        
-        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
-        
-        ! our work is done, exit routine
-        return
-
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      DO   ij  = 1, ngrid
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-
-      DO  ij   = 1, iim
-        ppn(ij) = aire(   ij   ) * pks(  ij     )
-        pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-      ENDDO
-      xpn      = SSUM(iim,ppn,1) /apoln
-      xps      = SSUM(iim,pps,1) /apols
-
-      DO ij   = 1, iip1
-        pks(   ij     )  =  xpn
-        pks( ij+ip1jm )  =  xps
-      ENDDO
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-        DO   ij   = 1, ngrid
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-      DO   ij   = 1, ngrid
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-
-
-c    calcul de pkf
-c    -------------
-      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-        DO   ij   = 1, ngrid
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-      ENDDO
-
-      RETURN
-      END
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu_m.F90	(revision 1302)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/exner_milieu_m.F90	(revision 1302)
@@ -0,0 +1,125 @@
+module exner_milieu_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_milieu ( ngrid, ps, p, pks, pk, pkf )
+    !
+    !     Auteurs :  F. Forget , Y. Wanherdrick
+    ! P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !     WARNING : CECI est une version speciale de exner_hyb originale
+    !               Utilise dans la version martienne pour pouvoir 
+    !               tourner avec des coordonnees verticales complexe
+    !              => Il ne verifie PAS la condition la proportionalite en 
+    !              energie totale/ interne / potentielle (F.Forget 2001)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    real, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    character(len=*),parameter :: modname="exner_milieu"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+
+       DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+       ENDDO
+
+       if (present(pkf)) then
+          pkf = pk
+          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    DO   ij  = 1, ngrid
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       DO   ij   = 1, ngrid
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+    ENDDO
+
+    !    .... Calcul de pk  pour la couche l = llm ..
+    !    (on met la meme distance (en log pression)  entre Pk(llm)
+    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+    DO   ij   = 1, ngrid
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+       pkf = pk
+       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu
+
+end module exner_milieu_m
+
Index: unk/LMDZ.COMMON/libf/dyn3d_common/grilles_gcm_netcdf.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/grilles_gcm_netcdf.F	(revision 1301)
+++ 	(revision )
@@ -1,305 +1,0 @@
-!
-! $Id: grilles_gcm_netcdf.F 1403 2010-07-01 09:02:53Z fairhead $
-!
-c
-c
-
-      PROGRAM create_fausse_var
-C
-      IMPLICIT NONE
-C
-C
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-
-      real temp(iim+1,jjm+1)
-#include "netcdf.inc"
-
-c Attributs netcdf sortie
-        character*64 fich_out
-        integer*4 ncid_out,rcode_out
-        integer*4 out_lonuid,out_lonvid,out_latuid,out_latvid
-        integer*4 out_varid
-        integer*4 out_lonudim,out_lonvdim
-        integer*4 out_latudim,out_latvdim,out_dim(3)
-
-      INTEGER         longcles
-      PARAMETER     ( longcles = 20 )
-      REAL  clesphy0( longcles )
-
-      integer start(4),count(4)
-
-	integer status,i,j
-        real rlatudeg(jjp1),rlatvdeg(jjm)
-        real rlonudeg(iip1),rlonvdeg(iip1)
-
-      real dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
-      real acoslat,dxkm,dykm,resol(iip1,jjp1)
-
-#include "serre.h"
-#include "fxyprim.h"
-
-      print*,'OK0'
-
-      rad = 6400000
-      omeg = 7.272205e-05
-      g = 9.8
-      kappa = 0.285716
-      daysec = 86400
-      cpp = 1004.70885
-
-      preff = 101325.
-      pa= 50000.
-
-c     open(99,file='run.def',status='old',form='formatted')
-c     CALL defrun_new( 99, .TRUE.,clesphy0 )
-c     close(99)
-
-      CALL conf_gcm( 99, .TRUE. , clesphy0 )
-      CALL iniconst
-      CALL inigeom
-
-
-      print*,'OK1'
-      do j=1,jjp1
-         rlatudeg(j)=rlatu(j)*180./pi
-      enddo
-      do j=1,jjm
-         rlatvdeg(j)=rlatv(j)*180./pi
-      enddo
-
-      do i=1,iip1
-         rlonudeg(i)=rlonu(i)*180./pi + 360.
-         rlonvdeg(i)=rlonv(i)*180./pi + 360.
-      enddo
-
-
-      print*,'OK2'
-c  2 ----- OUVERTURE DE LA SORTIE NETCDF
-c ---------------------------------------------------
-c CREATION OUTPUT
-c ouverture fichier netcdf de sortie out
-        fich_out='grilles_gcm.nc'
-
-        status=NF_CREATE(fich_out,NF_NOCLOBBER,ncid_out)
-        status=NF_DEF_DIM(ncid_out,'lonu',iim+1,out_lonudim)
-        status=NF_DEF_DIM(ncid_out,'lonv',iim+1,out_lonvdim)
-        status=NF_DEF_DIM(ncid_out,'latu',jjm+1,out_latudim)
-        status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
-
-
-      print*,'OK3'
-c   Longitudes en u
-        print *,'OUTID: ',ncid_out
-        status=NF_DEF_VAR(ncid_out,'lonu',NF_FLOAT,1,out_lonudim,
-     %  out_lonuid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'units',
-     %  12,'degrees_east')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'long_name',
-     %  9,'Longitude en u')
-
-c   Longitudes en v
-        print *,'OUTID: ',ncid_out
-        status=NF_DEF_VAR(ncid_out,'lonv',NF_FLOAT,1,out_lonvdim,
-     %  out_lonvid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'units',
-     %  12,'degrees_east')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'long_name',
-     %  9,'Longitude en v')
-
-c   Latitude en u
-        status=NF_DEF_VAR(ncid_out,'latu',NF_FLOAT,1,out_latudim,
-     %  out_latuid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'units',
-     %  13,'degrees_north')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'long_name',
-     %  8,'Latitude en u')
-
-c  Latitude en v
-        status=NF_DEF_VAR(ncid_out,'latv',NF_FLOAT,1,out_latvdim,
-     %  out_latvid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'units',
-     %  13,'degrees_north')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'long_name',
-     %  8,'Latitude en v')
-
-c   ecriture de la grille u
-        out_dim(1)=out_lonudim
-        out_dim(2)=out_latudim
-        status=NF_DEF_VAR(ncid_out,'grille_u',NF_FLOAT,2,out_dim,
-     %  out_varid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units',
-     %  6,'Kelvin')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',
-     %  16,'Grille aux point u')
-
-c   ecriture de la grille v
-        out_dim(1)=out_lonvdim
-        out_dim(2)=out_latvdim
-        status=NF_DEF_VAR(ncid_out,'grille_v',NF_FLOAT,2,out_dim,
-     %  out_varid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units',
-     %  6,'Kelvin')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',
-     %  16,'Grille aux point v')
-
-c   ecriture de la grille u
-        out_dim(1)=out_lonvdim
-        out_dim(2)=out_latudim
-        status=NF_DEF_VAR(ncid_out,'grille_s',NF_FLOAT,2,out_dim,
-     %  out_varid)
-        call handle_err(status)
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units',
-     %  6,'Kelvin')
-        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',
-     %  16,'Grille aux point u')
-
-
-      print*,'OK4'
-        status=NF_ENDDEF(ncid_out)
-c 5) ----- FERMETURE DES FICHIERS NETCDF------------------
-c --------------------------------------------------------
-c 3-b- Ecriture de la grille pour la sortie
-c rajoute l'ecriture de la grille
-
-#ifdef NC_DOUBLE
-      status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonuid,1,iim+1,rlonudeg)
-      status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
-      status=NF_PUT_VARA_DOUBLE(ncid_out,out_latuid,1,jjm+1,rlatudeg)
-      status=NF_PUT_VARA_DOUBLE(ncid_out,out_latvid,1,jjm,rlatvdeg)
-#else
-      status=NF_PUT_VARA_REAL(ncid_out,out_lonuid,1,iim+1,rlonudeg)
-      status=NF_PUT_VARA_REAL(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
-      status=NF_PUT_VARA_REAL(ncid_out,out_latuid,1,jjm+1,rlatudeg)
-      status=NF_PUT_VARA_REAL(ncid_out,out_latvid,1,jjm,rlatvdeg)
-#endif
-
-        start(1)=1
-        start(2)=1
-        start(3)=1
-        start(4)=1
-
-        count(1)=iim+1
-        count(2)=jjm+1
-        count(3)=1
-        count(4)=1
-
-        do j=1,jjm+1
-           do i=1,iim+1
-              temp(i,j)=mod(i,2)+mod(j,2)
-           enddo
-        enddo
-
-#ifdef NC_DOUBLE
-        status=NF_PUT_VARA_DOUBLE(ncid_out,out_varid,start,
-     s  count,temp)
-#else
-        status=NF_PUT_VARA_REAL(ncid_out,out_varid,start,
-     s  count,temp)
-#endif
-
-
-c fermeture du fichier netcdf
-        call ncclos(ncid_out,rcode_out)
-        write(*,*) 'Fermeture: ',fich_out
-
-
-      print*,'OK5'
-c   Ecriture grads
-      open (20,file='grille.dat',form='unformatted',access='direct'
-     s      ,recl=4*ip1jmp1)
-      write(20,rec=1) ((REAL(mod(i,2)+mod(j,2)),i=1,iip1),j=1,jjp1)
-      write(20,rec=2) ((REAL(mod(i,2)*mod(j,2)),i=1,iip1),j=1,jjp1)
-      do j=2,jjm
-         dlat1(j)=180.*(rlatv(j)-rlatv(j-1))/pi
-c        dlat2(j)=180.*fyprim(REAL(j))/pi
-      enddo
-      do i=2,iip1
-         dlon1(i)=180.*(rlonu(i)-rlonu(i-1))/pi
-c        dlon2(i)=180.*fxprim(REAL(i))/pi
-      enddo
-      do j=2,jjm
-         dykm=(rlatv(j)-rlatv(j-1))*6400.
-         acoslat=6400.*cos(rlatu(j))
-         do i=2,iip1
-            dxkm=acoslat*(rlonu(i)-rlonu(i-1))
-            resol(i,j)=sqrt(dykm*dykm+dxkm*dxkm)
-         enddo
-         resol(1,j)=resol(iip1,j)
-      enddo
-      write(20,rec=3) resol
-      dlon1(1)=dlon1(iip1)
-      dlon2(1)=dlon2(iip1)
-      write(20,rec=4) ((dlon1(i),i=1,iip1),j=1,jjp1)
-      write(20,rec=5) ((dlon1(i)*pi/180.*0.001*
-     s   cos(rlatu(j))*rad,i=1,iip1),j=1,jjp1)
-      write(20,rec=6) ((dlon2(i),i=1,iip1),j=1,jjp1)
-      write(20,rec=7) ((dlat1(j),i=1,iip1),j=1,jjp1)
-      write(20,rec=8) ((dlat1(j)*pi/180.*rad*0.001,i=1,iip1),j=1,jjp1)
-      write(20,rec=9) ((dlat2(j),i=1,iip1),j=1,jjp1)
-
-      print*,'I, LON, DX (km)'
-      do i=1,iip1
-         print*,i,rlonu(i)*180./pi,dlon1(i)*pi/180.*0.001*
-     s   cos(clat*pi/180.)*rad
-      enddo
-      print*,'J, LAT, DY (km)'
-      do j=1,jjp1
-         print*,j,rlatu(j)*180./pi,dlat1(j)*pi/180.*0.001*rad
-      enddo
-
-      open (21,file='grille.ctl',form='formatted')
-
-c   WARNING! on reecrase le fichier .ctl a chaque ecriture
-      write(21,'(a5,1x,a40)')
-     &       'DSET ','^grille.dat'
-
-      write(21,'(a12)') 'UNDEF 1.0E30'
-      write(21,'(a5,1x,a40)') 'TITLE ','grille'
-      call formcoord(21,iip1,rlonv,180./pi,.false.,'XDEF')
-      call formcoord(21,jjp1,rlatu,180./pi,.true.,'YDEF')
-      call formcoord(21,1,0.,1.,.false.,'ZDEF')
-      write(21,'(a4,i10,a30)')
-     &       'TDEF ',1,' LINEAR 23OCT1994 3hr '
-      write(21,'(a4,2x,i5)') 'VARS',9
-      write(21,'(a18)') 'grille 0 99 grille'
-      write(21,'(a18)') 'gril   0 99 gril  '
-      write(21,'(a29)') 'resol   0 99 resolution (km)  '
-      write(21,'(a18)') 'dlon1  0 99 dlon1 '
-      write(21,'(a20)') 'dx     0 99 dx (km) '
-      write(21,'(a18)') 'dlon2  0 99 dlon2 '
-      write(21,'(a18)') 'dlat1  0 99 dlat1 '
-      write(21,'(a20)') 'dy     0 99 dy (km) '
-      write(21,'(a18)') 'dlat2  0 99 dlat2 '
-      write(21,'(a7)') 'ENDVARS'
-
-
-
-
-
-      print*,'OK6'
-	end
-
-
-
-        subroutine handle_err(status)
-#include "netcdf.inc"
-
-
-        integer status
-        print *,'handle code err: ',NF_NOERR
-        IF (status.NE.nf_noerr) THEN
-                print *,NF_STRERROR(status)
-                stop 'stopped'
-        ENDIF
-        END
-
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90	(revision 1302)
@@ -14,4 +14,6 @@
 #endif
   USE Write_Field
+  use exner_hyb_m, only: exner_hyb
+  use exner_milieu_m, only: exner_milieu
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -54,5 +56,4 @@
   REAL pks(ip1jmp1)                      ! exner au  sol
   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
   REAL phi(ip1jmp1,llm)                  ! geopotentiel
   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
@@ -223,7 +224,7 @@
         CALL pression ( ip1jmp1, ap, bp, ps, p       )
         if (pressure_exner) then
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        else
-          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
+        else
+          call exner_milieu(ip1jmp1,ps,p,pks,pk)
         endif
         CALL massdair(p,masse)
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/iniconst.F90	(revision 1302)
@@ -73,5 +73,5 @@
   if (disvert_type==1) then
      ! standard case for Earth (automatic generation of levels)
-     call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig, scaleheight)
+     call disvert()
   else if (disvert_type==2) then
      ! standard case for planets (levels generated using z2sig.def file)
Index: unk/LMDZ.COMMON/libf/dyn3d_common/pres2lev_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/pres2lev_mod.F90	(revision 1301)
+++ 	(revision )
@@ -1,79 +1,0 @@
-! $Id: pres2lev.F 1179 2009-06-11 14:18:47Z jghattas $
-!
-MODULE pres2lev_mod
-
-CONTAINS 
-
-!******************************************************
-SUBROUTINE pres2lev(varo,varn,lmo,lmn,po,pn,ni,nj,ok_invertp)
-!
-! interpolation lineaire pour passer
-! a une nouvelle discretisation verticale pour
-! les variables de GCM
-! Francois Forget (01/1995)
-! MOdif remy roca 12/97 pour passer de pres2sig
-! Modif F.Codron 07/08 po en 3D
-!**********************************************************
-
-  IMPLICIT NONE
-
-!   Declarations:
-! ==============
-!
-!  ARGUMENTS
-!  """""""""
-  LOGICAL, INTENT(IN) :: ok_invertp
-  INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
-  INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
-  
-  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
-  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
-  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
-
-  REAL, INTENT(IN)  :: varo(ni*nj,lmo) ! var dans l'ancienne grille
-  REAL, INTENT(OUT) :: varn(ni*nj,lmn) ! var dans la nouvelle grille
-
-  REAL :: zvaro(ni*nj,lmo),zpo(ni*nj,lmo)
-
-! Autres variables
-! """"""""""""""""
-  INTEGER ::  ln ,lo, k
-  REAL    :: coef
-
-
-! Inversion de l'ordre des niveaux verticaux
-  IF (ok_invertp) THEN
-    DO lo=1,lmo
-      DO k=1,ni*nj
-        zpo(k,lo)=po(k,lmo+1-lo)
-        zvaro(k,lo)=varo(k,lmo+1-lo)
-      ENDDO
-    ENDDO
-  ELSE
-    DO lo=1,lmo
-      DO k=1,ni*nj
-        zpo(k,lo)=po(k,lo)
-        zvaro(k,lo)=varo(k,lo)
-      ENDDO
-    ENDDO
-  ENDIF 
-
-  DO ln=1,lmn
-    DO lo=1,lmo-1
-      DO k=1,ni*nj
-        IF (pn(k,ln) >= zpo(k,1) ) THEN
-          varn(k,ln) = zvaro(k,1)
-        ELSE IF (pn(k,ln) <= zpo(k,lmo)) THEN
-          varn(k,ln) = zvaro(k,lmo)
-        ELSE IF ( pn(k,ln) <= zpo(k,lo) .AND. pn(k,ln) > zpo(k,lo+1) ) THEN
-          coef = (pn(k,ln)-zpo(k,lo)) / (zpo(k,lo+1)-zpo(k,lo))
-          varn(k,ln) = zvaro(k,lo) + coef*(zvaro(k,lo+1)-zvaro(k,lo))
-        ENDIF
-         
-      ENDDO  
-    ENDDO
-  ENDDO                
-
-END SUBROUTINE pres2lev    
-
-END MODULE pres2lev_mod
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/q_sat.F	(revision 1302)
@@ -2,39 +2,39 @@
 ! $Header$
 !
-c
-c
+!
+!
 
       subroutine q_sat(np,temp,pres,qsat)
-c
+!
       IMPLICIT none
-c======================================================================
-c Autheur(s): Z.X. Li (LMD/CNRS)
-c  reecriture vectorisee par F. Hourdin.
-c Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
-c======================================================================
-c Arguments:
-c kelvin---input-R: temperature en Kelvin
-c millibar--input-R: pression en mb
-c
-c q_sat----output-R: vapeur d'eau saturante en kg/kg
-c======================================================================
-c
+!======================================================================
+! Autheur(s): Z.X. Li (LMD/CNRS)
+!  reecriture vectorisee par F. Hourdin.
+! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
+!======================================================================
+! Arguments:
+! kelvin---input-R: temperature en Kelvin
+! millibar--input-R: pression en mb
+!
+! q_sat----output-R: vapeur d'eau saturante en kg/kg
+!======================================================================
+!
       integer np
       REAL temp(np),pres(np),qsat(np)
-c
+!
       REAL r2es
       PARAMETER (r2es=611.14 *18.0153/28.9644)
-c
+!
       REAL r3les, r3ies, r3es
       PARAMETER (R3LES=17.269)
       PARAMETER (R3IES=21.875)
-c
+!
       REAL r4les, r4ies, r4es
       PARAMETER (R4LES=35.86)
       PARAMETER (R4IES=7.66)
-c
+!
       REAL rtt
       PARAMETER (rtt=273.16)
-c
+!
       REAL retv
       PARAMETER (retv=28.9644/18.0153 - 1.0)
@@ -42,14 +42,14 @@
       real zqsat
       integer ip
-c
-C     ------------------------------------------------------------------
-c
-c
+!
+!     ------------------------------------------------------------------
+!
+!
 
       do ip=1,np
 
-c      write(*,*)'kelvin,millibar=',kelvin,millibar
-c       write(*,*)'temp,pres=',temp(ip),pres(ip)
-c
+!      write(*,*)'kelvin,millibar=',kelvin,millibar
+!       write(*,*)'temp,pres=',temp(ip),pres(ip)
+!
          IF (temp(ip) .LE. rtt) THEN
             r3es = r3ies
@@ -59,14 +59,15 @@
             r4es = r4les
          ENDIF
-c
+!
          zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
          zqsat=MIN(0.5,ZQSAT)
          zqsat=zqsat/(1.-retv *zqsat)
-c
+!
          qsat(ip)= zqsat
-c      write(*,*)'qsat=',qsat(ip)
+!      write(*,*)'qsat=',qsat(ip)
 
       enddo
-c
+!
       RETURN
       END
+
Index: unk/LMDZ.COMMON/libf/dyn3d_common/tourabs.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/tourabs.F	(revision 1301)
+++ 	(revision )
@@ -1,98 +1,0 @@
-      SUBROUTINE tourabs ( ntetaSTD,vcov, ucov, vorabs )
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c   Modif:  I. Musat (28/10/04)
-c   -------
-c   adaptation du code tourpot.F pour le calcul de la vorticite absolue
-c   cf. P. Le Van
-c
-c   Objet: 
-c   ------
-c
-c    *******************************************************************
-c    .............  calcul de la vorticite absolue     .................
-c    *******************************************************************
-c
-c     ntetaSTD, vcov,ucov      sont des argum. d'entree pour le s-pg .
-c             vorabs            est  un argum.de sortie pour le s-pg .
-c
-c=======================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "logic.h"
-#include "comconst.h"
-c
-      INTEGER ntetaSTD
-      REAL vcov( ip1jm,ntetaSTD ), ucov( ip1jmp1,ntetaSTD )
-      REAL vorabs( ip1jm,ntetaSTD )
-c
-c variables locales
-      INTEGER l, ij, i, j
-      REAL  rot( ip1jm,ntetaSTD )
-
-
-
-c  ... vorabs = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) ..
-
-
-
-c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
-
-      DO 5 l = 1,ntetaSTD
-
-      DO 2 i = 1, iip1
-      DO 2 j = 1, jjm
-c
-       ij=i+(j-1)*iip1
-       IF(ij.LE.ip1jm - 1) THEN
-c
-        IF(cv(ij).EQ.0..OR.cv(ij+1).EQ.0..OR.
-     $     cu(ij).EQ.0..OR.cu(ij+iip1).EQ.0.) THEN
-         rot( ij,l ) = 0.
-         continue
-        ELSE
-         rot( ij,l ) = (vcov(ij+1,l)/cv(ij+1)-vcov(ij,l)/cv(ij))/
-     $                 (2.*pi*RAD*cos(rlatv(j)))*REAL(iim)
-     $                +(ucov(ij+iip1,l)/cu(ij+iip1)-ucov(ij,l)/cu(ij))/
-     $                 (pi*RAD)*(REAL(jjm)-1.)
-c
-        ENDIF
-       ENDIF !(ij.LE.ip1jm - 1) THEN
-   2  CONTINUE
-
-c    ....  correction pour  rot( iip1,j,l )  .....
-c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
-
-CDIR$ IVDEP
-
-      DO 3 ij = iip1, ip1jm, iip1
-      rot( ij,l ) = rot( ij -iim, l )
-   3  CONTINUE
-
-   5  CONTINUE
-
-
-      CALL  filtreg( rot, jjm, ntetaSTD, 2, 1, .FALSE., 1 )
-
-
-      DO 10 l = 1, ntetaSTD
-
-      DO 6 ij = 1, ip1jm - 1
-      vorabs( ij,l ) = ( rot(ij,l) + fext(ij)*unsairez(ij) )
-   6  CONTINUE
-
-c    ..... correction pour  vorabs( iip1,j,l)  .....
-c    ....   vorabs(iip1,j,l)= vorabs(1,j,l) ....
-CDIR$ IVDEP
-      DO 8 ij = iip1, ip1jm, iip1
-      vorabs( ij,l ) = vorabs( ij -iim,l )
-   8  CONTINUE
-
-  10  CONTINUE
-
-      RETURN
-      END
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1302)
@@ -225,10 +225,4 @@
       save unskap
 
-cIM diagnostique PVteta, Amip2
-      INTEGER,PARAMETER :: ntetaSTD=3
-      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
-      REAL PVteta(klon,ntetaSTD)
-      
-      
       REAL SSUM
 
@@ -267,6 +261,4 @@
       klon=klon_mpi
       
-      PVteta(:,:)=0.
-            
       IF ( firstcal )  THEN
         debut = .TRUE.
@@ -684,16 +676,4 @@
       endif
 
-
-      IF (is_sequential.and.(planet_type=="earth")) THEN
-#ifdef CPP_EARTH
-! PVtheta calls tetalevel, which is in the physics
-cIM calcul PV a teta=350, 380, 405K
-        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
-     $           ztfi,zplay,zplev,
-     $           ntetaSTD,rtetaSTD,PVteta)
-c
-#endif
-      ENDIF
-
 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
       CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
@@ -923,7 +903,5 @@
      .             zdqfi_omp,
      .             zdpsrf_omp,
-cIM diagnostique PVteta, Amip2          
-     .             pducov,
-     .             PVteta)
+     .             pducov)
 
       else if ( planet_type=="generic" ) then
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 1302)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1615 2012-02-10 15:42:26Z emillour $
+! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -115,12 +115,10 @@
   END IF
 
-  IF (grilles_gcm_netcdf) THEN
-     WRITE(lunout,'(//)')
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,'(//)')
-     CALL grilles_gcm_netcdf_sub(masque,phis)
-  END IF
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
   
 #ifdef CPP_MPI
@@ -137,2 +135,3 @@
 !
 !-------------------------------------------------------------------------------
+
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1302)
@@ -2,8 +2,8 @@
 ! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit )
-c
+!
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -20,18 +20,18 @@
       use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
       IMPLICIT NONE
-c-----------------------------------------------------------------------
-c     Auteurs :   L. Fairhead , P. Le Van  .
-c
-c     Arguments :
-c
-c     tapedef   :
-c     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
-c     -metres  du zoom  avec  celles lues sur le fichier start .
-c
+!-----------------------------------------------------------------------
+!     Auteurs :   L. Fairhead , P. Le Van  .
+!
+!     Arguments :
+!
+!     tapedef   :
+!     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
+!     -metres  du zoom  avec  celles lues sur le fichier start .
+!
        LOGICAL etatinit
        INTEGER tapedef
 
-c   Declarations :
-c   --------------
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -45,8 +45,8 @@
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -60,32 +60,32 @@
       integer,external :: OMP_GET_NUM_THREADS
 #endif      
-c
-c  -------------------------------------------------------------------
-c
-c       .........     Version  du 29/04/97       ..........
-c
-c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
-c      tetatemp   ajoutes  pour la dissipation   .
-c
-c   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
-c
-c  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
-c    Sinon , choix de fxynew  , a derivee sinusoidale  ..
-c
-c   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
-c         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
-c                de limit.dat ( dic)                        ...........
-c           Sinon  etatinit = . FALSE .
-c
-c   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
-c    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
-c   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
-c    lectba .  
-c   Ces parmetres definissant entre autres la grille et doivent etre
-c   pareils et coherents , sinon il y aura  divergence du gcm .
-c
-c-----------------------------------------------------------------------
-c   initialisations:
-c   ----------------
+!
+!  -------------------------------------------------------------------
+!
+!       .........     Version  du 29/04/97       ..........
+!
+!   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
+!      tetatemp   ajoutes  pour la dissipation   .
+!
+!   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
+!
+!  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
+!    Sinon , choix de fxynew  , a derivee sinusoidale  ..
+!
+!   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
+!         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
+!                de limit.dat ( dic)                        ...........
+!           Sinon  etatinit = . FALSE .
+!
+!   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
+!    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
+!   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
+!    lectba .  
+!   Ces parmetres definissant entre autres la grille et doivent etre
+!   pareils et coherents , sinon il y aura  divergence du gcm .
+!
+!-----------------------------------------------------------------------
+!   initialisations:
+!   ----------------
 
 !Config  Key  = lunout
@@ -290,6 +290,6 @@
        CALL getin('dissip_period',dissip_period)
 
-ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
-ccc
+!cc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
+!cc
 
 !Config  Key  = lstardis
@@ -456,5 +456,5 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!     ...............................................................
 
 !Config  Key  =  read_start
@@ -632,12 +632,6 @@
       CALL getin('ok_etat0',ok_etat0)
 
-!Config  Key  = grilles_gcm_netcdf
-!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
-!Config  Def  = n
-      grilles_gcm_netcdf = .FALSE.
-      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
-
-c----------------------------------------
-c Parameters for zonal averages in the case of Titan
+!----------------------------------------
+! Parameters for zonal averages in the case of Titan
       moyzon_mu = .false.
       moyzon_ch = .false.
@@ -646,13 +640,13 @@
        CALL getin('moyzon_ch', moyzon_ch)
       endif
-c----------------------------------------
-
-c----------------------------------------
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
-C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
-c
-c----------------------------------------
+!----------------------------------------
+
+!----------------------------------------
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
+! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
+!
+!----------------------------------------
       IF( etatinit ) then
 
@@ -707,8 +701,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -758,8 +752,8 @@
        ysinus = .TRUE.
        CALL getin('ysinus',ysinus)
-c
-c----------------------------------------
+!
+!----------------------------------------
        else ! etatinit=false
-c----------------------------------------
+!----------------------------------------
 
 !Config  Key  = clon
@@ -779,6 +773,6 @@
        CALL getin('clat',clatt)
 
-c
-c
+!
+!
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
         write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
@@ -834,8 +828,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -862,5 +856,5 @@
          ENDIF
       ENDIF
-c
+!
 !Config  Key  = dzoomx
 !Config  Desc = extension en longitude
@@ -925,5 +919,5 @@
       ENDIF
 
-cc
+!c
       IF( .NOT.fxyhypb  )  THEN
 
@@ -955,5 +949,5 @@
 
       endif ! etatinit
-c----------------------------------------
+!----------------------------------------
 
 
@@ -1009,5 +1003,4 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
       if (planet_type=="titan") then
        write(lunout,*)' moyzon_mu = ', moyzon_mu
Index: unk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 1301)
+++ 	(revision )
@@ -1,231 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
-c
-c     Auteurs :  P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c                                 -------- z                                   
-c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
-c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
-c    ( voir note de Fr.Hourdin )  ,
-c
-c    on determine successivement , du haut vers le bas des couches, les 
-c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
-c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
-c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
-c
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL unpl2k,dellta
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_hyb_p"
-c
-
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-c$OMP BARRIER
-
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
-
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-      unpl2k    = 1.+ 2.* kappa
-c
-      ijb=ij_begin
-      ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
-
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO     ij      = ijb,ije
-       alpha(ij,llm) = 0.
-       beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-c$OMP DO SCHEDULE(STATIC)
-        DO ij = ijb, ije
-        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
-        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
-        beta (ij,l)  =   p(ij,l  ) / dellta   
-        ENDDO
-c$OMP ENDDO NOWAIT
-c
-      ENDDO
-
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
-     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-c
-c
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-
-      RETURN
-      END
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 1302)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 1302)
@@ -0,0 +1,187 @@
+module exner_hyb_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
+
+    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+    REAL alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    INTEGER ije,ijb,jje,jjb
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_hyb_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    !$OMP BARRIER
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
+
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    ijb=ij_begin
+    ije=ij_end
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
+
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO     ij      = ijb,ije
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij = ijb, ije
+          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
+          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
+          beta (ij,l)  =   p(ij,l  ) / dellta   
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
+            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+       !$OMP ENDDO NOWAIT        
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT             
+       ENDDO
+
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb_p
+
+end module exner_hyb_p_m
+
Index: unk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F	(revision 1301)
+++ 	(revision )
@@ -1,218 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL dum1
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_milieu_p"
-
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-      
-c$OMP BARRIER
-
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-              
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
-
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      ijb=ij_begin
-      ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
-
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-c$OMP ENDDO NOWAIT        
-
-
-c    calcul de pkf
-c    -------------
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-      RETURN
-      END
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 1302)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 1302)
@@ -0,0 +1,160 @@
+module exner_milieu_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_milieu_p ( ngrid, ps, p, pks, pk, pkf )
+    !
+    !     Auteurs :  F. Forget , Y. Wanherdrick
+    ! P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !     WARNING : CECI est une version speciale de exner_hyb originale
+    !               Utilise dans la version martienne pour pouvoir 
+    !               tourner avec des coordonnees verticales complexe
+    !              => Il ne verifie PAS la condition la proportionalite en 
+    !              energie totale/ interne / potentielle (F.Forget 2001)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij,ijb,ije,jjb,jje
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_milieu_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    !$OMP BARRIER
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
+
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    ijb=ij_begin
+    ije=ij_end
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
+
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !    .... Calcul de pk  pour la couche l = llm ..
+    !    (on met la meme distance (en log pression)  entre Pk(llm)
+    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
+    !$OMP ENDDO NOWAIT        
+
+    if (present(pkf)) then
+       !    calcul de pkf
+
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT
+       ENDDO
+
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu_p
+
+end module exner_milieu_p_m
+
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1302)
@@ -516,5 +516,9 @@
 
 
-      day_end = day_ini + nday
+      if (nday>=0) then ! standard case
+        day_end=day_ini+nday
+      else ! special case when nday <0, run -nday dynamical steps
+        day_end=day_ini-nday/day_step
+      endif
       if (less1day) then
         day_end=day_ini+floor(time_0+fractday)
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 1302)
@@ -328,4 +328,6 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
     USE parallel_lmdz
     USE control_mod
@@ -349,6 +351,5 @@
     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
     ! Variables pour fonction Exner (P milieu couche)
-    REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-    REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+    REAL, DIMENSION (iip1,jjp1,llm)    :: pk
     REAL, DIMENSION (iip1,jjp1)        :: pks    
     REAL                               :: unskap
@@ -491,19 +492,19 @@
     IF (f_out) THEN
 !       Calcul niveaux pression milieu de couches 
-	CALL pression_p( ip1jmp1, ap, bp, ps, p )
-	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
-	else
-          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        CALL pression_p( ip1jmp1, ap, bp, ps, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
+        else
+          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
         endif
         unskap=1./kappa
-	DO l = 1, llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
-		    p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-		ENDDO
-	    ENDDO
-	ENDDO
-        CALL guide_out("P",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
+        DO l = 1, llm
+           DO j=jjb_u,jje_u
+              DO i =1, iip1
+                 p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+              ENDDO
+           ENDDO
+        ENDDO
+        CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
     ENDIF
     
@@ -517,5 +518,7 @@
         if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
-        IF (f_out) CALL guide_out("U",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
+        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(:,:)+tau*ugui2(:,:),factt)
+        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
         ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -529,5 +532,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
-        IF (f_out) CALL guide_out("T",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
         teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -541,5 +544,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
-        IF (f_out) CALL guide_out("SP",jjp1,1,f_add(1:ip1jmp1,1),factt)
+        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
         ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
         CALL pression_p(ip1jmp1,ap,bp,ps,p)
@@ -555,5 +558,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
-        IF (f_out) CALL guide_out("Q",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
         q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -568,5 +571,7 @@
         if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
-        IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:),factt)
+        IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
+        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1(:,:)+tau*vgui2(:,:),factt)
+        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
         vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
     endif
@@ -689,4 +694,6 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
   USE parallel_lmdz
   USE mod_hallo
@@ -713,6 +720,5 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-  REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: unskap
@@ -784,25 +790,25 @@
     IF (guide_plevs.EQ.1) THEN
         DO l=1,llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
                     pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
-		ENDDO
-	    ENDDO
+                ENDDO
+            ENDDO
         ENDDO
     ELSE
-	CALL pression_p( ip1jmp1, ap, bp, psi, p )
-	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+        CALL pression_p( ip1jmp1, ap, bp, psi, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
         else
-          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
         endif
-	unskap=1./kappa
-	DO l = 1, llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
-		    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-		ENDDO
-	    ENDDO
-	ENDDO
+        unskap=1./kappa
+        DO l = 1, llm
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
+                    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+                ENDDO
+            ENDDO
+        ENDDO
     ENDIF
 
@@ -1024,10 +1030,10 @@
         ! Calcul des nouvelles valeurs des niveaux de pression du guidage
         IF (guide_plevs.EQ.1) THEN
-	 CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
-	 CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
             DO l=1,nlevnc
                 DO j=jjb_v,jje_v
@@ -1041,10 +1047,10 @@
             ENDDO
         ELSE IF (guide_plevs.EQ.2) THEN
-	 CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
-	 CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
             DO l=1,nlevnc
                 DO j=jjb_v,jje_v
@@ -1807,5 +1813,5 @@
     
     ! Variables entree
-    CHARACTER, INTENT(IN)                          :: varname
+    CHARACTER*(*), INTENT(IN)                          :: varname
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
@@ -1817,7 +1823,10 @@
     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
+    INTEGER       :: vid_au,vid_av
+    INTEGER       :: l
     INTEGER, DIMENSION (3) :: dim3
     INTEGER, DIMENSION (4) :: dim4,count,start
     INTEGER                :: ierr, varid
+    REAL, DIMENSION (iip1,hsize,vsize) :: field2
     
     CALL gather_field(field,iip1*hsize,vsize,0)
@@ -1834,4 +1843,5 @@
 ! Definition des dimensions
         ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu) 
+        print*,'id_lonu 1 ',id_lonu
         ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv) 
         ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu) 
@@ -1842,4 +1852,5 @@
 ! Creation des variables dimensions
         ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
+        print*,'id_lonu 2 ',id_lonu
         ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
         ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
@@ -1848,4 +1859,6 @@
         ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
         ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
+        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
         
         ierr=NF_ENDDEF(nid)
@@ -1853,4 +1866,5 @@
 ! Enregistrement des variables dimensions
 #ifdef NC_DOUBLE
+        print*,'id_lonu DOUBLE ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
@@ -1860,5 +1874,8 @@
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
 #else
+        print*,'id_lonu 3 ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
@@ -1868,4 +1885,6 @@
         ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
 #endif
 ! --------------------------------------------------------------------
@@ -1875,5 +1894,5 @@
 ! Pressure (GCM)
         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = NF_DEF_VAR(nid,"P",NF_FLOAT,4,dim4,varid)
+        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
@@ -1883,5 +1902,8 @@
 ! Zonal wind
         IF (guide_u) THEN
+        print*,'id_lonu 4 ',id_lonu,varname
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1889,4 +1911,6 @@
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1912,61 +1936,41 @@
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
+    IF (varname=="SP") timestep=timestep+1
+
+    IF (varname=="SP") THEN
+      print*,'varname=SP=',varname
+    ELSE
+      print*,'varname diff SP=',varname
+    ENDIF
+
+
+    ierr = NF_INQ_VARID(nid,varname,varid)
     SELECT CASE (varname)
-    CASE ("P")
-        timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"P",varid)
+    CASE ("SP","ps")
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("SP")
-        ierr = NF_INQ_VARID(nid,"ps",varid)
-        start=(/1,1,timestep,0/)
-        count=(/iip1,jjp1,1,0/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("U")
-        ierr = NF_INQ_VARID(nid,"ucov",varid)
+    CASE ("v","va","vcov")
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjm,llm,1/)
+    CASE DEFAULT
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
+    END SELECT
+
+    SELECT CASE (varname)
+    CASE("u","ua")
+        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
+        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
+    CASE("v","va")
+        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE DEFAULT
+        field2=field
+    END SELECT
+
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
 #else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
 #endif
-    CASE ("V")
-        ierr = NF_INQ_VARID(nid,"vcov",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjm,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("T")
-        ierr = NF_INQ_VARID(nid,"teta",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("Q")
-        ierr = NF_INQ_VARID(nid,"q",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    END SELECT
  
     ierr = NF_CLOSE(nid)
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1302)
@@ -8,4 +8,8 @@
      &                    time_0)
 
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
+      use exner_hyb_p_m, only: exner_hyb_p
+      use exner_milieu_p_m, only: exner_milieu_p
        USE misc_mod
        USE parallel_lmdz
@@ -17,5 +21,5 @@
        USE vampir
        USE timer_filtre, ONLY : print_filtre_timer
-       USE infotrac
+       USE infotrac, ONLY: nqtot
        USE guide_p_mod, ONLY : guide_main
        USE getparam
@@ -165,5 +169,4 @@
       character*10 string10
 
-      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
       REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
 
@@ -236,5 +239,10 @@
       lafin=.false.
       
-      itaufin   = nday*day_step
+      if (nday>=0) then
+         itaufin   = nday*day_step
+      else
+         ! to run a given (-nday) number of dynamical steps
+         itaufin   = -nday
+      endif
       if (less1day) then
 c MODIF VENUS: to run less than one day:
@@ -292,7 +300,7 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 c$OMP END MASTER
@@ -830,8 +838,10 @@
 c$OMP BARRIER
          if (pressure_exner) then
-           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
          endif
+! Compute geopotential (physics might need it)
+         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 c$OMP BARRIER
            jD_cur = jD_ref + day_ini - day_ref
@@ -1040,7 +1050,7 @@
 c$OMP BARRIER
           if (pressure_exner) then
-            CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
           else
-            CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
           endif
 c$OMP BARRIER
@@ -1202,7 +1212,7 @@
 c$OMP BARRIER
         if (pressure_exner) then
-          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
         endif
 c$OMP BARRIER
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h	(revision 1302)
@@ -11,5 +11,5 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
+     &  ,ok_limit,ok_etat0,hybrid                                       &
      &  ,moyzon_mu,moyzon_ch
 
@@ -19,5 +19,5 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0
       logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                      ! (only used if disvert_type==2)
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 1302)
@@ -1,4 +1,4 @@
 ! 
-! $Id: mod_const_mpi.F90 1700 2012-12-20 14:43:19Z lguez $
+! $Id: mod_const_mpi.F90 2055 2014-06-04 12:33:27Z acaubel $
 !
 MODULE mod_const_mpi
@@ -17,5 +17,11 @@
     USE ioipsl_getincom, only: getin
 #endif
-
+! Use of Oasis-MCT coupler 
+#ifdef CPP_OMCT
+    USE mod_prism
+#endif
+#ifdef CPP_XIOS
+    USE wxios, only: wxios_init
+#endif
     IMPLICIT NONE
 #ifdef CPP_MPI
@@ -38,6 +44,10 @@
 #ifdef CPP_COUPLE
 !$OMP MASTER
-       CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
+#ifdef CPP_XIOS
+      CALL wxios_init("LMDZ", outcom=COMM_LMDZ, type_ocean=type_ocean)
+#else
+       CALL prism_init_comp_proto (comp_id, 'LMDZ', ierr)
        CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
+#endif
 !$OMP END MASTER
 #endif
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90	(revision 1301)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90	(revision 1302)
@@ -225,7 +225,10 @@
 #endif
 #ifdef CPP_COUPLE
+! Use of Oasis-MCT coupler 
+#if defined CPP_OMCT
+    use mod_prism
+#else
     use mod_prism_proto
 #endif
-#ifdef CPP_EARTH
 ! Ehouarn: surface_data module is in 'phylmd' ...
       use surface_data, only : type_ocean
@@ -252,4 +255,8 @@
 
       if (type_ocean == 'couple') then
+#ifdef CPP_XIOS
+    !Fermeture propre de XIOS
+      CALL wxios_close()
+#else
 #ifdef CPP_COUPLE
          call prism_terminate_proto(ierr)
@@ -258,4 +265,5 @@
          endif
 #endif 
+#endif
       else
 #ifdef CPP_XIOS
Index: /trunk/LMDZ.COMMON/makelmdz
===================================================================
--- /trunk/LMDZ.COMMON/makelmdz	(revision 1301)
+++ /trunk/LMDZ.COMMON/makelmdz	(revision 1302)
@@ -26,7 +26,7 @@
 cosp=false
 sisvat=false
+rrtm=false
 bands=""
 scatterers=""
-dyn=""
 full=''
 
@@ -35,15 +35,7 @@
 ## try to recognise machine and infer arch from it
 machine=`hostname`
-if [[ "$machine" == "brodie" ]]
-then
-  arch="SX8_BRODIE"
-fi
 if [[ "${machine:0:3}" == "ada" ]]
 then
   arch="X64_ADA"
-fi
-if [[ "${machine:0:6}" == "ciclad" ]]
-then
-  arch="AMD64_CICLAD"
 fi
 if [[ "${machine:0:7}" == "platine" ]]
@@ -87,5 +79,5 @@
 
 CPP_KEY="" 
-INCLUDE='-I$(LIBF)/grid -I$(LIBF)/bibio -I$(LIBF)/dyn3d_common -I$(LIBF)/filtrez -I. '
+INCLUDE='-I$(LIBF)/grid -I$(LIBF)/bibio -I$(LIBF)/filtrez -I. '
 LIB=""
 adjnt=""
@@ -184,9 +176,4 @@
 	  CPP_KEY="$CPP_KEY $2" ; shift ; shift ;;
 
-      "-adjnt")
-	  echo "not operational ... work to be done here ";exit
-	  opt_dep="$opt_dep adjnt" ; adjnt="-ladjnt -ldyn3d "
-	  optim="$optim -Dadj" ; shift ;;
-
       "-cosp")
           cosp="$2" ; shift ; shift ;;
@@ -194,4 +181,7 @@
       "-sisvat")
           sisvat="$2" ; shift ; shift ;;
+      
+      "-rrtm")
+          rrtm="$2" ; shift ; shift ;;
       
       "-mem")
@@ -457,4 +447,9 @@
 fi
 
+if [[ "$rrtm" == "true" ]]
+then
+   CPP_KEY="$CPP_KEY CPP_RRTM"
+fi
+
 
 INCLUDE="$INCLUDE ${NETCDF_INCDIR}"
@@ -536,67 +531,20 @@
 echo dimension $dimension
 
-if (( $dimension == 3 ))
-then
+if (( $dimension == 3 )) ; then
   cd $LIBFGCM/grid
   \rm fxyprim.h
   cp -p fxy_${grille}.h fxyprim.h
-#else
-#  echo "Probleme dans les dimensions de la dynamique !!"
-#  echo "Non reactive pour l'instant !!!"
-fi
-if (( $dimension == 1 ))
-then
-  echo pas de dynamique
-  dyn="DYN= L_DYN=-ldyn3d "
-fi
-
-######################################################################
-# Gestion du filtre qui n'existe qu'en 3d.
-######################################################################
-
-if (( `expr $dimc \> 2` == 1 ))
-then
-   filtre="FILTRE=$filtre"
+  filtre="FILTRE=$filtre"
+  INCLUDE="$INCLUDE "'-I$(LIBF)/dyn3d${FLAG_PARA} -I$(LIBF)/dyn3d_common '
+elif (( $dimension == 2 )) ; then
+   filtre="FILTRE= L_FILTRE= "
+   INCLUDE="$INCLUDE "'-I$(LIBF)/dyn2d'
+elif (( $dimension == 1 )) ; then
+   echo pas de dynamique
+   CPP_KEY="$CPP_KEY CPP_1D"
+   filtre="L_DYN= DYN= FILTRE= L_FILTRE= DIRMAIN=phy$physique "
+   INCLUDE="$INCLUDE "'-I$(LIBF)/dyn3d -I$(LIBF)/dyn3d_common ' # Pas tres propre
 else
-   filtre="FILTRE= L_FILTRE= "
-fi
-echo MACRO FILTRE $filtre
-
-echo $dimc
-
-######################################################################
-# Creation du suffixe de la configuration
-######################################################################
-
-
-SUFF_NAME=_${dim_full}
-SUFF_NAME=${SUFF_NAME}_phy${physique}
-
-if [[ "$parallel" != "none" ]]
-then
-  SUFF_NAME=${SUFF_NAME}_para
-  DYN=dyn${dimc}d${paramem}
-  if [[ "$paramem" == "mem" ]]
-  then
-   SUFF_NAME=${SUFF_NAME}_${paramem}
-  fi
-else
-  SUFF_NAME=${SUFF_NAME}_seq
-  DYN=dyn${dimc}d
-fi
-
-if [[ $veget != "false" ]]
-then
-  SUFF_NAME=${SUFF_NAME}_orch
-fi
-
-if [[ $couple != "false" ]]
-then
-  SUFF_NAME=${SUFF_NAME}_couple
-fi
-
-if [[ $chimie == "INCA" ]]
-then
-  SUFF_NAME=${SUFF_NAME}_inca
+   echo Dimension dimc=$dimc pas prevu ; exit
 fi
 
@@ -626,5 +574,5 @@
 find libf -name '*.[Fh]' -print | sort > tmp77
 find libf -name '*.[Fh]' -exec egrep -i " *use *ioipsl" {} \; -print | sort > tmp90
-find libf -name '*.[Fh90]' -print | sort >> tmp90
+find libf -name '*.F90' -print | sort >> tmp90
 ######### CAS PARTICULIER CHIMIE TITAN
 if [[ "$physique" == "titan" ]]
@@ -766,5 +714,4 @@
 INCLUDEC="$includec" \
 $filtre \
-$dyn \
 LIBO=${LIBOGCM}/${nomlib} \
 $phys \
@@ -798,5 +745,4 @@
 INCLUDEC="$includec" \
 $filtre \
-$dyn \
 LIBO=${LIBOGCM}/${nomlib} \
 $phys \
Index: /trunk/LMDZ.COMMON/makelmdz_fcm
===================================================================
--- /trunk/LMDZ.COMMON/makelmdz_fcm	(revision 1301)
+++ /trunk/LMDZ.COMMON/makelmdz_fcm	(revision 1302)
@@ -23,4 +23,5 @@
 veget=false
 sisvat=false
+rrtm=false
 chimie=false
 chemistry=false
@@ -40,4 +41,6 @@
 LIBFGCM=$LMDGCM/libf
 # path for optional packages, but default set to ".void_dir"
+RRTM_PATH=$LMDGCM/.void_dir
+SISVAT_PATH=$LMDGCM/.void_dir
 COSP_PATH=$LMDGCM/.void_dir
 CHEM_PATH=$LMDGCM/.void_dir
@@ -76,5 +79,5 @@
 [-p PHYS]                  : set of physical parametrizations (in libf/phyPHYS), (default: lmd)
 [-prod / -dev / -debug]    : compilation mode production (default) / developement / debug .
-[-c false/MPI1/MPI2]       : (Earth) coupling with ocean model : MPI1/MPI2/false (default: false)
+[-c false/MPI1/OMCT]       : coupling with ocean model : MPI1/OMCT/false (default: false)
 [-v false/orchidee2.0/orchidee1.9/true] : version of the vegetation model to include (default: false)
           false       : no vegetation model
@@ -138,4 +141,7 @@
 	  sisvat="$2" ; shift ; shift ;;
 
+      "-rrtm")
+	  rrtm="$2" ; shift ; shift ;;
+
       "-chimie")
 	  chimie="$2" ; shift ; shift ;;
@@ -250,7 +256,16 @@
 if [[ "$couple" != "false" ]]
 then
-   CPP_KEY="$CPP_KEY CPP_COUPLE"
-   INCLUDE="$INCLUDE -I${OASIS_INCDIR}"
-   LIB="$LIB -L${OASIS_LIBDIR} -lpsmile.${couple} -lmpp_io"
+    if [[ "$couple" == "MPI1" ]]
+    then
+	CPP_KEY="$CPP_KEY CPP_COUPLE"
+	export OASIS_INCDIR=$LMDGCM/../../prism/X64/build/lib/psmile.MPI1
+	export OASIS_LIBDIR=$LMDGCM/../../prism/X64/lib
+	INCLUDE="$INCLUDE -I${OASIS_INCDIR}"
+	LIB="$LIB -L${OASIS_LIBDIR} -lpsmile.MPI1 -lmpp_io"
+    else
+	CPP_KEY="$CPP_KEY CPP_COUPLE CPP_OMCT"
+	INCLUDE="$INCLUDE -I${OASIS_INCDIR}"
+	LIB="$LIB -L${OASIS_LIBDIR} -lpsmile.MPI1 -lscrip -lmct -lmpeu"
+    fi
 fi
 
@@ -310,6 +325,15 @@
 then
    CPP_KEY="$CPP_KEY CPP_SISVAT"
-   sed -e 's/^#src::sisvat/src::sisvat/' bld.cfg > bld.tmp
-   mv bld.tmp bld.cfg
+   SISVAT_PATH="$LIBFGCM/%PHYS/sisvat"
+   #sed -e 's/^#src::sisvat/src::sisvat/' bld.cfg > bld.tmp
+   #mv bld.tmp bld.cfg
+fi
+
+if [[ "$rrtm" == "true" ]]
+then
+   CPP_KEY="$CPP_KEY CPP_RRTM"
+   RRTM_PATH="$LIBFGCM/%PHYS/rrtm"
+   #sed -e 's/^#src::rrtm/src::rrtm/' bld.cfg > bld.tmp
+   #mv bld.tmp bld.cfg
 fi
 
@@ -536,4 +560,6 @@
 echo "%DYN           $DYN"           >> $config_fcm
 echo "%PHYS          phy${physique}" >> $config_fcm
+echo "%RRTM          $RRTM_PATH"     >> $config_fcm
+echo "%SISVAT        $SISVAT_PATH"   >> $config_fcm
 echo "%COSP          $COSP_PATH"     >> $config_fcm
 echo "%CHEM          $CHEM_PATH"     >> $config_fcm
@@ -567,6 +593,10 @@
 
 #eventual cleanup for SISVAT
-sed -e 's/^src::sisvat/#src::sisvat/' bld.cfg > bld.tmp
-mv bld.tmp bld.cfg
+#sed -e 's/^src::sisvat/#src::sisvat/' bld.cfg > bld.tmp
+#mv bld.tmp bld.cfg
+
+#eventual cleanup for RRTM
+#sed -e 's/^src::rrtm/#src::rrtm/' bld.cfg > bld.tmp
+#mv bld.tmp bld.cfg
 
 
