Index: trunk/DOC/chantiers/commit_importants.log
===================================================================
--- trunk/DOC/chantiers/commit_importants.log	(revision 978)
+++ trunk/DOC/chantiers/commit_importants.log	(revision 979)
@@ -1168,2 +1168,53 @@
   Now dynamics in seq/MPI/OpenMP/mixed modes yield identical results when
   running in Newtonian mode (i.e. dynamics only).
+
+*********************
+**** commit_v979 ****
+*********************
+Ehouarn: Updates in the dynamics (seq and //) to keep up with updates
+in LMDZ5 (up to LMDZ5 trunk, rev 1760):
+* General stuff (essentially to keep up with Earth model):
+- Updated makelmdz_fcm and makelmdz (more control on dimension.h, 
+  added option -mem, although it is not usefull for now)
+- Updated build_gcm with more control over fcm
+- Updated create_make_gcm (enable looking for code in subdirectories)
+
+* bibio:
+- updates (just renaming the files actually...) new_unit.F90 => new_unit_m.F90,
+  pchsp_95.F90 => pchsp_95_m.F90 and pchfe_95.F90 => pchfe_95_m.F90
+
+* filtrez:
+- mod_fft.F90: use more baseline CPP directives for preprocessor compatibility
+- mod_filtre_fft_loc.F90: added this new file
+- filtreg_mod.f90: added calls to init_..._loc
+- filtreg.F: fixed calls to DGEMM into SGEMM (preprocessing does the switch)
+
+* dyn3d:
+- removed obsolete files: etat0_netcdf.F90 limit_netcdf.F90
+                          pres2lev.F90
+- added new file : pres2lev_mod.F90 (module containing "old" pres2lev)
+- gcm.F: changed args to call to inidissip (added arg "vert_prof_dissip")
+- inidissip.F90: added arg "vert_prof_dissip" and the "earth model"
+  discterizations (flagged with "planet_type=="earth")
+- comdissnew.h: added 'vert_prof_dissip' to the common block
+- guide_mod.F90: added the "use pres2lev_mod"
+- conf_gcm.F: cosmetics, and evaluation of vert_prof_dissip,
+  (and also of dissip_* factors, for Earth model)
+- comconst.h : added dissip_factz,dissip_zref variables (for Earth mode dissip)
+
+* dyn3dpar:
+- removed obsolete files: etat0_netcdf.F90 limit_netcdf.F90
+                          pres2lev.F90 mod_const_para.F90
+- added new files: pres2lev_mod.F90 (module containing "old" pres2lev)
+                   mod_const_mpi
+- abort_gcm : better control of abort in parallel mode
+- gcm.F: changed args to call to inidissip (added arg "vert_prof_dissip")
+- inidissip.F90: added arg "vert_prof_dissip" and the "earth model"
+  discterizations (flagged with "planet_type=="earth")
+- comdissnew.h: added 'vert_prof_dissip' to the common block
+- filtreg_p.F : bug correction (array bounds)
+- guide_p_mod.F90 : added the "use pres2lev_mod"
+- conf_gcm.F : cosmetics (and evaluation of vert_prof_dissip , and
+               also of dissip_* factors, for Earth model)
+               plus check if "adjust" is indeed not used in OpenMP
+- comconst.h : add dissip_factz,dissip_zref variables (for Earth mode dissip)
Index: trunk/DOC/chantiers/meschantiers-Ehouarn.txt
===================================================================
--- trunk/DOC/chantiers/meschantiers-Ehouarn.txt	(revision 978)
+++ trunk/DOC/chantiers/meschantiers-Ehouarn.txt	(revision 979)
@@ -14,10 +14,10 @@
   
 - dernière mise à jour des dynamiques (par rapport au trunk de LMDZ5 
-  "Terrestre" rev. 1678) avec la rev 841.
+  "Terrestre" rev. 1760) avec la rev 979.
   Reste à faire l'exercice de validation des modifs en //...
   Avec rev 841 et 843 Cp(T) en // OK
   Reste toutefois à paufiner les choses au niveau de calfis_p et des
   domaines d'utilisation de tpot2t_p et t2tpot_p (toutes les taches MPI
-  font l'intégralité du travail!)
+  font l'intégralité du travail!). OK.
 
 - Faire qu'on puisse compiler un exécutable qui serait dans la physique
Index: trunk/LMDZ.COMMON/build_gcm
===================================================================
--- trunk/LMDZ.COMMON/build_gcm	(revision 978)
+++ trunk/LMDZ.COMMON/build_gcm	(revision 979)
@@ -1,19 +1,48 @@
 #!/bin/bash
+# Usage:
+# ./build_gcm [path]
+# where 'path' is an optional argument to this script specifying
+# the full path to where the 'fcm' command is to be found
 
-if [[ -f '.lock' ]] ; then
-  echo 'ATTENTION: vous etes sans doute en train de compiler le modele par ailleurs'
-  echo "Attendez que la premiere compilation soit terminee pour relancer la suivante."
-  echo "Si vous etes sur que vous ne compilez pas le modele par ailleurs,"
-  echo  vous pouvez continuer en repondant oui.
-  echo "Voulez-vous vraiment continuer?"
-  read -p "" answer
-  if [[ $answer != "oui" ]] ; then
+if test -f '.lock' 
+  then
+    echo 'ATTENTION: vous etes sans doute en train de compiler le modele par ailleurs'
+    echo "Attendez que la premiere compilation soit terminee pour relancer la suivante."
+    echo "Si vous etes sur que vous ne compilez pas le modele par ailleurs,"
+    echo "vous pouvez continuer en repondant oui."
+    echo "Voulez-vous vraiment continuer?"
+    echo ""
+    echo "WARNING: you are probably already compiling the model somewhere else."
+    echo "Wait until the first compilation is finished before launching this one."
+    echo "If you are sure that you are not compiling elsewhere, just answer "
+    echo "yes (or 'oui') to the question below to proceed."
+    echo "Do you wish to continue?"
+    read ouinon
+    if [[ $ouinon == "oui" || $ouinon == "yes" ]] 
+	then
+	echo OK 
+    else
+	exit
+    fi
+else
+    echo "compiling..." > '.lock'
+fi
+
+dirname="" #path to where the fcm command will be found
+if (( $# >= 1 )) ; then
+  dirname=$1
+  # check that "dirname" exists and is a directory
+  if [[ ! -d $dirname ]] ; then
+    echo "$0 error : $dirname is not a directory"
     exit
   fi
-else
-  echo "compilation en cours..." > '.lock'
+  # add a trailing "/" to $dirname
+  dirname=${dirname}"/"
 fi
 
-fcm build
+# run "fcm build" command
+${dirname}fcm build
 
+# cleanup
 \rm -f '.lock' 
+
Index: trunk/LMDZ.COMMON/build_gcm.csh
===================================================================
--- trunk/LMDZ.COMMON/build_gcm.csh	(revision 978)
+++ 	(revision )
@@ -1,24 +1,0 @@
-#!/bin/csh
-
-if ( -f '.lock' ) then
-  echo 'ATTENTION: vous etes sans doute en train de compiler le modele par ailleurs'
-  echo "Attendez que la premiere compilation soit terminee pour relancer la suivante."
-  echo "Si vous etes sur que vous ne compilez pas le modele par ailleurs,"
-  echo  vous pouvez continuer en repondant oui.
-  echo "Voulez-vous vraiment continuer?"
-
-  if ( $< == "oui" ) then
-    
-  else
-    exit
-  endif
-else
-  echo "compilation en cours..." > '.lock'
-endif
-
-#set arch=$1
-
-
-fcm build
-
-\rm -f '.lock' 
Index: trunk/LMDZ.COMMON/create_make_gcm
===================================================================
--- trunk/LMDZ.COMMON/create_make_gcm	(revision 978)
+++ trunk/LMDZ.COMMON/create_make_gcm	(revision 979)
@@ -4,4 +4,7 @@
 #
 #set -xv
+
+rm -f create_make_gcm.log
+
 machine=`hostname`
 os=`uname`
@@ -134,14 +137,14 @@
 do
    if [ -d $diri ] ; then
-   if [ "`ls $diri/*.F`" != "" ] || [ "`ls $diri/*.F90`" != "" ]  ; then 
-      cd $diri >/dev/null 2>&1
+   if [ "`ls $diri/*.F $diri/*/*.F`" != "" ] || [ "`ls $diri/*.F90 $diri/*/*.F90`" != "" ]  ; then 
+#      cd $diri >/dev/null 2>&1
       echo
       listlib=""
-      for i in `ls *.F`
-      do
-         fili=`basename $i .F`
+# Liste des fichiers .F et .F90 n'etant pas des programmes principaux
+      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
-         test=` (  head $i | grep -i '^ *program' ) `
+         test=` (  head $fili | grep -i '^ *program' ) `
          if [ "$test" = "" ] ; then 
             # if it is not a main program, add it to the list
@@ -149,10 +152,9 @@
          fi
       done
-      for i in `ls *.F90`
-      do
-         fili=`basename $i .F90`
+      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)
-         test=` (  head $i | grep -i '^ \{0,\}program' ) `
+         test=` (  head $fili | grep -i '^ \{0,\}program' ) `
          if [ "$test" = "" ] ; then
             # if it is not a main program, add it to the list
@@ -162,4 +164,5 @@
 #
       echo
+      echo $listlib >> create_make_gcm.log
       echo
       echo "#======================================================================="
@@ -167,6 +170,7 @@
       echo "#======================================================================="
       echo
-      for fili in $listlib
-      do
+      for trufile in $listlib
+      do
+         fili=`echo $trufile | awk -F/ ' { print $NF } ' | cut -d. -f1`
          echo '$(LIBO)/lib'$diri".a : " '$(LIBO)/lib'$diri".a("$fili".o)"
          echo
@@ -177,47 +181,53 @@
       echo "# Compilation des membres de la bibliotheque lib"$diri".a"
       echo
-      for fili in $listlib
-      do
-         if [ -f $fili.F90 ] ; then
-           trufile=$fili.F90
-         else
-           trufile=$fili.F
-         fi
+      for trufile in $listlib
+      do
+#         if [ -f $fili.F90 ] ; then
+#           trufile=$fili.F90
+#         else
+#           trufile=$fili.F
+#         fi
+         fili=`echo $trufile | awk -F/ ' { print $NF } ' | cut -d. -f1`
          F90=0 ; egrep -i '^ *use ' $trufile > /dev/null 2>&1 && F90=1
                  egrep -i '^ *module ' $trufile > /dev/null 2>&1 && F90=1
                  egrep -i '#include*.inc ' $trufile > /dev/null 2>&1 && F90=1
-         str1='$(LIBO)/lib'$diri'.a('$fili'.o) : $(LIBF)/'$diri/$trufile
+         str1='$(LIBO)/lib'$diri'.a('$fili'.o) : $(LIBF)/'$trufile
          [ "$fili" = "chem.subs" ] && str1=$str1' $(LIBF)/'$diri/chem.mods.F
-         for stri in ` ( sed -n "/\#include/s/\#include//p" $trufile | sed 's/\"//g' ; egrep -i '^ *use ' $trufile | sed -e 's/,/ /' | awk ' { print $2 } ' ) `
-         do
-
-
+
+##########################################################################
+# TRAITEMENT DES DEPENDANCES
 # Differents cas de dependance correspondant a des include ou des
 # use module.
 # soit dans le repertoire local soit dans un autre.
-
-            stri=`echo $stri | tr [A-Z] [a-z]`
-            if [ -f $stri ] ; then
+# Pour dyn3d, il faudrait remplacer la chaine de caractere par $DYN
+# ou quelque chose du genre.
+##########################################################################
+
+         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
+            strj=`echo $stri | tr [A-Z] [a-z]`
+            str2=""
+            for dirinc in filtrez bibio grid dyn3d $diri $diri/*/ ; do
+# Recherche dans l'ordre hierarchique inverse car seule la derniere
+# ligne est conservee
+               if [ -f $dirinc/$stri ] ; then
+                  str2='$(LIBF)/'$dirinc/$stri
+               elif [ -f $dirinc/$strj ] ; then
+                  str2='$(LIBF)/'$dirinc/$stri
+               elif [ -f $dirinc/$strj.F90 ]  || [ -f $dirinc/$strj.F ]  ; then 
+                  strlib=`echo $dirinc | awk -F/ ' { print $1 } '`
+                  str2='$(LIBO)/lib'$strlib'.a('$strj'.o)'
+               fi
+            done
+            if [ "$str2" = "" ] ; then
+               echo Warn: Include $diri $trufile $fili $stri pas trouve >> $gcm/create_make_gcm.log
+            elif [ "$str2" = "$str1" ] ; then
+               echo Warn: Include $diri $trufile $fili $stri $str1 trouve 2x >> $gcm/create_make_gcm.log
+            else
                echo $str1 \\
-               str1='$(LIBF)/'$diri'/'$stri
-            else
-               if [ -f $stri.F ] || [ -f $stri.F90 ] ; then
-                  echo $str1 \\
-                  str1='$(LIBO)/lib'$diri'.a('$stri'.o)'
-               else
-                  for dirinc in dyn3d grid bibio filtrez
-                  do
-                     if [ -f ../$dirinc/$stri ] ; then
-                        echo $str1 \\
-                        str1='$(LIBF)/'`cd .. ; ls */$stri | head -1`
-                     fi
-                     if [ -f ../$dirinc/$stri.F90 ] ; then
-                        echo $str1 \\
-                        str1='$(LIBO)/lib'$dirinc'.a('$stri'.o)'
-                     fi
-                  done
-               fi
+               str1=$str2
             fi
-         done
+         done # Fin de la boucle sur les dependances et ecriture de la derniere
          echo $str1
          # Compile in LIBO directory; and before compiling, remove
@@ -227,12 +237,12 @@
 	 if [ "$F90" -eq '0' ] ; then
          ## Fixed Form Fortran 77
-	   echo '	$(COMPILE) $(LIBF)/'$diri'/'$trufile' ; \'
+	   echo '	$(COMPILE) $(LIBF)/'$trufile' ; \'
 	 else
          ## Fortran 90
            if [ -f $fili.F90 ] ; then
 	      ## Free Form
-              echo '	$(COMPTRU90) $(LIBF)/'$diri'/'$trufile' ; \'
+              echo '	$(COMPTRU90) $(LIBF)/'$trufile' ; \'
            else
-	      echo '	$(COMPILE90) $(LIBF)/'$diri'/'$trufile' ; \'
+	      echo '	$(COMPILE90) $(LIBF)/'$trufile' ; \'
            fi
 	 fi
Index: trunk/LMDZ.COMMON/libf/bibio/new_unit.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/new_unit.F90	(revision 978)
+++ 	(revision )
@@ -1,25 +1,0 @@
-module new_unit_m
-
-  implicit none
-
-contains
-
-  subroutine new_unit(unit)
-
-    integer, intent(out):: unit
-
-    ! Variables local to the procedure:
-    logical opened, exist
-
-    !------------------------------------------------------
-
-    unit = 0
-    do
-       inquire(unit=unit, opened=opened, exist=exist)
-       if (exist .and. .not. opened) exit
-       unit = unit + 1
-    end do
-
-  end subroutine new_unit
-
-end module new_unit_m
Index: trunk/LMDZ.COMMON/libf/bibio/new_unit_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/new_unit_m.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/bibio/new_unit_m.F90	(revision 979)
@@ -0,0 +1,25 @@
+module new_unit_m
+
+  implicit none
+
+contains
+
+  subroutine new_unit(unit)
+
+    integer, intent(out):: unit
+
+    ! Variables local to the procedure:
+    logical opened, exist
+
+    !------------------------------------------------------
+
+    unit = 0
+    do
+       inquire(unit=unit, opened=opened, exist=exist)
+       if (exist .and. .not. opened) exit
+       unit = unit + 1
+    end do
+
+  end subroutine new_unit
+
+end module new_unit_m
Index: trunk/LMDZ.COMMON/libf/bibio/pchfe_95.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/pchfe_95.F90	(revision 978)
+++ 	(revision )
@@ -1,76 +1,0 @@
-module PCHFE_95_m
-
-  implicit none
-
-contains
-
-  SUBROUTINE PCHFE_95(X, F, D, SKIP, XE, FE, IERR)
-
-    ! PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
-    !            points.  May be used by itself for Hermite interpolation,
-    !            or as an evaluator for PCHIM or PCHIC.
-    ! CATEGORY  E3
-    ! KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
-    !             PIECEWISE CUBIC EVALUATION
-
-    !          PCHFE:  Piecewise Cubic Hermite Function Evaluator
-    ! Evaluates the cubic Hermite function defined by  X, F, D  at
-    ! the points  XE.
-
-    use assert_eq_m, only: assert_eq
-
-    REAL, intent(in):: X(:) ! real array of independent variable values
-    ! The elements of X must be strictly increasing.
-
-    REAL, intent(in):: F(:) ! real array of function values
-    ! F(I) is the value corresponding to X(I).
-
-    REAL, intent(in):: D(:) ! real array of derivative values
-    ! D(I) is the value corresponding to X(I).
-
-    LOGICAL, intent(inout):: SKIP 
-    ! request to skip checks for validity of "x"
-    ! If "skip" is false then "pchfe" will check that size(x) >= 2 and
-    ! "x" is in strictly ascending order.
-    ! Setting "skip" to true will save time in case these checks have
-    ! already been performed (say, in "PCHIM" or "PCHIC").
-    ! "SKIP" will be set to TRUE on normal return.
-
-    real, intent(in):: XE(:) ! points at which the function is to be evaluated
-    ! NOTES:
-    ! 1. The evaluation will be most efficient if the elements of XE
-    ! are increasing relative to X.
-    ! That is,   XE(J) .GE. X(I)
-    ! implies    XE(K) .GE. X(I),  all K.GE.J
-    ! 2. If any of the XE are outside the interval [X(1),X(N)], values
-    ! are extrapolated from the nearest extreme cubic, and a warning
-    ! error is returned.
-
-    real, intent(out):: FE(:) ! values of the cubic Hermite function
-    ! defined by X, F, D at the points XE
-
-    integer, intent(out):: IERR ! error flag
-    ! Normal return:
-    ! IERR = 0  no error
-    ! Warning error:
-    ! IERR > 0  means that extrapolation was performed at IERR points
-    ! "Recoverable" errors:
-    !              IERR = -1  if N < 2
-    !              IERR = -3  if the X-array is not strictly increasing
-    !              IERR = -4  if NE < 1
-    ! NOTE: The above errors are checked in the order listed, and
-    ! following arguments have **NOT** been validated.
-
-    ! Variables local to the procedure:
-
-    INTEGER  N, NE
-
-    !---------------------------------------
-
-    n = assert_eq(size(x), size(f), size(d), "PCHFE_95 n")
-    ne = assert_eq(size(xe), size(fe), "PCHFE_95 ne")
-    call PCHFE(N, X, F, D, 1, SKIP, NE, XE, FE, IERR)
-
-  end SUBROUTINE PCHFE_95
-
-end module PCHFE_95_m
Index: trunk/LMDZ.COMMON/libf/bibio/pchfe_95_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/pchfe_95_m.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/bibio/pchfe_95_m.F90	(revision 979)
@@ -0,0 +1,76 @@
+module PCHFE_95_m
+
+  implicit none
+
+contains
+
+  SUBROUTINE PCHFE_95(X, F, D, SKIP, XE, FE, IERR)
+
+    ! PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
+    !            points.  May be used by itself for Hermite interpolation,
+    !            or as an evaluator for PCHIM or PCHIC.
+    ! CATEGORY  E3
+    ! KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
+    !             PIECEWISE CUBIC EVALUATION
+
+    !          PCHFE:  Piecewise Cubic Hermite Function Evaluator
+    ! Evaluates the cubic Hermite function defined by  X, F, D  at
+    ! the points  XE.
+
+    use assert_eq_m, only: assert_eq
+
+    REAL, intent(in):: X(:) ! real array of independent variable values
+    ! The elements of X must be strictly increasing.
+
+    REAL, intent(in):: F(:) ! real array of function values
+    ! F(I) is the value corresponding to X(I).
+
+    REAL, intent(in):: D(:) ! real array of derivative values
+    ! D(I) is the value corresponding to X(I).
+
+    LOGICAL, intent(inout):: SKIP 
+    ! request to skip checks for validity of "x"
+    ! If "skip" is false then "pchfe" will check that size(x) >= 2 and
+    ! "x" is in strictly ascending order.
+    ! Setting "skip" to true will save time in case these checks have
+    ! already been performed (say, in "PCHIM" or "PCHIC").
+    ! "SKIP" will be set to TRUE on normal return.
+
+    real, intent(in):: XE(:) ! points at which the function is to be evaluated
+    ! NOTES:
+    ! 1. The evaluation will be most efficient if the elements of XE
+    ! are increasing relative to X.
+    ! That is,   XE(J) .GE. X(I)
+    ! implies    XE(K) .GE. X(I),  all K.GE.J
+    ! 2. If any of the XE are outside the interval [X(1),X(N)], values
+    ! are extrapolated from the nearest extreme cubic, and a warning
+    ! error is returned.
+
+    real, intent(out):: FE(:) ! values of the cubic Hermite function
+    ! defined by X, F, D at the points XE
+
+    integer, intent(out):: IERR ! error flag
+    ! Normal return:
+    ! IERR = 0  no error
+    ! Warning error:
+    ! IERR > 0  means that extrapolation was performed at IERR points
+    ! "Recoverable" errors:
+    !              IERR = -1  if N < 2
+    !              IERR = -3  if the X-array is not strictly increasing
+    !              IERR = -4  if NE < 1
+    ! NOTE: The above errors are checked in the order listed, and
+    ! following arguments have **NOT** been validated.
+
+    ! Variables local to the procedure:
+
+    INTEGER  N, NE
+
+    !---------------------------------------
+
+    n = assert_eq(size(x), size(f), size(d), "PCHFE_95 n")
+    ne = assert_eq(size(xe), size(fe), "PCHFE_95 ne")
+    call PCHFE(N, X, F, D, 1, SKIP, NE, XE, FE, IERR)
+
+  end SUBROUTINE PCHFE_95
+
+end module PCHFE_95_m
Index: trunk/LMDZ.COMMON/libf/bibio/pchsp_95.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/pchsp_95.F90	(revision 978)
+++ 	(revision )
@@ -1,116 +1,0 @@
-module pchsp_95_m
-
-  implicit none
-
-contains
-
-  function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)
-
-    ! PURPOSE: Set derivatives needed to determine the Hermite
-    ! representation of the cubic spline interpolant to given data,
-    ! with specified boundary conditions.
-
-    ! Part of the "pchip" package.
-
-    ! CATEGORY: E1A
-
-    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
-    ! interpolation, spline interpolation
-
-    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
-    ! Computes the Hermite representation of the cubic spline
-    ! interpolant to the data given in X and F satisfying the boundary
-    ! conditions specified by Ibeg, iend, vc_beg and VC_end.
-
-    ! The resulting piecewise cubic Hermite function may be evaluated
-    ! by "pchfe" or "pchfd".
-
-    ! NOTE: This is a modified version of C. de Boor's cubic spline
-    ! routine "cubspl".
-
-    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
-    ! 2001, pages 43-47
-
-    use assert_eq_m, only: assert_eq
-
-    real, intent(in):: x(:)
-    ! independent variable values
-    ! The elements of X must be strictly increasing:
-    !                X(I-1) < X(I),  I = 2...N.
-    !           (Error return if not.)
-    ! (error if size(x) < 2)
-
-    real, intent(in):: f(:)
-    !     dependent variable values to be interpolated
-    !  F(I) is value corresponding to X(I).
-
-    INTEGER, intent(in):: ibeg
-    !     desired boundary condition at beginning of data
-
-    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
-    !              tinuous at X(2).  This is the "not a knot" condition
-    !              provided by de Boor's cubic spline routine CUBSPL.
-    !              This is the default boundary condition.
-    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
-    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
-    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
-    !              (Reverts to the default boundary condition if size(x) < 3 .)
-    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
-    !              (Reverts to the default boundary condition if size(x) < 4 .)
-
-    !          NOTES:
-    !           1. An error return is taken if IBEG is out of range.
-    !           2. For the "natural" boundary condition, use IBEG=2 and
-    !              VC_BEG=0.
-
-    INTEGER, intent(in):: iend
-    !           IC(2) = IEND, desired condition at end of data.
-    !  IEND may take on the same values as IBEG, but applied to
-    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.
-
-    !          NOTES:
-    !           1. An error return is taken if IEND is out of range.
-    !           2. For the "natural" boundary condition, use IEND=2 and
-    !              VC_END=0.
-
-    REAL, intent(in), optional:: vc_beg
-    ! desired boundary value, as indicated above.
-    !           VC_BEG need be set only if IBEG = 1 or 2 .
-
-    REAL, intent(in), optional:: vc_end
-    ! desired boundary value, as indicated above.
-    !           VC_END need be set only if Iend = 1 or 2 .
-
-    real pchsp_95(size(x))
-    ! derivative values at the data points
-    !           These values will determine the cubic spline interpolant
-    !           with the requested boundary conditions.
-    !           The value corresponding to X(I) is stored in
-    !                PCHSP_95(I),  I=1...N.
-
-    ! LOCAL VARIABLES:
-    real wk(2, size(x)) ! real array of working storage
-    INTEGER n ! number of data points
-    integer ierr, ic(2)
-    real vc(2)
-
-    !-------------------------------------------------------------------
-
-    n = assert_eq(size(x), size(f), "pchsp_95 n")
-    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
-       print *, "vc_beg required for IBEG = 1 or 2"
-       stop 1
-    end if
-    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
-       print *, "vc_end required for IEND = 1 or 2"
-       stop 1
-    end if
-    ic = (/ibeg, iend/)
-    if (present(vc_beg)) vc(1) = vc_beg
-    if (present(vc_end)) vc(2) = vc_end
-    call PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
-    if (ierr /= 0) stop 1
-
-  END function pchsp_95
-
-end module pchsp_95_m
Index: trunk/LMDZ.COMMON/libf/bibio/pchsp_95_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/bibio/pchsp_95_m.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/bibio/pchsp_95_m.F90	(revision 979)
@@ -0,0 +1,116 @@
+module pchsp_95_m
+
+  implicit none
+
+contains
+
+  function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)
+
+    ! PURPOSE: Set derivatives needed to determine the Hermite
+    ! representation of the cubic spline interpolant to given data,
+    ! with specified boundary conditions.
+
+    ! Part of the "pchip" package.
+
+    ! CATEGORY: E1A
+
+    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
+    ! interpolation, spline interpolation
+
+    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
+    ! Computes the Hermite representation of the cubic spline
+    ! interpolant to the data given in X and F satisfying the boundary
+    ! conditions specified by Ibeg, iend, vc_beg and VC_end.
+
+    ! The resulting piecewise cubic Hermite function may be evaluated
+    ! by "pchfe" or "pchfd".
+
+    ! NOTE: This is a modified version of C. de Boor's cubic spline
+    ! routine "cubspl".
+
+    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
+    ! 2001, pages 43-47
+
+    use assert_eq_m, only: assert_eq
+
+    real, intent(in):: x(:)
+    ! independent variable values
+    ! The elements of X must be strictly increasing:
+    !                X(I-1) < X(I),  I = 2...N.
+    !           (Error return if not.)
+    ! (error if size(x) < 2)
+
+    real, intent(in):: f(:)
+    !     dependent variable values to be interpolated
+    !  F(I) is value corresponding to X(I).
+
+    INTEGER, intent(in):: ibeg
+    !     desired boundary condition at beginning of data
+
+    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
+    !              tinuous at X(2).  This is the "not a knot" condition
+    !              provided by de Boor's cubic spline routine CUBSPL.
+    !              This is the default boundary condition.
+    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
+    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
+    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
+    !              (Reverts to the default boundary condition if size(x) < 3 .)
+    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
+    !              (Reverts to the default boundary condition if size(x) < 4 .)
+
+    !          NOTES:
+    !           1. An error return is taken if IBEG is out of range.
+    !           2. For the "natural" boundary condition, use IBEG=2 and
+    !              VC_BEG=0.
+
+    INTEGER, intent(in):: iend
+    !           IC(2) = IEND, desired condition at end of data.
+    !  IEND may take on the same values as IBEG, but applied to
+    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.
+
+    !          NOTES:
+    !           1. An error return is taken if IEND is out of range.
+    !           2. For the "natural" boundary condition, use IEND=2 and
+    !              VC_END=0.
+
+    REAL, intent(in), optional:: vc_beg
+    ! desired boundary value, as indicated above.
+    !           VC_BEG need be set only if IBEG = 1 or 2 .
+
+    REAL, intent(in), optional:: vc_end
+    ! desired boundary value, as indicated above.
+    !           VC_END need be set only if Iend = 1 or 2 .
+
+    real pchsp_95(size(x))
+    ! derivative values at the data points
+    !           These values will determine the cubic spline interpolant
+    !           with the requested boundary conditions.
+    !           The value corresponding to X(I) is stored in
+    !                PCHSP_95(I),  I=1...N.
+
+    ! LOCAL VARIABLES:
+    real wk(2, size(x)) ! real array of working storage
+    INTEGER n ! number of data points
+    integer ierr, ic(2)
+    real vc(2)
+
+    !-------------------------------------------------------------------
+
+    n = assert_eq(size(x), size(f), "pchsp_95 n")
+    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
+       print *, "vc_beg required for IBEG = 1 or 2"
+       stop 1
+    end if
+    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
+       print *, "vc_end required for IEND = 1 or 2"
+       stop 1
+    end if
+    ic = (/ibeg, iend/)
+    if (present(vc_beg)) vc(1) = vc_beg
+    if (present(vc_end)) vc(2) = vc_end
+    call PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
+    if (ierr /= 0) stop 1
+
+  END function pchsp_95
+
+end module pchsp_95_m
Index: trunk/LMDZ.COMMON/libf/dyn3d/comconst.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/comconst.h	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/comconst.h	(revision 979)
@@ -29,4 +29,7 @@
       REAL g ! (m/s2) gravity
       REAL omeg ! (rad/s) rotation rate of the planet
+! Dissipation factors, for Earth model:
+      REAL dissip_factz,dissip_zref !dissip_deltaz
+! Dissipation factors, for other planets:
       REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
       REAL dissip_pupstart
Index: trunk/LMDZ.COMMON/libf/dyn3d/comdissnew.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/comdissnew.h	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/comdissnew.h	(revision 979)
@@ -12,8 +12,15 @@
 
       COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv,  &
-     &                   tetagrot,tetatemp,coefdis 
+     &                   tetagrot,tetatemp,coefdis, vert_prof_dissip
 
       LOGICAL lstardis
       INTEGER nitergdiv, nitergrot, niterh
+
+! For the Earth model:
+      integer vert_prof_dissip ! vertical profile of horizontal dissipation
+!     Allowed values:
+!     0: rational fraction, function of pressure
+!     1: tanh of altitude
+
       REAL     tetagdiv, tetagrot,  tetatemp, coefdis
 
Index: trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 979)
@@ -14,4 +14,6 @@
 #endif
       USE infotrac, ONLY : type_trac
+      use assert_m, only: assert
+
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -88,5 +90,6 @@
       CALL getin('lunout', lunout)
       IF (lunout /= 5 .and. lunout /= 6) THEN
-        OPEN(lunout,FILE='lmdz.out')
+        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write', 
+     &          STATUS='unknown',FORM='formatted')
       ENDIF
 
@@ -311,4 +314,18 @@
        CALL getin('tetatemp',tetatemp )
 
+! For Earth model only:
+! Parametres controlant la variation sur la verticale des constantes de
+! dissipation.
+! Pour le moment actifs uniquement dans la version a 39 niveaux
+! avec ok_strato=y
+
+       dissip_factz=4.
+       dissip_deltaz=10.
+       dissip_zref=30.
+       CALL getin('dissip_factz',dissip_factz )
+       CALL getin('dissip_deltaz',dissip_deltaz )
+       CALL getin('dissip_zref',dissip_zref )
+
+! For other planets:
 ! Parametres controlant la variation sur la verticale des constantes de
 ! dissipation.
@@ -474,5 +491,5 @@
         write(lunout,*)'STOP !!!'
         write(lunout,*)'use_filtre_fft n est pas implemente dans dyn3d'
-        STOP
+        STOP 1
       ENDIF
       
@@ -484,4 +501,14 @@
       ok_strato=.TRUE.
       CALL getin('ok_strato',ok_strato)
+
+! NB: vert_prof_dissip is Earth-specific; should not impact other models
+      if (planet_type=="earth") then
+       vert_prof_dissip = merge(1, 0, ok_strato .and. llm==39)
+       CALL getin('vert_prof_dissip', vert_prof_dissip)
+       call assert(vert_prof_dissip == 0 .or. vert_prof_dissip ==  1,
+     $     "bad value for vert_prof_dissip")
+      else
+       vert_prof_dissip=0
+      endif
 
 !Config  Key  = ok_gradsfile
@@ -834,4 +861,5 @@
       write(lunout,*)' day_step = ', day_step
       write(lunout,*)' iperiod = ', iperiod
+      write(lunout,*)' nsplit_phys = ', nsplit_phys
       write(lunout,*)' iconser = ', iconser
       write(lunout,*)' iecri = ', iecri
Index: trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90	(revision 978)
+++ 	(revision )
@@ -1,530 +1,0 @@
-!
-! $Id: etat0_netcdf.F90 1625 2012-05-09 13:14:48Z lguez $
-!
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
-!
-!-------------------------------------------------------------------------------
-! Purpose: Creates initial states
-!-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-  USE startvar
-  USE ioipsl
-  USE dimphy
-  USE infotrac
-  USE fonte_neige_mod
-  USE pbl_surface_mod
-  USE phys_state_var_mod
-  USE filtreg_mod
-  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
-  USE conf_phys_m,            ONLY: conf_phys
-! For parameterization of ozone chemistry:
-  use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
-  use press_coefoz_m, only: press_coefoz
-  use regr_pr_o3_m, only: regr_pr_o3
-  USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
-#endif
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-  LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
-  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
-  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
-  LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
-#ifndef CPP_EARTH
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "comgeom2.h"
-#include "comvert.h"
-#include "comconst.h"
-#include "indicesol.h"
-#include "dimsoil.h"
-#include "temps.h"
-  REAL,    DIMENSION(klon)                 :: tsol, qsol
-  REAL,    DIMENSION(klon)                 :: sn, rugmer, run_off_lic_0
-  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol
-  REAL,    DIMENSION(iip1,jjp1,llm+1)      :: p3d
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: uvent, t3d, tpot, qsat, qd
-  REAL,    DIMENSION(iip1,jjm ,llm)        :: vvent
-  REAL,    DIMENSION(:,:,:,:), ALLOCATABLE :: q3d
-  REAL,    DIMENSION(klon,nbsrf)           :: qsolsrf, snsrf, evap
-  REAL,    DIMENSION(klon,nbsrf)           :: frugs, agesno
-  REAL,    DIMENSION(klon,nsoilmx,nbsrf)   :: tsoil
-
-!--- Local variables for sea-ice reading:
-  INTEGER                                  :: iml_lic, jml_lic, llm_tmp
-  INTEGER                                  :: ttm_tmp, iret, fid
-  INTEGER, DIMENSION(1)                    :: itaul
-  REAL,    DIMENSION(1)                    :: lev
-  REAL                                     :: date
-  REAL,    DIMENSION(:,:),   ALLOCATABLE   ::  lon_lic,  lat_lic, fraclic
-  REAL,    DIMENSION(:),     ALLOCATABLE   :: dlon_lic, dlat_lic
-  REAL,    DIMENSION(iip1,jjp1)            :: flic_tmp
-
-!--- Misc
-  CHARACTER(LEN=80)                        :: x, fmt
-  INTEGER                                  :: i, j, l, ji
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: alpha, beta, pk, pls, y
-  REAL,    DIMENSION(ip1jmp1)              :: pks
-
-#include "comdissnew.h"
-#include "serre.h"
-#include "clesphys.h"
-
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: masse
-  INTEGER :: itau, iday
-  REAL    :: xpn, xps, time, phystep
-  REAL,    DIMENSION(iim)                  :: xppn, xpps
-  REAL,    DIMENSION(ip1jmp1,llm)          :: pbaru, phi, w
-  REAL,    DIMENSION(ip1jm  ,llm)          :: pbarv
-  REAL,    DIMENSION(klon)                 :: fder
-
-!--- Local variables for ocean mask reading:
-  INTEGER :: nid_o2a, iml_omask, jml_omask
-  LOGICAL :: couple=.FALSE.
-  REAL,    DIMENSION(:,:), ALLOCATABLE ::  lon_omask, lat_omask, ocemask, ocetmp
-  REAL,    DIMENSION(:),   ALLOCATABLE :: dlon_omask,dlat_omask
-  REAL,    DIMENSION(klon)             :: ocemask_fi
-  INTEGER, DIMENSION(klon-2)           :: isst
-  REAL,    DIMENSION(iim,jjp1)         :: zx_tmp_2d
-  REAL    :: dummy
-  LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
-  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
-  INTEGER :: iflag_radia, flag_aerosol
-  REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
-  REAL    :: tau_ratqs
-  INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake
-  INTEGER :: iflag_thermals, nsplit_thermals
-  INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
-  REAL    :: tau_thermals, solarlong0,  seuil_inversion
-  INTEGER :: read_climoz ! read ozone climatology
-!  Allowed values are 0, 1 and 2
-!     0: do not read an ozone climatology
-!     1: read a single ozone climatology that will be used day and night
-!     2: read two ozone climatologies, the average day and night
-!     climatology and the daylight climatology
-!-------------------------------------------------------------------------------
-  REAL    :: alp_offset
-  logical found
-
-!--- Constants
-  pi     = 4. * ATAN(1.)
-  rad    = 6371229.
-  daysec = 86400.
-  omeg   = 2.*pi/daysec
-  g      = 9.8
-  kappa  = 0.2857143
-  cpp    = 1004.70885
-  preff  = 101325.
-  pa     = 50000.
-  jmp1   = jjm + 1
-
-!--- CONSTRUCT A GRID
-  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
-                   callstats,                                           &
-                   solarlong0,seuil_inversion,                          &
-                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
-                   iflag_cldcon,                                        &
-                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, aerosol_couple,                      &
-                   flag_aerosol, new_aod,                               &
-                   bl95_b0, bl95_b1,                                    &
-                   read_climoz,                                         &
-                   alp_offset)
-
-! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
-  co2_ppm0 = co2_ppm
-
-  dtvr   = daysec/FLOAT(day_step)
-  WRITE(lunout,*)'dtvr',dtvr
-
-  CALL iniconst()
-  CALL inigeom()
-
-!--- Initializations for tracers
-  CALL infotrac_init
-  ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
-
-  CALL inifilr()
-  CALL phys_state_var_init(read_climoz)
-
-  rlat(1) = ASIN(1.0)
-  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
-  rlat(klon) = - ASIN(1.0)
-  rlat(:)=rlat(:)*(180.0/pi)
-
-  rlon(1) = 0.0
-  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
-  rlon(klon) = 0.0
-  rlon(:)=rlon(:)*(180.0/pi)
-
-! For a coupled simulation, the ocean mask from ocean model is used to compute
-! the weights an to insure ocean fractions are the same for atmosphere and ocean
-! Otherwise, mask is created using Relief file.
-
-  WRITE(lunout,*)'Essai de lecture masque ocean'
-  iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)
-  IF(iret/=NF90_NOERR) THEN
-    WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve'
-    WRITE(lunout,*)'Run force'
-    x='masque'
-    masque(:,:)=0.0
-    CALL startget_phys2d(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, &
-   &              rlonu, rlatv, ib)
-    WRITE(lunout,*)'MASQUE construit : Masque'
-    WRITE(lunout,'(97I1)') nINT(masque)
-    CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
-    WHERE(   zmasq(:)<EPSFRA) zmasq(:)=0.
-    WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.
-  ELSE
-    WRITE(lunout,*)'ATTENTION!! fichier o2a.nc trouve'
-    WRITE(lunout,*)'Run couple'
-    couple=.true.
-    iret=NF90_CLOSE(nid_o2a)
-    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
-    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
-      WRITE(lunout,*)'Dimensions non compatibles pour masque ocean'
-      WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask
-      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
-      CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc&
-     &ean',1)
-    END IF
-    ALLOCATE(   ocemask(iml_omask,jml_omask),   ocetmp(iml_omask,jml_omask))
-    ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
-    ALLOCATE(dlon_omask(iml_omask),         dlat_omask(jml_omask))
-    CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask,&
-   &              lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
-    CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, &
-   &              1, 1, ocetmp)
-    CALL flinclo(fid)
-    dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1)
-    dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask)
-    ocemask = ocetmp
-    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
-      DO j=1,jml_omask
-        ocemask(:,j) = ocetmp(:,jml_omask-j+1)
-      END DO
-    END IF
-!
-! Ocean mask to physical grid
-!*******************************************************************************
-    WRITE(lunout,*)'ocemask '
-    WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt)
-    WRITE(lunout,fmt)int(ocemask)
-    ocemask_fi(1)=ocemask(1,1)
-    DO j=2,jjm; ocemask_fi((j-2)*iim+2:(j-1)*iim+1)=ocemask(1:iim,j); END DO
-    ocemask_fi(klon)=ocemask(1,jjp1)
-    zmasq=1.-ocemask_fi
-  END IF
-
-  CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque)
-
-  ! The startget calls need to be replaced by a call to restget to get the
-  ! values in the restart file
-  x = 'relief'; orog(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,&
- &               masque)
-
-  x = 'rugosite'; rugo(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib)
-!  WRITE(lunout,'(49I1)') INT(orog(:,:)*10)
-!  WRITE(lunout,'(49I1)') INT(rugo(:,:)*10)
-
-! Sub-surfaces initialization
-!*******************************************************************************
-  pctsrf=0.
-  x = 'psol'; psol(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib)
-!  WRITE(lunout,*) 'PSOL :', psol(10,20)
-!  WRITE(lunout,*) ap(:), bp(:)
-
-! Mid-levels pressure computation
-!*******************************************************************************
-  CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  if (pressure_exner) then
-    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
-  else
-    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
-  endif
-  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
-!  WRITE(lunout,*) 'P3D :', p3d(10,20,:)
-!  WRITE(lunout,*) 'PK:',    pk(10,20,:)
-!  WRITE(lunout,*) 'PLS :', pls(10,20,:)
-
-  x = 'surfgeo'; phis(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,phis, 0.0,jjm, rlonu,rlatv,ib)
-
-  x = 'u';    uvent(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonu,rlatu,pls,y,uvent,0.0,  &
- &                  rlonv,rlatv,ib)
-
-  x = 'v';    vvent(:,:,:) = 0.0
-  CALL startget_dyn(x, rlonv,rlatv,pls(:, :jjm, :),y(:, :jjm, :),vvent,0.0, &
- &                  rlonu,rlatu(:jjm),ib)
-
-  x = 't';    t3d(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,y,t3d,0.0,    &
- &                  rlonu,rlatv,ib)
-
-  x = 'tpot'; tpot(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,pk,tpot,0.0,  &
- &                  rlonu,rlatv,ib)
-
-  WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:))
-  WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:))
-
-! Humidity at saturation computation
-!*******************************************************************************
-  WRITE(lunout,*) 'avant q_sat'
-  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
-  WRITE(lunout,*) 'apres q_sat'
-  WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:))
-!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
-
-  x = 'q';    qd (:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,qsat,qd,0.0, rlonu,rlatv,ib)
-  q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
-
-! Parameterization of ozone chemistry:
-! Look for ozone tracer:
-  i = 1
-  DO
-    found = tname(i)=="O3" .OR. tname(i)=="o3"
-    if (found .or. i == nqtot) exit
-    i = i + 1
-  end do
-  if (found) then
-    call regr_lat_time_coefoz
-    call press_coefoz
-    call regr_pr_o3(p3d, q3d(:, :, :, i))
-!   Convert from mole fraction to mass fraction:
-    q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
-  end if
-
-!--- OZONE CLIMATOLOGY
-  IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz)
-
-  x = 'tsol'; tsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,tsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'qsol';  qsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,qsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'snow';  sn(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,sn,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'rads';  radsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'rugmer'; rugmer(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zmea';  zmea(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zstd';  zstd(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zsig';  zsig(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zgam';  zgam(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zthe';  zthe(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zpic';  zpic(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zval';  zval(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,jjm,rlonu,rlatv,ib)
-
-!  WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273)
-
-! Soil ice file reading for soil fraction and soil ice fraction
-!*******************************************************************************
-  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
-  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
-  ALLOCATE( fraclic(iml_lic,jml_lic))
-  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
- &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
-  CALL flinclo(fid)
-
-! Interpolation on model T-grid
-!*******************************************************************************
-  WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic
-! conversion if coordinates are in degrees
-  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180.
-  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
-  dlon_lic(:)=lon_lic(:,1)
-  dlat_lic(:)=lat_lic(1,:)
-  CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1,   &
- &               rlonv, rlatu, flic_tmp(1:iim,:) )
-  flic_tmp(iip1,:)=flic_tmp(1,:)
-
-!--- To the physical grid
-  CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
-
-!--- Adequation with soil/sea mask
-  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
-  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
-  pctsrf(:,is_ter)=zmasq(:)
-  DO ji=1,klon
-    IF(zmasq(ji)>EPSFRA) THEN 
-      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
-        pctsrf(ji,is_lic)=zmasq(ji)
-        pctsrf(ji,is_ter)=0.
-      ELSE
-        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
-        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
-          pctsrf(ji,is_ter)=0.
-          pctsrf(ji,is_lic)=zmasq(ji)
-        END IF 
-      END IF 
-    END IF 
-  END DO 
-
-! sub-surface ocean and sea ice (sea ice set to zero for start)
-!*******************************************************************************
-  pctsrf(:,is_oce)=(1.-zmasq(:))
-  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
-  IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
-  isst=0
-  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
-
-! It is checked that the sub-surfaces sum is equal to 1
-!*******************************************************************************
-  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
-  IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points'
-  CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d)
-!  WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt)
-!  WRITE(lunout,*)'zmasq = '
-!  WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d)
-  CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
-  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
-  WRITE(lunout,*) 'MASQUE construit : Masque'
-  WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
-
-! Intermediate computation
-!*******************************************************************************
-  CALL massdair(p3d,masse)
-  WRITE(lunout,*)' ALPHAX ',alphax
-  DO l=1,llm
-    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
-    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
-    xpn=SUM(xppn)/apoln
-    xps=SUM(xpps)/apols
-    masse(:,1   ,l)=xpn
-    masse(:,jjp1,l)=xps
-  END DO
-  q3d(iip1,:,:,:)=q3d(1,:,:,:)
-  phis(iip1,:) = phis(1,:)
-
-  IF(.NOT.letat0) RETURN
-
-! Writing
-!*******************************************************************************
-  CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp)
-  WRITE(lunout,*)'sortie inidissip'
-  itau=0
-  itau_dyn=0
-  itau_phy=0
-  iday=dayref+itau/day_step
-  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
-  IF(time>1.) THEN
-   time=time-1
-   iday=iday+1
-  END IF
-  day_ref=dayref
-  annee_ref=anneeref
-
-  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
-  WRITE(lunout,*)'sortie geopot'
-
-  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
-                phi,  w, pbaru, pbarv, time+iday-dayref)
-  WRITE(lunout,*)'sortie caldyn0'     
-  CALL dynredem0( "start.nc", dayref, phis)
-  WRITE(lunout,*)'sortie dynredem0'
-  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
-  WRITE(lunout,*)'sortie dynredem1' 
-
-! Physical initial state writting
-!*******************************************************************************
-  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
-  phystep   = dtvr * FLOAT(iphysiq)
-  radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
-  WRITE(lunout,*)'phystep =', phystep, radpas
-
-! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
-!*******************************************************************************
-  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
-  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
-  falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6
-  falb1(:,is_oce) = 0.5;  falb1(:,is_sic) = 0.6
-  falb2 = falb1
-  evap(:,:) = 0.
-  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
-  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
-  rain_fall = 0.; snow_fall = 0.
-  solsw = 165.;   sollw = -53.
-  t_ancien = 273.15
-  q_ancien = 0.
-  agesno = 0.
-  frugs(:,is_oce) = rugmer(:)
-  frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  frugs(:,is_sic) = 0.001
-  fder = 0.0
-  clwcon = 0.0
-  rnebcon = 0.0
-  ratqs = 0.0
-  run_off_lic_0 = 0.0 
-  rugoro = 0.0
-
-! Before phyredem calling, surface modules and values to be saved in startphy.nc
-! are initialized
-!*******************************************************************************
-  dummy = 1.0
-  pbl_tke(:,:,:) = 1.e-8 
-  zmax0(:) = 40.
-  f0(:) = 1.e-5
-  ema_work1(:,:) = 0.
-  ema_work2(:,:) = 0.
-  wake_deltat(:,:) = 0.
-  wake_deltaq(:,:) = 0.
-  wake_s(:) = 0.
-  wake_cstar(:) = 0.
-  wake_fip(:) = 0.
-  wake_pe = 0.
-  fm_therm = 0.
-  entr_therm = 0.
-  detr_therm = 0.
-
-  CALL fonte_neige_init(run_off_lic_0)
-  CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
-  CALL phyredem( "startphy.nc" )
-
-!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
-!  WRITE(lunout,*)'entree histclo'
-  CALL histclo()
-
-#endif 
-!#endif of #ifdef CPP_EARTH
-  RETURN
-
-END SUBROUTINE etat0_netcdf
-!
-!-------------------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/gcm.F	(revision 979)
@@ -428,5 +428,5 @@
 
       CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
-     *                tetagdiv, tetagrot , tetatemp              )
+     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90	(revision 979)
@@ -12,4 +12,5 @@
   USE Write_Field
   use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
+  use pres2lev_mod
 
   IMPLICIT NONE
Index: trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90	(revision 979)
@@ -3,5 +3,5 @@
 !
 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
-     tetagdiv,tetagrot,tetatemp             )
+     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
   !=======================================================================
   !   initialisation de la dissipation horizontale
@@ -11,5 +11,5 @@
   !   -------------
 
-  USE control_mod, only : dissip_period,iperiod
+  USE control_mod, only : dissip_period,iperiod,planet_type
 
   IMPLICIT NONE
@@ -26,4 +26,10 @@
   REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
 
+  integer, INTENT(in):: vert_prof_dissip ! for the Earth model !!
+  ! Vertical profile of horizontal dissipation
+  ! Allowed values:
+  ! 0: rational fraction, function of pressure
+  ! 1: tanh of altitude
+
 ! Local variables:
   REAL fact,zvert(llm),zz
@@ -35,4 +41,5 @@
   INTEGER l,ij,idum,ii
   REAL tetamin
+  REAL pseudoz
   REAL Pup
   character (len=80) :: abort_message
@@ -166,15 +173,37 @@
   !   variation verticale du coefficient de dissipation:
   !   --------------------------------------------------
-
+  
+  if (planet_type.eq."earth") then
+
+   if (vert_prof_dissip == 1) then
+     do l=1,llm
+        pseudoz=8.*log(preff/presnivs(l))
+        zvert(l)=1+ &
+             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
+             *(dissip_factz-1.)
+     enddo
+   else
+     DO l=1,llm
+        zvert(l)=1.
+     ENDDO
+     fact=2.
+     DO l = 1, llm
+        zz      = 1. - preff/presnivs(l)
+        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
+     ENDDO
+   endif ! of if (vert_prof_dissip == 1)
+
+  else ! other planets
+  
 ! First step: going from 1 to dissip_fac_mid (in gcm.def)
 !============
-  DO l=1,llm
+   DO l=1,llm
      zz      = 1. - preff/presnivs(l)
      zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
-  ENDDO
-
-  write(lunout,*) 'Dissipation : '
-  write(lunout,*) 'Multiplication de la dissipation en altitude :'
-  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
+   ENDDO
+
+   write(lunout,*) 'Dissipation : '
+   write(lunout,*) 'Multiplication de la dissipation en altitude :'
+   write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
 
 ! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
@@ -188,5 +217,5 @@
 !   atmosphere par celui-ci.
 
-  if (ok_strato) then
+   if (ok_strato) then
 
     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
@@ -201,5 +230,7 @@
                    dissip_pupstart,'Pa', dissip_deltaz , '(km)'
 
-  endif
+   endif ! of if (ok_strato)
+  
+  endif ! of if (planet_type.eq."earth")
 
 
Index: trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90	(revision 978)
+++ 	(revision )
@@ -1,642 +1,0 @@
-!
-! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
-!
-!-------------------------------------------------------------------------------
-! Author : L. Fairhead, 27/01/94
-!-------------------------------------------------------------------------------
-! Purpose: Boundary conditions files building for new model using climatologies.
-!          Both grids have to be regular.
-!-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
-!-------------------------------------------------------------------------------
-! Modification history:
-!  * 23/03/1994: Z. X. Li
-!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
-!  *    07/2001: P. Le Van
-!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
-!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-  USE dimphy
-  USE ioipsl,             ONLY : ioget_year_len
-  USE phys_state_var_mod, ONLY : pctsrf
-  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
-                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
-                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
-		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
-  USE inter_barxy_m, only: inter_barxy
-#endif
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-  LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
-  LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
-  LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
-  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
-#ifndef CPP_EARTH
-  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "logic.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "comconst.h"
-#include "indicesol.h"
-
-!--- INPUT NETCDF FILES NAMES --------------------------------------------------
-  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
-  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
-                                  famipsic='amipbc_sic_1x1.nc        ',        &
-                                  fcpldsst='cpl_atm_sst.nc           ',        &
-                                  fcpldsic='cpl_atm_sic.nc           ',        &
-                                  fhistsst='histmth_sst.nc           ',        &
-                                  fhistsic='histmth_sic.nc           ',        &
-                                  frugo   ='Rugos.nc                 ',        &
-                                  falbe   ='Albedo.nc                '
-  CHARACTER(LEN=10) :: varname
-!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
-  REAL,   DIMENSION(klon)                :: fi_ice, verif
-  REAL,   DIMENSION(:,:),   POINTER      :: phy_rug=>NULL(), phy_ice=>NULL()
-  REAL,   DIMENSION(:,:),   POINTER      :: phy_sst=>NULL(), phy_alb=>NULL()
-  REAL,   DIMENSION(:,:),   ALLOCATABLE  :: phy_bil
-  REAL,   DIMENSION(:,:,:), ALLOCATABLE  :: pctsrf_t
-  INTEGER                                :: nbad
-
-!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
-  INTEGER :: ierr, nid, ndim, ntim, k
-  INTEGER, DIMENSION(2) :: dims
-  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
-  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
-  INTEGER :: NF90_FORMAT
-  INTEGER :: ndays                   !--- Depending on the output calendar
-
-!--- INITIALIZATIONS -----------------------------------------------------------
-#ifdef NC_DOUBLE
-  NF90_FORMAT=NF90_DOUBLE
-#else
-  NF90_FORMAT=NF90_FLOAT
-#endif
-
-  pi    = 4.*ATAN(1.)
-  rad   = 6371229.
-  daysec= 86400.
-  omeg  = 2.*pi/daysec
-  g     = 9.8
-  kappa = 0.2857143
-  cpp   = 1004.70885
-  dtvr  = daysec/REAL(day_step)
-  CALL inigeom
-
-!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
-  ndays=ioget_year_len(anneeref)
-
-!--- RUGOSITY TREATMENT --------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
-  varname='RUGOS'
-  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
-
-!--- OCEAN TREATMENT -----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
-
-! Input SIC file selection
-! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(famipsic)
-     varname='sicbcs'
-  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fcpldsic)
-     varname='SIICECOV'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fhistsic)
-     varname='pourc_sic'
-  ELSE
-     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
-     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
-  END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
-
-  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
-
-  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
-  DO k=1,ndays
-     fi_ice=phy_ice(:,k)
-     WHERE(fi_ice>=1.0  ) fi_ice=1.0
-     WHERE(fi_ice<EPSFRA) fi_ice=0.0
-     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
-     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
-        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
-        pctsrf_t(:,is_sic,k)=fi_ice(:)
-     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
-        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-     END IF
-     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
-     WHERE(1.0-zmasq<EPSFRA)
-        pctsrf_t(:,is_sic,k)=0.0
-        pctsrf_t(:,is_oce,k)=0.0
-     ELSEWHERE
-        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
-           pctsrf_t(:,is_sic,k)=1.0-zmasq
-           pctsrf_t(:,is_oce,k)=0.0
-        ELSEWHERE
-           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
-           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
-              pctsrf_t(:,is_oce,k)=0.0
-              pctsrf_t(:,is_sic,k)=1.0-zmasq
-           END WHERE
-        END WHERE
-     END WHERE
-     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
-     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
-     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
-  END DO
-  DEALLOCATE(phy_ice)
-
-!--- SST TREATMENT -------------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
-
-! Input SST file selection
-! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(famipsst)
-     varname='tosbcs'
-  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fcpldsst)
-     varname='SISUTESW'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fhistsst)
-     varname='tsol_oce'
-  ELSE
-     WRITE(lunout,*) 'ERROR! No sst input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
-     CALL abort_gcm('limit_netcdf','No sst file was found',1)
-  END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
-
-  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
-
-!--- ALBEDO TREATMENT ----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
-  varname='ALBEDO'
-  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
-
-!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
-  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
-
-!--- OUTPUT FILE WRITING -------------------------------------------------------
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
-
-  !--- File creation
-  ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
-  ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
-
-  !--- Dimensions creation
-  ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
-  ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
-
-  dims=(/ndim,ntim/)
-
-  !--- Variables creation
-  ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,(/ntim/),id_tim)
-  ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)
-  ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)
-  ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)
-  ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)
-  ierr=NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST)
-  ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)
-  ierr=NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB)
-  ierr=NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG)
-
-  !--- Attributes creation
-  ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")
-  ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")
-  ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")
-  ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")
-  ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")
-  ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")
-  ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")
-  ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")
-  ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")
-
-  ierr=NF90_ENDDEF(nid)
-
-  !--- Variables saving
-  ierr=NF90_PUT_VAR(nid,id_tim,(/(REAL(k),k=1,ndays)/))
-  ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
-
-  ierr=NF90_CLOSE(nid)
-
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
-
-  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
-
-
-!===============================================================================
-!
-  CONTAINS
-!
-!===============================================================================
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
-!
-!-----------------------------------------------------------------------------
-! Comments:
-!   There are two assumptions concerning the NetCDF files, that are satisfied
-!   with files that are conforming NC convention:
-!     1) The last dimension of the variables used is the time record.
-!     2) Dimensional variables have the same names as corresponding dimensions.
-!-----------------------------------------------------------------------------
-  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
-       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
-       NF90_GET_ATT
-  USE dimphy, ONLY : klon
-  USE phys_state_var_mod, ONLY : pctsrf
-  USE control_mod
-  use pchsp_95_m, only: pchsp_95
-  use pchfe_95_m, only: pchfe_95
-  use arth_m, only: arth
-
-  IMPLICIT NONE
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-#include "indicesol.h"
-#include "iniprint.h"
-!-----------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
-  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
-  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
-  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
-  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
-  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
-  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
-!------------------------------------------------------------------------------
-! Local variables:
-!--- NetCDF
-  INTEGER :: ncid, varid                  ! NetCDF identifiers
-  CHARACTER(LEN=30)               :: dnam       ! dimension name
-!--- dimensions
-  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
-  REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini   ! initial longitudes vector
-  REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini   ! initial latitudes  vector
-  REAL, POINTER,     DIMENSION(:) :: dlon, dlat ! reordered lon/lat  vectors
-!--- fields
-  INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
-  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
-  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
-  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
-  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
-  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
-!--- input files
-  CHARACTER(LEN=20)                 :: cal_in   ! calendar
-  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
-  INTEGER                           :: ndays_in ! number of days
-!--- misc
-  INTEGER :: i, j, k, l                         ! loop counters
-  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
-  CHARACTER(LEN=25)                 :: title    ! for messages
-  LOGICAL                           :: extrp    ! flag for extrapolation
-  LOGICAL                           :: oldice   ! flag for old way ice computation 
-  REAL                              :: chmin, chmax
-  INTEGER ierr
-  integer n_extrap ! number of extrapolated points
-  logical skip
-
-!------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' -------------------------------------
-  NULLIFY(champo)
-
-  SELECT CASE(mode)
-  CASE('RUG'); title='Rugosite'
-  CASE('SIC'); title='Sea-ice'
-  CASE('SST'); title='SST'
-  CASE('ALB'); title='Albedo'
-  END SELECT
-  
-
-  extrp=.FALSE. 
-  oldice=.FALSE.
-  IF ( PRESENT(flag) ) THEN 
-    IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
-    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 
-  END IF
-
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
-  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
-  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
-  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
-
-!--- Read unit for sea-ice variable only
-  IF (mode=='SIC') THEN
-     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
-     IF(ierr/=NF90_NOERR) THEN
-        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
-        unit_sic='X'
-     ELSE
-        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
-     END IF
-  END IF
-
-!--- Longitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
-
-!--- Latitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
-
-!--- Time (variable is not needed - it is rebuilt - but calendar is)
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  cal_in=' '
-  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
-  IF(ierr/=NF90_NOERR) THEN
-    SELECT CASE(mode)
-      CASE('RUG', 'ALB'); cal_in='360d'
-      CASE('SIC', 'SST'); cal_in='gregorian'
-    END SELECT
-    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
-         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
-  END IF
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
-       cal_in
-
-  
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
-  !--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref, cal_in)
-
-!--- Time vector reconstruction (time vector from file is not trusted)
-!--- If input records are not monthly, time sampling has to be constant !
-  timeyear=mid_months(anneeref, cal_in, lmdep)
-  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
-       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
-
-!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
-  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
-  IF(extrp) ALLOCATE(work(imdep, jmdep))
-
-  IF (prt_level>5) WRITE(lunout, *)
-  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
-  DO l=1, lmdep
-    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
-    CALL ncerr(ierr, fnam)
-    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
-         champ, ibar)
-
-    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
-         work)
-
-    IF(ibar .AND. .NOT.oldice) THEN
-      IF(l==1 .AND. prt_level>5) THEN
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
-        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
-      END IF
-      IF(mode=='RUG') champ=LOG(champ)
-      CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim),     &
-                         rlatv, champint)
-      IF(mode=='RUG') THEN
-        champint=EXP(champint)
-        WHERE(NINT(mask)/=1) champint=0.001
-      END IF
-    ELSE
-      SELECT CASE(mode)
-        CASE('RUG');       CALL rugosite(imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint, mask)
-        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint)
-        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint)
-      END SELECT
-    END IF
-    champtime(:, :, l)=champint
-  END DO
-  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
-
-  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
-  IF(extrp) DEALLOCATE(work)
-
-!--- TIME INTERPOLATION ------------------------------------------------------
-  IF (prt_level>5) THEN
-     WRITE(lunout, *)
-     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
-     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
-  END IF
-
-  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
-  skip = .false.
-  n_extrap = 0
-  DO j=1, jjp1
-    DO i=1, iim
-      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
-           vc_beg=0., vc_end=0.)
-      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
-           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
-      if (ierr < 0) stop 1
-      n_extrap = n_extrap + ierr
-    END DO
-  END DO
-  if (n_extrap /= 0) then
-     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
-  end if
-  champan(iip1, :, :)=champan(1, :, :)
-  DEALLOCATE(yder, champtime, timeyear)
-
-!--- Checking the result
-  DO j=1, jjp1
-    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
-  END DO
-
-!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
-  IF(mode=='SST') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
-    WHERE(champan<271.38) champan=271.38
-  END IF
-
-!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
-  IF(mode=='SIC') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-
-    IF (unit_sic=='1') THEN
-       ! Nothing to be done for sea-ice field is already in fraction of 1
-       ! This is the case for sea-ice in file cpl_atm_sic.nc
-       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
-    ELSE
-       ! Convert sea ice from percentage to fraction of 1
-       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 
-       champan(:, :, :)=champan(:, :, :)/100.
-    END IF
-
-    champan(iip1, :, :)=champan(1, :, :)
-    WHERE(champan>1.0) champan=1.0
-    WHERE(champan<0.0) champan=0.0
- END IF
-
-!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
-  ALLOCATE(champo(klon, ndays))
-  DO k=1, ndays
-    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
-  END DO
-  DEALLOCATE(champan)
-
-END SUBROUTINE get_2Dfield
-!
-!-------------------------------------------------------------------------------
-
-
-
-!-------------------------------------------------------------------------------
-!
-FUNCTION year_len(y,cal_in)
-!
-!-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER                       :: year_len
-  INTEGER,           INTENT(IN) :: y
-  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
-!-------------------------------------------------------------------------------
-!--- Getting the input calendar to reset at the end of the function
-  CALL ioget_calendar(cal_out)
-
-!--- Unlocking calendar and setting it to wanted one
-  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
-
-!--- Getting the number of days in this year
-  year_len=ioget_year_len(y)
-
-!--- Back to original calendar
-  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-END FUNCTION year_len
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-FUNCTION mid_months(y,cal_in,nm)
-!
-!-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,                INTENT(IN) :: y               ! year
-  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
-  INTEGER,                INTENT(IN) :: nm              ! months/year number
-  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=99)                  :: mess            ! error message
-  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
-  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
-  INTEGER                            :: m               ! months counter
-  INTEGER                            :: nd              ! number of days
-!-------------------------------------------------------------------------------
-  nd=year_len(y,cal_in)
-
-  IF(nm==12) THEN
-
-  !--- Getting the input calendar to reset at the end of the function
-    CALL ioget_calendar(cal_out)
-
-  !--- Unlocking calendar and setting it to wanted one
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
-
-  !--- Getting the length of each month
-    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
-
-  !--- Back to original calendar
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-  ELSE IF(MODULO(nd,nm)/=0) THEN
-    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
-      nm,' months/year. Months number should divide days number.'
-    CALL abort_gcm('mid_months',TRIM(mess),1)
-
-  ELSE
-    mnth=(/(m,m=1,nm,nd/nm)/)
-  END IF
-
-!--- Mid-months times
-  mid_months(1)=0.5*REAL(mnth(1))
-  DO k=2,nm
-    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
-  END DO
-
-END FUNCTION mid_months
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE ncerr(ncres,fnam)
-!
-!-------------------------------------------------------------------------------
-! Purpose: NetCDF errors handling.
-!-------------------------------------------------------------------------------
-  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,          INTENT(IN) :: ncres
-  CHARACTER(LEN=*), INTENT(IN) :: fnam
-!-------------------------------------------------------------------------------
-#include "iniprint.h"
-  IF(ncres/=NF90_NOERR) THEN
-    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
-    CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
-  END IF
-
-END SUBROUTINE ncerr
-!
-!-------------------------------------------------------------------------------
-
-#endif
-! of #ifdef CPP_EARTH
-
-END SUBROUTINE limit_netcdf
Index: trunk/LMDZ.COMMON/libf/dyn3d/pres2lev.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/pres2lev.F90	(revision 978)
+++ 	(revision )
@@ -1,74 +1,0 @@
-! $Id: pres2lev.F90 1418 2010-07-19 15:11:24Z jghattas $
-!
-!******************************************************
-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
-  
-  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
-  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
-
-  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
-
-  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    
Index: trunk/LMDZ.COMMON/libf/dyn3d/pres2lev_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/pres2lev_mod.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/dyn3d/pres2lev_mod.F90	(revision 979)
@@ -0,0 +1,80 @@
+! $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
+  
+  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
+  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
+
+  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
+
+  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/dyn3dpar/abort_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/abort_gcm.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/abort_gcm.F	(revision 979)
@@ -1,4 +1,4 @@
 !
-! $Id: abort_gcm.F 1475 2011-01-21 14:41:03Z emillour $
+! $Id: abort_gcm.F 1748 2013-04-24 14:18:40Z emillour $
 !
 c
@@ -24,5 +24,5 @@
 
       character(len=*) modname
-      integer ierr
+      integer ierr, ierror_mpi
       character(len=*) message
 
@@ -47,5 +47,11 @@
       else
         write(lunout,*) 'Houston, we have a problem ', ierr
+#ifdef CPP_MPI
+C$OMP CRITICAL (MPI_ABORT_GCM)
+        call MPI_ABORT(COMM_LMDZ, 1, ierror_mpi)
+C$OMP END CRITICAL (MPI_ABORT_GCM)
+#else
         stop 1
+#endif          
       endif
       END
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/comconst.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/comconst.h	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/comconst.h	(revision 979)
@@ -29,4 +29,7 @@
       REAL g ! (m/s2) gravity
       REAL omeg ! (rad/s) rotation rate of the planet
+! Dissipation factors, for Earth model:
+      REAL dissip_factz,dissip_zref !dissip_deltaz
+! Dissipation factors, for other planets:
       REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
       REAL dissip_pupstart
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/comdissnew.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/comdissnew.h	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/comdissnew.h	(revision 979)
@@ -12,8 +12,15 @@
 
       COMMON/comdissnew/ lstardis,nitergdiv,nitergrot,niterh,tetagdiv,  &
-     &                   tetagrot,tetatemp,coefdis 
+     &                   tetagrot,tetatemp,coefdis, vert_prof_dissip
 
       LOGICAL lstardis
       INTEGER nitergdiv, nitergrot, niterh
+
+! For the Earth model:
+      integer vert_prof_dissip ! vertical profile of horizontal dissipation
+!     Allowed values:
+!     0: rational fraction, function of pressure
+!     1: tanh of altitude
+
       REAL     tetagdiv, tetagrot,  tetatemp, coefdis
 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 979)
@@ -18,4 +18,5 @@
       USE control_mod
       USE infotrac, ONLY : type_trac
+      use assert_m, only: assert
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -54,5 +55,9 @@
       LOGICAL  fxyhypbb, ysinuss
       INTEGER i
-      
+      character(len=*),parameter :: modname="conf_gcm"
+      character (len=80) :: abort_message
+#ifdef CPP_OMP
+      integer,external :: OMP_GET_NUM_THREADS
+#endif      
 c
 c  -------------------------------------------------------------------
@@ -83,11 +88,5 @@
 c   initialisations:
 c   ----------------
-      adjust=.false.
-      call getin('adjust',adjust)
-      
-      itaumax=0
-      call getin('itaumax',itaumax);
-      if (itaumax<=0) itaumax=HUGE(itaumax)
-      
+
 !Config  Key  = lunout
 !Config  Desc = unite de fichier pour les impressions
@@ -103,4 +102,23 @@
       ENDIF
 
+      adjust=.false.
+      call getin('adjust',adjust)
+      
+#ifdef CPP_OMP
+      ! adjust=y not implemented in case of OpenMP threads...
+!$OMP PARALLEL
+      if ((OMP_GET_NUM_THREADS()>1).and.adjust) then
+        write(lunout,*)'conf_gcm: Error, adjust should be set to n'
+     &,' when running with OpenMP threads'
+        abort_message = 'Wrong value for adjust'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!$OMP END PARALLEL         
+#endif
+
+      itaumax=0
+      call getin('itaumax',itaumax);
+      if (itaumax<=0) itaumax=HUGE(itaumax)
+      
 !Config  Key  = prt_level
 !Config  Desc = niveau d'impressions de débogage
@@ -323,4 +341,18 @@
        CALL getin('tetatemp',tetatemp )
 
+! For Earth model only:
+! Parametres controlant la variation sur la verticale des constantes de
+! dissipation.
+! Pour le moment actifs uniquement dans la version a 39 niveaux
+! avec ok_strato=y
+
+       dissip_factz=4.
+       dissip_deltaz=10.
+       dissip_zref=30.
+       CALL getin('dissip_factz',dissip_factz )
+       CALL getin('dissip_deltaz',dissip_deltaz )
+       CALL getin('dissip_zref',dissip_zref )
+
+! For other planets:
 ! Parametres controlant la variation sur la verticale des constantes de
 ! dissipation.
@@ -494,5 +526,5 @@
         write(lunout,*)"Le zoom en longitude est incompatible",
      &                 " avec l'utilisation du filtre FFT ",
-     &                 "---> filtre FFT dÃ©sactivÃ© "
+     &                 "---> FFT filter not active"
        use_filtre_fft=.FALSE.
       ENDIF
@@ -525,4 +557,14 @@
       ok_strato=.TRUE.
       CALL getin('ok_strato',ok_strato)
+
+! NB: vert_prof_dissip is Earth-specific; should not impact other models
+      if (planet_type=="earth") then
+       vert_prof_dissip = merge(1, 0, ok_strato .and. llm==39)
+       CALL getin('vert_prof_dissip', vert_prof_dissip)
+       call assert(vert_prof_dissip == 0 .or. vert_prof_dissip ==  1,
+     $     "bad value for vert_prof_dissip")
+      else
+       vert_prof_dissip=0
+      endif
 
 !Config  Key  = ok_gradsfile
@@ -875,4 +917,5 @@
       write(lunout,*)' day_step = ', day_step
       write(lunout,*)' iperiod = ', iperiod
+      write(lunout,*)' nsplit_phys = ', nsplit_phys
       write(lunout,*)' iconser = ', iconser
       write(lunout,*)' iecri = ', iecri
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/etat0_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/etat0_netcdf.F90	(revision 978)
+++ 	(revision )
@@ -1,530 +1,0 @@
-!
-! $Id: etat0_netcdf.F90 1625 2012-05-09 13:14:48Z lguez $
-!
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
-!
-!-------------------------------------------------------------------------------
-! Purpose: Creates initial states
-!-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-  USE startvar
-  USE ioipsl
-  USE dimphy
-  USE infotrac
-  USE fonte_neige_mod
-  USE pbl_surface_mod
-  USE phys_state_var_mod
-  USE filtreg_mod
-  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
-  USE conf_phys_m,            ONLY: conf_phys
-! For parameterization of ozone chemistry:
-  use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
-  use press_coefoz_m, only: press_coefoz
-  use regr_pr_o3_m, only: regr_pr_o3
-  USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
-#endif
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-  LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
-  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
-  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
-  LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
-#ifndef CPP_EARTH
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "comgeom2.h"
-#include "comvert.h"
-#include "comconst.h"
-#include "indicesol.h"
-#include "dimsoil.h"
-#include "temps.h"
-  REAL,    DIMENSION(klon)                 :: tsol, qsol
-  REAL,    DIMENSION(klon)                 :: sn, rugmer, run_off_lic_0
-  REAL,    DIMENSION(iip1,jjp1)            :: orog, rugo, psol
-  REAL,    DIMENSION(iip1,jjp1,llm+1)      :: p3d
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: uvent, t3d, tpot, qsat, qd
-  REAL,    DIMENSION(iip1,jjm ,llm)        :: vvent
-  REAL,    DIMENSION(:,:,:,:), ALLOCATABLE :: q3d
-  REAL,    DIMENSION(klon,nbsrf)           :: qsolsrf, snsrf, evap
-  REAL,    DIMENSION(klon,nbsrf)           :: frugs, agesno
-  REAL,    DIMENSION(klon,nsoilmx,nbsrf)   :: tsoil
-
-!--- Local variables for sea-ice reading:
-  INTEGER                                  :: iml_lic, jml_lic, llm_tmp
-  INTEGER                                  :: ttm_tmp, iret, fid
-  INTEGER, DIMENSION(1)                    :: itaul
-  REAL,    DIMENSION(1)                    :: lev
-  REAL                                     :: date
-  REAL,    DIMENSION(:,:),   ALLOCATABLE   ::  lon_lic,  lat_lic, fraclic
-  REAL,    DIMENSION(:),     ALLOCATABLE   :: dlon_lic, dlat_lic
-  REAL,    DIMENSION(iip1,jjp1)            :: flic_tmp
-
-!--- Misc
-  CHARACTER(LEN=80)                        :: x, fmt
-  INTEGER                                  :: i, j, l, ji
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: alpha, beta, pk, pls, y
-  REAL,    DIMENSION(ip1jmp1)              :: pks
-
-#include "comdissnew.h"
-#include "serre.h"
-#include "clesphys.h"
-
-  REAL,    DIMENSION(iip1,jjp1,llm)        :: masse
-  INTEGER :: itau, iday
-  REAL    :: xpn, xps, time, phystep
-  REAL,    DIMENSION(iim)                  :: xppn, xpps
-  REAL,    DIMENSION(ip1jmp1,llm)          :: pbaru, phi, w
-  REAL,    DIMENSION(ip1jm  ,llm)          :: pbarv
-  REAL,    DIMENSION(klon)                 :: fder
-
-!--- Local variables for ocean mask reading:
-  INTEGER :: nid_o2a, iml_omask, jml_omask
-  LOGICAL :: couple=.FALSE.
-  REAL,    DIMENSION(:,:), ALLOCATABLE ::  lon_omask, lat_omask, ocemask, ocetmp
-  REAL,    DIMENSION(:),   ALLOCATABLE :: dlon_omask,dlat_omask
-  REAL,    DIMENSION(klon)             :: ocemask_fi
-  INTEGER, DIMENSION(klon-2)           :: isst
-  REAL,    DIMENSION(iim,jjp1)         :: zx_tmp_2d
-  REAL    :: dummy
-  LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
-  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
-  INTEGER :: iflag_radia, flag_aerosol
-  REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
-  REAL    :: tau_ratqs
-  INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake
-  INTEGER :: iflag_thermals, nsplit_thermals
-  INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
-  REAL    :: tau_thermals, solarlong0,  seuil_inversion
-  INTEGER :: read_climoz ! read ozone climatology
-!  Allowed values are 0, 1 and 2
-!     0: do not read an ozone climatology
-!     1: read a single ozone climatology that will be used day and night
-!     2: read two ozone climatologies, the average day and night
-!     climatology and the daylight climatology
-!-------------------------------------------------------------------------------
-  REAL    :: alp_offset
-  logical found
-
-!--- Constants
-  pi     = 4. * ATAN(1.)
-  rad    = 6371229.
-  daysec = 86400.
-  omeg   = 2.*pi/daysec
-  g      = 9.8
-  kappa  = 0.2857143
-  cpp    = 1004.70885
-  preff  = 101325.
-  pa     = 50000.
-  jmp1   = jjm + 1
-
-!--- CONSTRUCT A GRID
-  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
-                   callstats,                                           &
-                   solarlong0,seuil_inversion,                          &
-                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
-                   iflag_cldcon,                                        &
-                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, aerosol_couple,                      &
-                   flag_aerosol, new_aod,                               &
-                   bl95_b0, bl95_b1,                                    &
-                   read_climoz,                                         &
-                   alp_offset)
-
-! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
-  co2_ppm0 = co2_ppm
-
-  dtvr   = daysec/FLOAT(day_step)
-  WRITE(lunout,*)'dtvr',dtvr
-
-  CALL iniconst()
-  CALL inigeom()
-
-!--- Initializations for tracers
-  CALL infotrac_init
-  ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
-
-  CALL inifilr()
-  CALL phys_state_var_init(read_climoz)
-
-  rlat(1) = ASIN(1.0)
-  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
-  rlat(klon) = - ASIN(1.0)
-  rlat(:)=rlat(:)*(180.0/pi)
-
-  rlon(1) = 0.0
-  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
-  rlon(klon) = 0.0
-  rlon(:)=rlon(:)*(180.0/pi)
-
-! For a coupled simulation, the ocean mask from ocean model is used to compute
-! the weights an to insure ocean fractions are the same for atmosphere and ocean
-! Otherwise, mask is created using Relief file.
-
-  WRITE(lunout,*)'Essai de lecture masque ocean'
-  iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)
-  IF(iret/=NF90_NOERR) THEN
-    WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve'
-    WRITE(lunout,*)'Run force'
-    x='masque'
-    masque(:,:)=0.0
-    CALL startget_phys2d(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, &
-   &              rlonu, rlatv, ib)
-    WRITE(lunout,*)'MASQUE construit : Masque'
-    WRITE(lunout,'(97I1)') nINT(masque)
-    CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
-    WHERE(   zmasq(:)<EPSFRA) zmasq(:)=0.
-    WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.
-  ELSE
-    WRITE(lunout,*)'ATTENTION!! fichier o2a.nc trouve'
-    WRITE(lunout,*)'Run couple'
-    couple=.true.
-    iret=NF90_CLOSE(nid_o2a)
-    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
-    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
-      WRITE(lunout,*)'Dimensions non compatibles pour masque ocean'
-      WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask
-      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
-      CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc&
-     &ean',1)
-    END IF
-    ALLOCATE(   ocemask(iml_omask,jml_omask),   ocetmp(iml_omask,jml_omask))
-    ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
-    ALLOCATE(dlon_omask(iml_omask),         dlat_omask(jml_omask))
-    CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask,&
-   &              lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
-    CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, &
-   &              1, 1, ocetmp)
-    CALL flinclo(fid)
-    dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1)
-    dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask)
-    ocemask = ocetmp
-    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
-      DO j=1,jml_omask
-        ocemask(:,j) = ocetmp(:,jml_omask-j+1)
-      END DO
-    END IF
-!
-! Ocean mask to physical grid
-!*******************************************************************************
-    WRITE(lunout,*)'ocemask '
-    WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt)
-    WRITE(lunout,fmt)int(ocemask)
-    ocemask_fi(1)=ocemask(1,1)
-    DO j=2,jjm; ocemask_fi((j-2)*iim+2:(j-1)*iim+1)=ocemask(1:iim,j); END DO
-    ocemask_fi(klon)=ocemask(1,jjp1)
-    zmasq=1.-ocemask_fi
-  END IF
-
-  CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque)
-
-  ! The startget calls need to be replaced by a call to restget to get the
-  ! values in the restart file
-  x = 'relief'; orog(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,&
- &               masque)
-
-  x = 'rugosite'; rugo(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib)
-!  WRITE(lunout,'(49I1)') INT(orog(:,:)*10)
-!  WRITE(lunout,'(49I1)') INT(rugo(:,:)*10)
-
-! Sub-surfaces initialization
-!*******************************************************************************
-  pctsrf=0.
-  x = 'psol'; psol(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib)
-!  WRITE(lunout,*) 'PSOL :', psol(10,20)
-!  WRITE(lunout,*) ap(:), bp(:)
-
-! Mid-levels pressure computation
-!*******************************************************************************
-  CALL pression(ip1jmp1, ap, bp, psol, p3d)
-  if (pressure_exner) then
-    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
-  else
-    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
-  endif
-  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
-!  WRITE(lunout,*) 'P3D :', p3d(10,20,:)
-!  WRITE(lunout,*) 'PK:',    pk(10,20,:)
-!  WRITE(lunout,*) 'PLS :', pls(10,20,:)
-
-  x = 'surfgeo'; phis(:,:) = 0.0
-  CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,phis, 0.0,jjm, rlonu,rlatv,ib)
-
-  x = 'u';    uvent(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonu,rlatu,pls,y,uvent,0.0,  &
- &                  rlonv,rlatv,ib)
-
-  x = 'v';    vvent(:,:,:) = 0.0
-  CALL startget_dyn(x, rlonv,rlatv,pls(:, :jjm, :),y(:, :jjm, :),vvent,0.0, &
- &                  rlonu,rlatu(:jjm),ib)
-
-  x = 't';    t3d(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,y,t3d,0.0,    &
- &                  rlonu,rlatv,ib)
-
-  x = 'tpot'; tpot(:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,pk,tpot,0.0,  &
- &                  rlonu,rlatv,ib)
-
-  WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:))
-  WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:))
-
-! Humidity at saturation computation
-!*******************************************************************************
-  WRITE(lunout,*) 'avant q_sat'
-  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
-  WRITE(lunout,*) 'apres q_sat'
-  WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:))
-!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
-
-  x = 'q';    qd (:,:,:) = 0.0
-  CALL startget_dyn(x,rlonv,rlatu,pls,qsat,qd,0.0, rlonu,rlatv,ib)
-  q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
-
-! Parameterization of ozone chemistry:
-! Look for ozone tracer:
-  i = 1
-  DO
-    found = tname(i)=="O3" .OR. tname(i)=="o3"
-    if (found .or. i == nqtot) exit
-    i = i + 1
-  end do
-  if (found) then
-    call regr_lat_time_coefoz
-    call press_coefoz
-    call regr_pr_o3(p3d, q3d(:, :, :, i))
-!   Convert from mole fraction to mass fraction:
-    q3d(:, :, :, i) = q3d(:, :, :, i)  * 48. / 29.
-  end if
-
-!--- OZONE CLIMATOLOGY
-  IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz)
-
-  x = 'tsol'; tsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,tsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'qsol';  qsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,qsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'snow';  sn(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,sn,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'rads';  radsol(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'rugmer'; rugmer(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zmea';  zmea(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zstd';  zstd(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zsig';  zsig(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zgam';  zgam(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zthe';  zthe(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zpic';  zpic(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,jjm,rlonu,rlatv,ib)
-
-  x = 'zval';  zval(:) = 0.0
-  CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,jjm,rlonu,rlatv,ib)
-
-!  WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273)
-
-! Soil ice file reading for soil fraction and soil ice fraction
-!*******************************************************************************
-  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
-  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
-  ALLOCATE( fraclic(iml_lic,jml_lic))
-  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
- &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
-  CALL flinclo(fid)
-
-! Interpolation on model T-grid
-!*******************************************************************************
-  WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic
-! conversion if coordinates are in degrees
-  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180.
-  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
-  dlon_lic(:)=lon_lic(:,1)
-  dlat_lic(:)=lat_lic(1,:)
-  CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1,   &
- &               rlonv, rlatu, flic_tmp(1:iim,:) )
-  flic_tmp(iip1,:)=flic_tmp(1,:)
-
-!--- To the physical grid
-  CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
-
-!--- Adequation with soil/sea mask
-  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
-  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
-  pctsrf(:,is_ter)=zmasq(:)
-  DO ji=1,klon
-    IF(zmasq(ji)>EPSFRA) THEN 
-      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
-        pctsrf(ji,is_lic)=zmasq(ji)
-        pctsrf(ji,is_ter)=0.
-      ELSE
-        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
-        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
-          pctsrf(ji,is_ter)=0.
-          pctsrf(ji,is_lic)=zmasq(ji)
-        END IF 
-      END IF 
-    END IF 
-  END DO 
-
-! sub-surface ocean and sea ice (sea ice set to zero for start)
-!*******************************************************************************
-  pctsrf(:,is_oce)=(1.-zmasq(:))
-  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
-  IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
-  isst=0
-  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
-
-! It is checked that the sub-surfaces sum is equal to 1
-!*******************************************************************************
-  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
-  IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points'
-  CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d)
-!  WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt)
-!  WRITE(lunout,*)'zmasq = '
-!  WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d)
-  CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
-  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
-  WRITE(lunout,*) 'MASQUE construit : Masque'
-  WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
-
-! Intermediate computation
-!*******************************************************************************
-  CALL massdair(p3d,masse)
-  WRITE(lunout,*)' ALPHAX ',alphax
-  DO l=1,llm
-    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
-    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
-    xpn=SUM(xppn)/apoln
-    xps=SUM(xpps)/apols
-    masse(:,1   ,l)=xpn
-    masse(:,jjp1,l)=xps
-  END DO
-  q3d(iip1,:,:,:)=q3d(1,:,:,:)
-  phis(iip1,:) = phis(1,:)
-
-  IF(.NOT.letat0) RETURN
-
-! Writing
-!*******************************************************************************
-  CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp)
-  WRITE(lunout,*)'sortie inidissip'
-  itau=0
-  itau_dyn=0
-  itau_phy=0
-  iday=dayref+itau/day_step
-  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
-  IF(time>1.) THEN
-   time=time-1
-   iday=iday+1
-  END IF
-  day_ref=dayref
-  annee_ref=anneeref
-
-  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
-  WRITE(lunout,*)'sortie geopot'
-
-  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
-                phi,  w, pbaru, pbarv, time+iday-dayref)
-  WRITE(lunout,*)'sortie caldyn0'     
-  CALL dynredem0( "start.nc", dayref, phis)
-  WRITE(lunout,*)'sortie dynredem0'
-  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
-  WRITE(lunout,*)'sortie dynredem1' 
-
-! Physical initial state writting
-!*******************************************************************************
-  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
-  phystep   = dtvr * FLOAT(iphysiq)
-  radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
-  WRITE(lunout,*)'phystep =', phystep, radpas
-
-! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
-!*******************************************************************************
-  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
-  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
-  falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6
-  falb1(:,is_oce) = 0.5;  falb1(:,is_sic) = 0.6
-  falb2 = falb1
-  evap(:,:) = 0.
-  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
-  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
-  rain_fall = 0.; snow_fall = 0.
-  solsw = 165.;   sollw = -53.
-  t_ancien = 273.15
-  q_ancien = 0.
-  agesno = 0.
-  frugs(:,is_oce) = rugmer(:)
-  frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  frugs(:,is_sic) = 0.001
-  fder = 0.0
-  clwcon = 0.0
-  rnebcon = 0.0
-  ratqs = 0.0
-  run_off_lic_0 = 0.0 
-  rugoro = 0.0
-
-! Before phyredem calling, surface modules and values to be saved in startphy.nc
-! are initialized
-!*******************************************************************************
-  dummy = 1.0
-  pbl_tke(:,:,:) = 1.e-8 
-  zmax0(:) = 40.
-  f0(:) = 1.e-5
-  ema_work1(:,:) = 0.
-  ema_work2(:,:) = 0.
-  wake_deltat(:,:) = 0.
-  wake_deltaq(:,:) = 0.
-  wake_s(:) = 0.
-  wake_cstar(:) = 0.
-  wake_fip(:) = 0.
-  wake_pe = 0.
-  fm_therm = 0.
-  entr_therm = 0.
-  detr_therm = 0.
-
-  CALL fonte_neige_init(run_off_lic_0)
-  CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
-  CALL phyredem( "startphy.nc" )
-
-!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
-!  WRITE(lunout,*)'entree histclo'
-  CALL histclo()
-
-#endif 
-!#endif of #ifdef CPP_EARTH
-  RETURN
-
-END SUBROUTINE etat0_netcdf
-!
-!-------------------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 979)
@@ -443,5 +443,5 @@
 
       CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
-     *                tetagdiv, tetagrot , tetatemp              )
+     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 979)
@@ -11,5 +11,6 @@
   USE getparam
   USE Write_Field_p
-  use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
+  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
+  USE pres2lev_mod
 
   IMPLICIT NONE
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90	(revision 979)
@@ -3,5 +3,5 @@
 !
 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
-     tetagdiv,tetagrot,tetatemp             )
+     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
   !=======================================================================
   !   initialisation de la dissipation horizontale
@@ -11,5 +11,5 @@
   !   -------------
 
-  USE control_mod, only : dissip_period,iperiod
+  USE control_mod, only : dissip_period,iperiod,planet_type
 
   IMPLICIT NONE
@@ -26,4 +26,10 @@
   REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
 
+  integer, INTENT(in):: vert_prof_dissip ! for the Earth model !!
+  ! Vertical profile of horizontal dissipation
+  ! Allowed values:
+  ! 0: rational fraction, function of pressure
+  ! 1: tanh of altitude
+
 ! Local variables:
   REAL fact,zvert(llm),zz
@@ -35,4 +41,5 @@
   INTEGER l,ij,idum,ii
   REAL tetamin
+  REAL pseudoz
   REAL Pup
   character (len=80) :: abort_message
@@ -166,15 +173,37 @@
   !   variation verticale du coefficient de dissipation:
   !   --------------------------------------------------
-
+  
+  if (planet_type.eq."earth") then
+
+   if (vert_prof_dissip == 1) then
+     do l=1,llm
+        pseudoz=8.*log(preff/presnivs(l))
+        zvert(l)=1+ &
+             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
+             *(dissip_factz-1.)
+     enddo
+   else
+     DO l=1,llm
+        zvert(l)=1.
+     ENDDO
+     fact=2.
+     DO l = 1, llm
+        zz      = 1. - preff/presnivs(l)
+        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
+     ENDDO
+   endif ! of if (vert_prof_dissip == 1)
+
+  else ! other planets
+  
 ! First step: going from 1 to dissip_fac_mid (in gcm.def)
 !============
-  DO l=1,llm
+   DO l=1,llm
      zz      = 1. - preff/presnivs(l)
      zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
-  ENDDO
-
-  write(lunout,*) 'Dissipation : '
-  write(lunout,*) 'Multiplication de la dissipation en altitude :'
-  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
+   ENDDO
+
+   write(lunout,*) 'Dissipation : '
+   write(lunout,*) 'Multiplication de la dissipation en altitude :'
+   write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
 
 ! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
@@ -188,5 +217,5 @@
 !   atmosphere par celui-ci.
 
-  if (ok_strato) then
+   if (ok_strato) then
 
     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
@@ -201,5 +230,7 @@
                    dissip_pupstart,'Pa', dissip_deltaz , '(km)'
 
-  endif
+   endif ! of if (ok_strato)
+  
+  endif ! of if (planet_type.eq."earth")
 
 
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90	(revision 978)
+++ 	(revision )
@@ -1,642 +1,0 @@
-!
-! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
-!
-!-------------------------------------------------------------------------------
-! Author : L. Fairhead, 27/01/94
-!-------------------------------------------------------------------------------
-! Purpose: Boundary conditions files building for new model using climatologies.
-!          Both grids have to be regular.
-!-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
-!-------------------------------------------------------------------------------
-! Modification history:
-!  * 23/03/1994: Z. X. Li
-!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
-!  *    07/2001: P. Le Van
-!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
-!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-  USE dimphy
-  USE ioipsl,             ONLY : ioget_year_len
-  USE phys_state_var_mod, ONLY : pctsrf
-  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
-                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
-                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
-		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
-  USE inter_barxy_m, only: inter_barxy
-#endif
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-  LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
-  LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
-  LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
-  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
-#ifndef CPP_EARTH
-  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "logic.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "comconst.h"
-#include "indicesol.h"
-
-!--- INPUT NETCDF FILES NAMES --------------------------------------------------
-  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
-  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
-                                  famipsic='amipbc_sic_1x1.nc        ',        &
-                                  fcpldsst='cpl_atm_sst.nc           ',        &
-                                  fcpldsic='cpl_atm_sic.nc           ',        &
-                                  fhistsst='histmth_sst.nc           ',        &
-                                  fhistsic='histmth_sic.nc           ',        &
-                                  frugo   ='Rugos.nc                 ',        &
-                                  falbe   ='Albedo.nc                '
-  CHARACTER(LEN=10) :: varname
-!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
-  REAL,   DIMENSION(klon)                :: fi_ice, verif
-  REAL,   DIMENSION(:,:),   POINTER      :: phy_rug=>NULL(), phy_ice=>NULL()
-  REAL,   DIMENSION(:,:),   POINTER      :: phy_sst=>NULL(), phy_alb=>NULL()
-  REAL,   DIMENSION(:,:),   ALLOCATABLE  :: phy_bil
-  REAL,   DIMENSION(:,:,:), ALLOCATABLE  :: pctsrf_t
-  INTEGER                                :: nbad
-
-!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
-  INTEGER :: ierr, nid, ndim, ntim, k
-  INTEGER, DIMENSION(2) :: dims
-  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
-  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
-  INTEGER :: NF90_FORMAT
-  INTEGER :: ndays                   !--- Depending on the output calendar
-
-!--- INITIALIZATIONS -----------------------------------------------------------
-#ifdef NC_DOUBLE
-  NF90_FORMAT=NF90_DOUBLE
-#else
-  NF90_FORMAT=NF90_FLOAT
-#endif
-
-  pi    = 4.*ATAN(1.)
-  rad   = 6371229.
-  daysec= 86400.
-  omeg  = 2.*pi/daysec
-  g     = 9.8
-  kappa = 0.2857143
-  cpp   = 1004.70885
-  dtvr  = daysec/REAL(day_step)
-  CALL inigeom
-
-!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
-  ndays=ioget_year_len(anneeref)
-
-!--- RUGOSITY TREATMENT --------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
-  varname='RUGOS'
-  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
-
-!--- OCEAN TREATMENT -----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
-
-! Input SIC file selection
-! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(famipsic)
-     varname='sicbcs'
-  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fcpldsic)
-     varname='SIICECOV'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fhistsic)
-     varname='pourc_sic'
-  ELSE
-     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
-     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
-  END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
-
-  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
-
-  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
-  DO k=1,ndays
-     fi_ice=phy_ice(:,k)
-     WHERE(fi_ice>=1.0  ) fi_ice=1.0
-     WHERE(fi_ice<EPSFRA) fi_ice=0.0
-     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
-     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
-        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
-        pctsrf_t(:,is_sic,k)=fi_ice(:)
-     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
-        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-     END IF
-     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
-     WHERE(1.0-zmasq<EPSFRA)
-        pctsrf_t(:,is_sic,k)=0.0
-        pctsrf_t(:,is_oce,k)=0.0
-     ELSEWHERE
-        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
-           pctsrf_t(:,is_sic,k)=1.0-zmasq
-           pctsrf_t(:,is_oce,k)=0.0
-        ELSEWHERE
-           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
-           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
-              pctsrf_t(:,is_oce,k)=0.0
-              pctsrf_t(:,is_sic,k)=1.0-zmasq
-           END WHERE
-        END WHERE
-     END WHERE
-     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
-     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
-     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
-     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
-  END DO
-  DEALLOCATE(phy_ice)
-
-!--- SST TREATMENT -------------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
-
-! Input SST file selection
-! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(famipsst)
-     varname='tosbcs'
-  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fcpldsst)
-     varname='SISUTESW'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fhistsst)
-     varname='tsol_oce'
-  ELSE
-     WRITE(lunout,*) 'ERROR! No sst input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
-     CALL abort_gcm('limit_netcdf','No sst file was found',1)
-  END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
-
-  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
-
-!--- ALBEDO TREATMENT ----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
-  varname='ALBEDO'
-  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
-
-!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
-  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
-
-!--- OUTPUT FILE WRITING -------------------------------------------------------
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
-
-  !--- File creation
-  ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
-  ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
-
-  !--- Dimensions creation
-  ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
-  ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
-
-  dims=(/ndim,ntim/)
-
-  !--- Variables creation
-  ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,(/ntim/),id_tim)
-  ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)
-  ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)
-  ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)
-  ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)
-  ierr=NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST)
-  ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)
-  ierr=NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB)
-  ierr=NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG)
-
-  !--- Attributes creation
-  ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")
-  ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")
-  ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")
-  ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")
-  ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")
-  ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")
-  ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")
-  ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")
-  ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")
-
-  ierr=NF90_ENDDEF(nid)
-
-  !--- Variables saving
-  ierr=NF90_PUT_VAR(nid,id_tim,(/(REAL(k),k=1,ndays)/))
-  ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
-  ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
-
-  ierr=NF90_CLOSE(nid)
-
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
-
-  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
-
-
-!===============================================================================
-!
-  CONTAINS
-!
-!===============================================================================
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
-!
-!-----------------------------------------------------------------------------
-! Comments:
-!   There are two assumptions concerning the NetCDF files, that are satisfied
-!   with files that are conforming NC convention:
-!     1) The last dimension of the variables used is the time record.
-!     2) Dimensional variables have the same names as corresponding dimensions.
-!-----------------------------------------------------------------------------
-  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
-       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
-       NF90_GET_ATT
-  USE dimphy, ONLY : klon
-  USE phys_state_var_mod, ONLY : pctsrf
-  USE control_mod
-  use pchsp_95_m, only: pchsp_95
-  use pchfe_95_m, only: pchfe_95
-  use arth_m, only: arth
-
-  IMPLICIT NONE
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-#include "indicesol.h"
-#include "iniprint.h"
-!-----------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
-  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
-  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
-  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
-  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
-  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
-  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
-!------------------------------------------------------------------------------
-! Local variables:
-!--- NetCDF
-  INTEGER :: ncid, varid                  ! NetCDF identifiers
-  CHARACTER(LEN=30)               :: dnam       ! dimension name
-!--- dimensions
-  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
-  REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini   ! initial longitudes vector
-  REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini   ! initial latitudes  vector
-  REAL, POINTER,     DIMENSION(:) :: dlon, dlat ! reordered lon/lat  vectors
-!--- fields
-  INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
-  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
-  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
-  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
-  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
-  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
-!--- input files
-  CHARACTER(LEN=20)                 :: cal_in   ! calendar
-  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
-  INTEGER                           :: ndays_in ! number of days
-!--- misc
-  INTEGER :: i, j, k, l                         ! loop counters
-  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
-  CHARACTER(LEN=25)                 :: title    ! for messages
-  LOGICAL                           :: extrp    ! flag for extrapolation
-  LOGICAL                           :: oldice   ! flag for old way ice computation 
-  REAL                              :: chmin, chmax
-  INTEGER ierr
-  integer n_extrap ! number of extrapolated points
-  logical skip
-
-!------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' -------------------------------------
-  NULLIFY(champo)
-
-  SELECT CASE(mode)
-  CASE('RUG'); title='Rugosite'
-  CASE('SIC'); title='Sea-ice'
-  CASE('SST'); title='SST'
-  CASE('ALB'); title='Albedo'
-  END SELECT
-  
-
-  extrp=.FALSE. 
-  oldice=.FALSE.
-  IF ( PRESENT(flag) ) THEN 
-    IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 
-    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 
-  END IF
-
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
-  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
-  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
-  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
-
-!--- Read unit for sea-ice variable only
-  IF (mode=='SIC') THEN
-     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
-     IF(ierr/=NF90_NOERR) THEN
-        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
-        unit_sic='X'
-     ELSE
-        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
-     END IF
-  END IF
-
-!--- Longitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
-
-!--- Latitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
-
-!--- Time (variable is not needed - it is rebuilt - but calendar is)
-  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
-  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
-  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
-  cal_in=' '
-  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
-  IF(ierr/=NF90_NOERR) THEN
-    SELECT CASE(mode)
-      CASE('RUG', 'ALB'); cal_in='360d'
-      CASE('SIC', 'SST'); cal_in='gregorian'
-    END SELECT
-    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
-         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
-  END IF
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
-       cal_in
-
-  
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
-  !--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref, cal_in)
-
-!--- Time vector reconstruction (time vector from file is not trusted)
-!--- If input records are not monthly, time sampling has to be constant !
-  timeyear=mid_months(anneeref, cal_in, lmdep)
-  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
-       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
-
-!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
-  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
-  IF(extrp) ALLOCATE(work(imdep, jmdep))
-
-  IF (prt_level>5) WRITE(lunout, *)
-  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
-  DO l=1, lmdep
-    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
-    CALL ncerr(ierr, fnam)
-    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
-         champ, ibar)
-
-    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
-         work)
-
-    IF(ibar .AND. .NOT.oldice) THEN
-      IF(l==1 .AND. prt_level>5) THEN
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
-        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
-      END IF
-      IF(mode=='RUG') champ=LOG(champ)
-      CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim),     &
-                         rlatv, champint)
-      IF(mode=='RUG') THEN
-        champint=EXP(champint)
-        WHERE(NINT(mask)/=1) champint=0.001
-      END IF
-    ELSE
-      SELECT CASE(mode)
-        CASE('RUG');       CALL rugosite(imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint, mask)
-        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint)
-        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
-                                    iim, jjp1, rlonv, rlatu, champint)
-      END SELECT
-    END IF
-    champtime(:, :, l)=champint
-  END DO
-  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
-
-  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
-  IF(extrp) DEALLOCATE(work)
-
-!--- TIME INTERPOLATION ------------------------------------------------------
-  IF (prt_level>5) THEN
-     WRITE(lunout, *)
-     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
-     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
-     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
-  END IF
-
-  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
-  skip = .false.
-  n_extrap = 0
-  DO j=1, jjp1
-    DO i=1, iim
-      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
-           vc_beg=0., vc_end=0.)
-      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
-           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
-      if (ierr < 0) stop 1
-      n_extrap = n_extrap + ierr
-    END DO
-  END DO
-  if (n_extrap /= 0) then
-     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
-  end if
-  champan(iip1, :, :)=champan(1, :, :)
-  DEALLOCATE(yder, champtime, timeyear)
-
-!--- Checking the result
-  DO j=1, jjp1
-    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
-  END DO
-
-!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
-  IF(mode=='SST') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
-    WHERE(champan<271.38) champan=271.38
-  END IF
-
-!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
-  IF(mode=='SIC') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-
-    IF (unit_sic=='1') THEN
-       ! Nothing to be done for sea-ice field is already in fraction of 1
-       ! This is the case for sea-ice in file cpl_atm_sic.nc
-       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
-    ELSE
-       ! Convert sea ice from percentage to fraction of 1
-       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 
-       champan(:, :, :)=champan(:, :, :)/100.
-    END IF
-
-    champan(iip1, :, :)=champan(1, :, :)
-    WHERE(champan>1.0) champan=1.0
-    WHERE(champan<0.0) champan=0.0
- END IF
-
-!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
-  ALLOCATE(champo(klon, ndays))
-  DO k=1, ndays
-    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
-  END DO
-  DEALLOCATE(champan)
-
-END SUBROUTINE get_2Dfield
-!
-!-------------------------------------------------------------------------------
-
-
-
-!-------------------------------------------------------------------------------
-!
-FUNCTION year_len(y,cal_in)
-!
-!-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER                       :: year_len
-  INTEGER,           INTENT(IN) :: y
-  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
-!-------------------------------------------------------------------------------
-!--- Getting the input calendar to reset at the end of the function
-  CALL ioget_calendar(cal_out)
-
-!--- Unlocking calendar and setting it to wanted one
-  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
-
-!--- Getting the number of days in this year
-  year_len=ioget_year_len(y)
-
-!--- Back to original calendar
-  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-END FUNCTION year_len
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-FUNCTION mid_months(y,cal_in,nm)
-!
-!-------------------------------------------------------------------------------
-  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,                INTENT(IN) :: y               ! year
-  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
-  INTEGER,                INTENT(IN) :: nm              ! months/year number
-  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=99)                  :: mess            ! error message
-  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
-  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
-  INTEGER                            :: m               ! months counter
-  INTEGER                            :: nd              ! number of days
-!-------------------------------------------------------------------------------
-  nd=year_len(y,cal_in)
-
-  IF(nm==12) THEN
-
-  !--- Getting the input calendar to reset at the end of the function
-    CALL ioget_calendar(cal_out)
-
-  !--- Unlocking calendar and setting it to wanted one
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
-
-  !--- Getting the length of each month
-    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
-
-  !--- Back to original calendar
-    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
-
-  ELSE IF(MODULO(nd,nm)/=0) THEN
-    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
-      nm,' months/year. Months number should divide days number.'
-    CALL abort_gcm('mid_months',TRIM(mess),1)
-
-  ELSE
-    mnth=(/(m,m=1,nm,nd/nm)/)
-  END IF
-
-!--- Mid-months times
-  mid_months(1)=0.5*REAL(mnth(1))
-  DO k=2,nm
-    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
-  END DO
-
-END FUNCTION mid_months
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE ncerr(ncres,fnam)
-!
-!-------------------------------------------------------------------------------
-! Purpose: NetCDF errors handling.
-!-------------------------------------------------------------------------------
-  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,          INTENT(IN) :: ncres
-  CHARACTER(LEN=*), INTENT(IN) :: fnam
-!-------------------------------------------------------------------------------
-#include "iniprint.h"
-  IF(ncres/=NF90_NOERR) THEN
-    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
-    CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
-  END IF
-
-END SUBROUTINE ncerr
-!
-!-------------------------------------------------------------------------------
-
-#endif
-! of #ifdef CPP_EARTH
-
-END SUBROUTINE limit_netcdf
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 979)
@@ -0,0 +1,77 @@
+! 
+! $Id: mod_const_mpi.F90 1700 2012-12-20 14:43:19Z lguez $
+!
+MODULE mod_const_mpi
+
+  INTEGER,SAVE :: COMM_LMDZ
+  INTEGER,SAVE :: MPI_REAL_LMDZ
+ 
+
+CONTAINS 
+
+  SUBROUTINE Init_const_mpi
+#ifdef CPP_IOIPSL
+    USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+    USE ioipsl_getincom
+#endif
+
+    IMPLICIT NONE
+#ifdef CPP_MPI
+    INCLUDE 'mpif.h'
+#endif
+    INTEGER             :: ierr
+    INTEGER             :: comp_id
+    INTEGER             :: thread_required
+    INTEGER             :: thread_provided
+    CHARACTER(len = 6)  :: type_ocean
+
+!$OMP MASTER
+    type_ocean = 'force '
+    CALL getin('type_ocean', type_ocean)
+!$OMP END MASTER
+!$OMP BARRIER
+
+    IF (type_ocean=='couple') THEN
+#ifdef CPP_COUPLE
+!$OMP MASTER
+       CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
+       CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
+!$OMP END MASTER
+#endif
+#ifdef CPP_MPI
+      MPI_REAL_LMDZ=MPI_REAL8
+#endif
+    ELSE
+      CALL init_mpi
+    ENDIF
+
+  END SUBROUTINE Init_const_mpi
+  
+  SUBROUTINE Init_mpi
+  IMPLICIT NONE
+#ifdef CPP_MPI
+     INCLUDE 'mpif.h'
+#endif
+    INTEGER             :: ierr
+    INTEGER             :: thread_required
+    INTEGER             :: thread_provided
+
+#ifdef CPP_MPI
+!$OMP MASTER
+      thread_required=MPI_THREAD_SERIALIZED
+
+      CALL MPI_INIT_THREAD(thread_required,thread_provided,ierr)
+      IF (thread_provided < thread_required) THEN
+        PRINT *,'Warning : The multithreaded level of MPI librairy do not provide the requiered level',  &
+                ' in mod_const_mpi::Init_const_mpi'
+      ENDIF
+      COMM_LMDZ=MPI_COMM_WORLD
+      MPI_REAL_LMDZ=MPI_REAL8
+!$OMP END MASTER
+#endif
+
+   END SUBROUTINE Init_mpi
+    
+END MODULE mod_const_mpi
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_para.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_para.F90	(revision 978)
+++ 	(revision )
@@ -1,77 +1,0 @@
-! 
-! $Id: mod_const_para.F90 1279 2009-12-10 09:02:56Z fairhead $
-!
-MODULE mod_const_mpi
-
-  INTEGER,SAVE :: COMM_LMDZ
-  INTEGER,SAVE :: MPI_REAL_LMDZ
- 
-
-CONTAINS 
-
-  SUBROUTINE Init_const_mpi
-#ifdef CPP_IOIPSL
-    USE IOIPSL
-#else
-! if not using IOIPSL, we still need to use (a local version of) getin
-    USE ioipsl_getincom
-#endif
-
-    IMPLICIT NONE
-#ifdef CPP_MPI
-    INCLUDE 'mpif.h'
-#endif
-    INTEGER             :: ierr
-    INTEGER             :: comp_id
-    INTEGER             :: thread_required
-    INTEGER             :: thread_provided
-    CHARACTER(len = 6)  :: type_ocean
-
-!$OMP MASTER
-    type_ocean = 'force '
-    CALL getin('type_ocean', type_ocean)
-!$OMP END MASTER
-!$OMP BARRIER
-
-    IF (type_ocean=='couple') THEN
-#ifdef CPP_COUPLE
-!$OMP MASTER
-       CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
-       CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
-!$OMP END MASTER
-#endif
-#ifdef CPP_MPI
-      MPI_REAL_LMDZ=MPI_REAL8
-#endif
-    ELSE
-      CALL init_mpi
-    ENDIF
-
-  END SUBROUTINE Init_const_mpi
-  
-  SUBROUTINE Init_mpi
-  IMPLICIT NONE
-#ifdef CPP_MPI
-     INCLUDE 'mpif.h'
-#endif
-    INTEGER             :: ierr
-    INTEGER             :: thread_required
-    INTEGER             :: thread_provided
-
-#ifdef CPP_MPI
-!$OMP MASTER
-      thread_required=MPI_THREAD_SERIALIZED
-
-      CALL MPI_INIT_THREAD(thread_required,thread_provided,ierr)
-      IF (thread_provided < thread_required) THEN
-        PRINT *,'Warning : The multithreaded level of MPI librairy do not provide the requiered level',  &
-                ' in mod_const_mpi::Init_const_mpi'
-      ENDIF
-      COMM_LMDZ=MPI_COMM_WORLD
-      MPI_REAL_LMDZ=MPI_REAL8
-!$OMP END MASTER
-#endif
-
-   END SUBROUTINE Init_mpi
-    
-END MODULE mod_const_mpi
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/pres2lev.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/pres2lev.F90	(revision 978)
+++ 	(revision )
@@ -1,74 +1,0 @@
-! $Id: pres2lev.F 1179 2009-06-11 14:18:47Z jghattas $
-!
-!******************************************************
-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
-  
-  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
-  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
-
-  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
-
-  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    
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/pres2lev_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/pres2lev_mod.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/pres2lev_mod.F90	(revision 979)
@@ -0,0 +1,80 @@
+! $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
+  
+  REAL, INTENT(IN) :: po(ni*nj,lmo) ! niveau de pression ancienne grille
+  REAL, INTENT(IN) :: pn(ni*nj,lmn) ! niveau de pression nouvelle grille
+
+  INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontal
+
+  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/filtrez/filtreg.F
===================================================================
--- trunk/LMDZ.COMMON/libf/filtrez/filtreg.F	(revision 978)
+++ trunk/LMDZ.COMMON/libf/filtrez/filtreg.F	(revision 979)
@@ -185,5 +185,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matrinvn(1,1,j),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
@@ -199,5 +199,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matriceun(1,1,j),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
@@ -213,5 +213,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matricevn(1,1,j),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
@@ -231,5 +231,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matrinvs(1,1,j-jfiltsu+1),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
@@ -247,5 +247,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matriceus(1,1,j-jfiltsu+1),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
@@ -262,5 +262,5 @@
                DO j = jdfil,jffil
 #ifdef BLAS
-                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+                  CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, 
      &                 matricevs(1,1,j-jfiltsv+1),
      &                 iim, champ(1,j,1), iip1*nlat, 0.0,
Index: trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/filtrez/filtreg_mod.F90	(revision 979)
@@ -10,6 +10,6 @@
 
   SUBROUTINE inifilr
-  USE mod_filtre_fft
-    !
+  USE mod_filtre_fft, ONLY : use_filtre_fft,Init_filtre_fft
+  USE mod_filtre_fft_loc, ONLY : Init_filtre_fft_loc=>Init_filtre_fft    !
     !    ... H. Upadhyaya, O.Sharma   ...
     !
@@ -541,4 +541,6 @@
        CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
                            coefilv,modfrstv,jfiltnv,jfiltsv)
+       CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
+                           coefilv,modfrstv,jfiltnv,jfiltsv)
     ENDIF
 
Index: trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90	(revision 978)
+++ trunk/LMDZ.COMMON/libf/filtrez/mod_fft.F90	(revision 979)
@@ -3,11 +3,15 @@
 #ifdef FFT_MATHKEISAN
   USE mod_fft_mathkeisan
-#elif FFT_FFTW
+#else
+#ifdef FFT_FFTW
   USE mod_fft_fftw
-#elif FFT_MKL
+#else
+#ifdef FFT_MKL
   USE mod_fft_mkl
 #else
   USE mod_fft_wrapper
 #endif
+#endif
+#endif
 
 END MODULE mod_fft
Index: trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90	(revision 979)
+++ trunk/LMDZ.COMMON/libf/filtrez/mod_filtre_fft_loc.F90	(revision 979)
@@ -0,0 +1,343 @@
+MODULE mod_filtre_fft_loc
+
+  LOGICAL,SAVE :: use_filtre_fft
+  REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
+  REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
+  REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
+
+CONTAINS
+  
+  SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
+    USE mod_fft
+    IMPLICIT NONE
+    include 'dimensions.h'
+    REAL,   INTENT(IN) :: coeffu(iim,jjm)
+    INTEGER,INTENT(IN) :: modfrstu(jjm)
+    INTEGER,INTENT(IN) :: jfiltnu
+    INTEGER,INTENT(IN) :: jfiltsu
+    REAL,   INTENT(IN) :: coeffv(iim,jjm)
+    INTEGER,INTENT(IN) :: modfrstv(jjm)
+    INTEGER,INTENT(IN) :: jfiltnv
+    INTEGER,INTENT(IN) :: jfiltsv
+    
+    INTEGER            :: index_vp(iim)
+    INTEGER            :: i,j
+    
+    index_vp(1)=1
+    DO i=1,iim/2
+      index_vp(i+1)=i*2
+    ENDDO
+    
+    DO i=1,iim/2-1
+      index_vp(iim/2+i+1)=iim-2*i+1
+    ENDDO
+    
+    ALLOCATE(Filtre_u(iim,jjm))
+    ALLOCATE(Filtre_v(iim,jjm))
+    ALLOCATE(Filtre_inv(iim,jjm))
+  
+    
+    DO j=2,jfiltnu
+      DO i=1,iim
+        IF (index_vp(i) < modfrstu(j)) THEN
+          Filtre_u(i,j)=0
+        ELSE
+          Filtre_u(i,j)=coeffu(index_vp(i),j)
+        ENDIF
+      ENDDO
+    ENDDO
+    
+    DO j=jfiltsu,jjm
+      DO i=1,iim
+        IF (index_vp(i) < modfrstu(j)) THEN
+          Filtre_u(i,j)=0
+        ELSE
+          Filtre_u(i,j)=coeffu(index_vp(i),j)
+        ENDIF
+      ENDDO
+    ENDDO
+ 
+    DO j=1,jfiltnv
+      DO i=1,iim
+        IF (index_vp(i) < modfrstv(j)) THEN
+          Filtre_v(i,j)=0
+        ELSE
+          Filtre_v(i,j)=coeffv(index_vp(i),j)
+        ENDIF
+      ENDDO
+    ENDDO
+   
+    DO j=jfiltsv,jjm
+      DO i=1,iim
+        IF (index_vp(i) < modfrstv(j)) THEN
+          Filtre_v(i,j)=0
+        ELSE
+          Filtre_v(i,j)=coeffv(index_vp(i),j)
+        ENDIF
+      ENDDO
+    ENDDO
+         
+    DO j=2,jfiltnu
+      DO i=1,iim
+        IF (index_vp(i) < modfrstu(j)) THEN
+          Filtre_inv(i,j)=0
+        ELSE
+          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
+        ENDIF
+      ENDDO
+    ENDDO
+
+    DO j=jfiltsu,jjm
+      DO i=1,iim
+        IF (index_vp(i) < modfrstu(j)) THEN
+          Filtre_inv(i,j)=0
+        ELSE
+          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
+        ENDIF
+      ENDDO
+    ENDDO
+    
+    
+!    CALL Init_FFT(iim,(jjm+1)*(llm+1))
+        
+    
+  END SUBROUTINE Init_filtre_fft
+  
+  SUBROUTINE Filtre_u_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
+    USE mod_fft
+#ifdef CPP_PARA
+    USE parallel,ONLY : OMP_CHUNK
+#endif
+    IMPLICIT NONE
+    include 'dimensions.h'
+    INTEGER,INTENT(IN) :: jjb
+    INTEGER,INTENT(IN) :: jje
+    INTEGER,INTENT(IN) :: jj_begin
+    INTEGER,INTENT(IN) :: jj_end
+    INTEGER,INTENT(IN) :: nbniv
+    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
+
+    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
+    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
+    INTEGER            :: nb_vect
+    INTEGER :: i,j,l
+    INTEGER :: ll_nb
+!    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
+    
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+
+    nb_vect=(jj_end-jj_begin+1)*ll_nb
+
+!    vect_tmp=vect
+
+    CALL FFT_forward(vect,TF_vect,nb_vect)
+
+!    CALL FFT_forward(vect,TF_vect_test,nb_vect)
+!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
+!      DO j=1,jj_end-jj_begin+1
+!      DO i=1,iim/2+1
+!         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
+!       ENDDO
+!      ENDDO
+
+    DO l=1,ll_nb
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim/2+1
+          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
+        ENDDO
+      ENDDO
+    ENDDO
+       
+    CALL FFT_backward(TF_vect,vect,nb_vect)
+!    CALL FFT_backward(TF_vect_test,vect_test,nb_vect)
+          
+!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
+!      DO j=1,jj_end-jj_begin+1
+!         DO i=1,iim
+!           PRINT *,"====",i,j,"----->",vect_test(i,j,1)
+!         ENDDO
+!      ENDDO
+
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+
+  END SUBROUTINE Filtre_u_fft
+  
+
+  SUBROUTINE Filtre_v_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
+    USE mod_fft
+#ifdef CPP_PARA
+    USE parallel,ONLY : OMP_CHUNK
+#endif
+    IMPLICIT NONE
+    INCLUDE 'dimensions.h'
+    INTEGER,INTENT(IN) :: jjb
+    INTEGER,INTENT(IN) :: jje
+    INTEGER,INTENT(IN) :: jj_begin
+    INTEGER,INTENT(IN) :: jj_end
+    INTEGER,INTENT(IN) :: nbniv
+    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
+
+    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
+    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
+    INTEGER            :: nb_vect
+    INTEGER :: i,j,l
+    INTEGER :: ll_nb
+    
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+
+    
+    nb_vect=(jj_end-jj_begin+1)*ll_nb
+
+    CALL FFT_forward(vect,TF_vect,nb_vect)
+  
+    DO l=1,ll_nb
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim/2+1
+          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
+        ENDDO
+      ENDDO
+    ENDDO
+  
+    CALL FFT_backward(TF_vect,vect,nb_vect)
+    
+    
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+  
+  END SUBROUTINE Filtre_v_fft
+
+
+  SUBROUTINE Filtre_inv_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
+    USE mod_fft
+#ifdef CPP_PARA
+    USE parallel,ONLY : OMP_CHUNK
+#endif
+    IMPLICIT NONE
+    INCLUDE 'dimensions.h'
+    INTEGER,INTENT(IN) :: jjb
+    INTEGER,INTENT(IN) :: jje
+    INTEGER,INTENT(IN) :: jj_begin
+    INTEGER,INTENT(IN) :: jj_end
+    INTEGER,INTENT(IN) :: nbniv
+    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
+
+    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
+    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
+    INTEGER            :: nb_vect
+    INTEGER :: i,j,l
+    INTEGER :: ll_nb
+    
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+
+    nb_vect=(jj_end-jj_begin+1)*ll_nb
+
+    CALL FFT_forward(vect,TF_vect,nb_vect)
+  
+    DO l=1,ll_nb
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim/2+1
+          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
+        ENDDO
+      ENDDO
+    ENDDO
+  
+    CALL FFT_backward(TF_vect,vect,nb_vect)
+
+    ll_nb=0
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+    DO l=1,nbniv
+      ll_nb=ll_nb+1
+      DO j=1,jj_end-jj_begin+1
+        DO i=1,iim+1
+          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP END DO NOWAIT
+
+  END SUBROUTINE Filtre_inv_fft  
+  
+  
+!  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
+!  IMPLICIT NONE
+!    INTEGER,INTENT(IN)  :: nbniv
+!    INTEGER,INTENT(OUT) :: ll_index(nbniv)
+!    INTEGER,INTENT(OUT) :: ll_nb
+!
+!    INTEGER :: l,ll_begin, ll_end
+!   INTEGER :: omp_rank,omp_size
+!   INTEGER :: OMP_GET_NUM_THREADS
+!   INTEGER :: omp_chunk
+!   EXTERNAL OMP_GET_NUM_THREADS
+!   INTEGER :: OMP_GET_THREAD_NUM
+!   EXTERNAL OMP_GET_THREAD_NUM
+!
+!   
+!   omp_size=OMP_GET_NUM_THREADS()
+!   omp_rank=OMP_GET_THREAD_NUM()    
+!   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
+!   
+!   ll_begin=omp_rank*OMP_CHUNK+1
+!   ll_nb=0
+!   DO WHILE (ll_begin<=nbniv)
+!     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
+!     DO l=ll_begin,ll_end
+!       ll_nb=ll_nb+1
+!       ll_index(ll_nb)=l
+!     ENDDO
+!     ll_begin=ll_begin+omp_size*OMP_CHUNK
+!   ENDDO
+!  
+!  END SUBROUTINE get_ll_index
+   
+END MODULE mod_filtre_fft_loc
+ 
Index: trunk/LMDZ.COMMON/makelmdz
===================================================================
--- trunk/LMDZ.COMMON/makelmdz	(revision 978)
+++ trunk/LMDZ.COMMON/makelmdz	(revision 979)
@@ -19,4 +19,5 @@
 chimie=false
 parallel=none
+paramem="par"
 compil_mod=prod
 io=ioipsl
@@ -29,5 +30,5 @@
 
 # guess a default 'arch'
-arch="g95" # start with assuming we're on a Linux/Unix machine with g95
+arch="local" # start with assuming we're on a local machine with local arch file
 ## try to recognise machine and infer arch from it
 machine=`hostname`
@@ -126,4 +127,5 @@
 [-cpp CPP_KEY]             : additional preprocessing definitions
 [-adjnt]                   : adjoint model, not operational ...
+[-mem]                     : reduced memory dynamics (if in parallel mode)
 [-filtre NOMFILTRE]        : use filtre from libf/NOMFILTRE (default: filtrez)
 [-link LINKS]              : additional links with other libraries
@@ -138,5 +140,5 @@
       
       "-p")
-	  physique="$2" ;  shift ; shift ;;
+	  physique="$2" ; shift ; shift ;;
 
       "-s")
@@ -186,4 +188,7 @@
       "-cosp")
           cosp="$2" ; shift ; shift ;;
+      
+      "-mem")
+          paramem="mem" ; shift ;;
 
       "-filtre")
@@ -424,4 +429,28 @@
 ########################################################################
 
+cd $LIBFGCM/grid
+if [[ -f dimensions.h ]]
+then
+  echo 'ATTENTION: vous etes sans doute en train de compiler le modele par ailleurs'
+  echo "Attendez que la premiere compilation soit terminee pour relancer la suivante."
+  echo "Si vous etes sur que vous ne compilez pas le modele par ailleurs,"
+  echo  vous pouvez continuer en repondant oui.
+  echo "Voulez-vous vraiment continuer?"
+  echo ""
+  echo "WARNING: you are probably already compiling the model somewhere else."
+  echo "Wait until the first compilation is finished before launching this one."
+  echo "If you are sure that you are not compiling elsewhere, just answer "
+  echo "yes (or 'oui') to the question below to proceed."
+  echo "Do you wish to continue?"
+  read reponse
+  if [[ $reponse == "oui" || $reponse == "yes" ]]
+  then
+    \rm -f $LIBFGCM/grid/dimensions.h
+  else
+    exit
+  fi
+fi
+
+
 cd $LIBFGCM/grid/dimension
 ./makdim $dim
@@ -496,6 +525,12 @@
 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
 
@@ -634,4 +669,8 @@
 then
   INCLUDE="$INCLUDE"' -I$(LIBF)/phy'"$physique"
+  phys="PHYS=$physique"
+else
+  # trick to avoid having a physics-related library in the makefile
+  phys="L_PHY= LIBPHY="
 fi
 if (( $dimension == 1 ))
@@ -674,5 +713,5 @@
 $dyn \
 LIBO=${LIBOGCM}/${nomlib} \
-PHYS=$physique \
+$phys \
 DIM=$dimc \
 FLAG_PARA=$FLAG_PARA \
@@ -705,5 +744,5 @@
 $dyn \
 LIBO=${LIBOGCM}/${nomlib} \
-PHYS=$physique \
+$phys \
 DIM=$dimc \
 FLAG_PARA=$FLAG_PARA \
@@ -725,7 +764,7 @@
 PROG=$code
 
-if [[ -r $libf/grid/dimensions.h ]]
+if [[ -r $LIBFGCM/grid/dimensions.h ]]
 then
   # Cleanup: remove dimension.h file
-  \rm -f $libf/grid/dimensions.h
-fi
+  \rm -f $LIBFGCM/grid/dimensions.h
+fi
Index: trunk/LMDZ.COMMON/makelmdz_fcm
===================================================================
--- trunk/LMDZ.COMMON/makelmdz_fcm	(revision 978)
+++ trunk/LMDZ.COMMON/makelmdz_fcm	(revision 979)
@@ -25,4 +25,5 @@
 chemistry=false
 parallel=none
+paramem="par"
 compil_mod=prod
 io=ioipsl
@@ -80,4 +81,5 @@
 [-cpp CPP_KEY]             : additional preprocessing definitions
 [-adjnt]                   : adjoint model, not operational ...
+[-mem]                     : reduced memory dynamics (if in parallel mode)
 [-filtre NOMFILTRE]        : use filtre from libf/NOMFILTRE (default: filtrez)
 [-link LINKS]              : additional links with other libraries
@@ -142,4 +144,6 @@
           cosp="$2" ; shift ; shift ;;
 
+      "-mem")
+          paramem="mem" ; shift ;;
 
       "-filtre")
@@ -164,10 +168,18 @@
 
 ###############################################################
-# mettre le chemin du fcm dans le path
+# path to fcm 
 ###############################################################
+# handle case when provided path to fcm was given as a relative
+# path (from makelmdz_fcm script directory) and not an absolute path
+if [[ ${fcm_path:0:1} != "/" ]] ; then
+  # prepend with makelmdz_fcm location
+  fcm_path=$(cd $(dirname $0) ; pwd)"/"${fcm_path}
+fi
+
+# add fcm_path to PATH
 export PATH=${fcm_path}:${PATH}
 
-echo "Chemin du fcm utlise :" 
-which fcm
+echo "Path to fcm:"
+echo ${fcm_path}
 
 ###############################################################
@@ -295,4 +307,28 @@
 # on cree ou remplace le fichier des dimensions
 ########################################################################
+
+cd $LIBFGCM/grid
+if [[ -f dimensions.h ]]
+then
+  echo 'ATTENTION: vous etes sans doute en train de compiler le modele par ailleurs'
+  echo "Attendez que la premiere compilation soit terminee pour relancer la suivante."
+  echo "Si vous etes sur que vous ne compilez pas le modele par ailleurs,"
+  echo  vous pouvez continuer en repondant oui.
+  echo "Voulez-vous vraiment continuer?"
+  echo ""
+  echo "WARNING: you are probably already compiling the model somewhere else."
+  echo "Wait until the first compilation is finished before launching this one."
+  echo "If you are sure that you are not compiling elsewhere, just answer "
+  echo "yes (or 'oui') to the question below to proceed."
+  echo "Do you wish to continue?"
+  read reponse
+  if [[ $reponse == "oui" || $reponse == "yes" ]]
+  then
+    \rm -f $LIBFGCM/grid/dimensions.h
+  else
+    exit
+  fi
+fi
+
 
 cd $LIBFGCM/grid/dimension
@@ -393,5 +429,9 @@
 then
   SUFF_NAME=${SUFF_NAME}_para
-  DYN=dyn${dimc}dpar
+  DYN=dyn${dimc}d${paramem}
+  if [[ "$paramem" == "mem" ]]
+  then
+   SUFF_NAME=${SUFF_NAME}_${paramem}
+  fi
 else
   SUFF_NAME=${SUFF_NAME}_seq
@@ -463,2 +503,8 @@
 ln -s $LIBOGCM/${arch}${SUFF_NAME}/.config config
 ln -s $LIBOGCM/${arch}${SUFF_NAME}/.config/tmp tmp_src
+
+if [[ -r $LIBFGCM/grid/dimensions.h ]]
+then
+  # Cleanup: remove dimension.h file
+  \rm -f $LIBFGCM/grid/dimensions.h
+fi
