Index: trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/NS_tridag.F90	(revision 3498)
@@ -27,4 +27,4 @@
       enddo
 !      return
-      END
+      END SUBROUTINE tridag
 !  (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31&#5(+.
Index: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3498)
@@ -40,5 +40,5 @@
 
 !=======================================================================
-SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
+SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
                                m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
 
@@ -47,5 +47,5 @@
 ! Inputs
 integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
-real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! tendancies on the glaciers [1]
+real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! water ice at the surface [kg/m^2]
 real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! co2 ice at the surface [kg/m^2]
@@ -66,7 +66,7 @@
 ! -------------
 ! Compute H2O adsorption, then CO2 adsorption
-call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
                             theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg)
-call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
+call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
                             tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
 
@@ -74,5 +74,5 @@
 
 !=======================================================================
-SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
                                   theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg)
 
@@ -100,5 +100,5 @@
 ! Inputs
 integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
-real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
+real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
 real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
@@ -151,6 +151,6 @@
     do ig = 1,ngrid
         do islope = 1,nslope
-            if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
-            if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
+            if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
+            if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
         enddo
     enddo
@@ -224,5 +224,5 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
+SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
 
 use comsoil_h_PEM,         only: layer_PEM, index_breccia, index_breccia
@@ -244,5 +244,5 @@
 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Mean Soil Temperature [K]
 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Mean Thermal Inertia [USI]
-real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
+real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
 real,    dimension(ngrid,timelen),          intent(in) :: q_co2, q_h2o                       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
@@ -292,6 +292,6 @@
     do ig = 1,ngrid
         do islope = 1,nslope
-            if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
-            if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
+            if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
+            if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
         enddo
     enddo
@@ -317,5 +317,5 @@
 
     ! 1. Compute the fraction of the pores occupied by H2O
-    call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+    call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
                                 theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
 
Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3498)
@@ -472,2 +472,7 @@
     - The execution command line in the job scripts that should be modified by the user according to the set-up is now given as an argument at the beginning to be more identifiable and adaptable;
     - Making the job scripts more robust to detect a successful end.
+
+== 07/11/2024 == JBC
+- Correction of the variable name for the ice table depth in "pemetat0.F90". So it is now got as intended from the "startpem.nc" file;
+- Renaming of the tendencies in the PEM with the prefix 'd_' instead of 'tend_';
+- Modification of the PEM time step type from integer to real. As a consequence, all time variables are now of real type. This change adds the possibility to consider fractions of year as time step.
Index: trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90	(revision 3498)
@@ -7,5 +7,5 @@
 !=======================================================================
 
-SUBROUTINE compute_tend(ngrid,nslope,min_ice,tendencies_ice)
+SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice)
 
 implicit none
@@ -25,14 +25,14 @@
 
 !   OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! Difference between the minima = evolution of perennial ice
+real, dimension(ngrid,nslope), intent(out) :: d_ice ! Difference between the minima = evolution of perennial ice
 !=======================================================================
 ! We compute the difference
-tendencies_ice = min_ice(:,:,2) - min_ice(:,:,1)
+d_ice = min_ice(:,:,2) - min_ice(:,:,1)
 
 ! If the difference is too small, then there is no evolution
-where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
+where (abs(d_ice) < 1.e-10) d_ice = 0.
 
 ! If the minimum over the last year is 0, then we have no perennial ice
-where (abs(min_ice(:,:,2)) < 1.e-10) tendencies_ice = 0.
+where (abs(min_ice(:,:,2)) < 1.e-10) d_ice = 0.
 
 END SUBROUTINE compute_tend
Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3498)
@@ -23,6 +23,6 @@
 #endif
 
-use time_evol_mod,         only: year_bp_ini, dt_pem, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
-                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years, ecritpem
+use time_evol_mod,         only: year_bp_ini, dt, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
+                                 evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years, ecritpem
 use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
 use adsorption_mod,        only: adsorption_pem
@@ -34,5 +34,5 @@
 implicit none
 
-integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
+real, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
 
 character(20), parameter :: modname ='conf_pem'
@@ -64,5 +64,5 @@
 year_bp_ini = 0.
 call getin('year_earth_bp_ini',year_earth_bp_ini)
-year_bp_ini = int(year_earth_bp_ini/convert_years)
+year_bp_ini = year_earth_bp_ini/convert_years
 
 var_obl = .true.
@@ -91,6 +91,6 @@
 call getin('ps_criterion',ps_criterion)
 
-dt_pem = 1
-call getin('dt_pem', dt_pem)
+dt = 1.
+call getin('dt',dt)
 
 !#---------- Subsurface parameters ----------#
@@ -128,4 +128,8 @@
     call abort_physic(modname,"Ice table must be used when soil_pem = T",1)
 endif
+if (icetable_equilibrium .or. icetable_dynamic) then
+    write(*,*) 'Ice table is asked to be computed both by the equilibrium and dynamic method.'
+    write(*,*) 'The dynamic method is then chosen.'
+endif
 
 if ((.not. soil_pem) .and. adsorption_pem) then
@@ -144,8 +148,8 @@
 endif
 
-if (evol_orbit_pem .and. year_bp_ini == 0.) then
-    write(*,*) 'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
+if (evol_orbit_pem .and. abs(year_bp_ini) < 1.e-10) then
+    write(*,*) 'You want to follow the file "obl_ecc_lsp.asc" for changing orbital parameters,'
     write(*,*) 'but you did not specify from which year to start.'
-    call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
+    call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini = 0",1)
 endif
 
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/README
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/README	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/README	(revision 3498)
@@ -7,5 +7,5 @@
       (ii)  nPCM_ini -> the number of initial PCM runs (at least 2);
       (iii) nPCM -> the number of PCM runs between each PEM run (usually 2);
-      (iv)  mode -> the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on supercomputer.
+      (iv)  mode -> the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on a supercomputer.
   The script can take an argument:
       - If there is no argument, then the script initiates a PEM simulation from scratch.
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh	(revision 3498)
@@ -14,6 +14,6 @@
 ###############################################
 # Set the number of years to be simulated, either Martian or Earth years:
-n_mars_years=100
-#n_earth_years=300
+n_mars_years=100.
+#n_earth_years=300.
 
 # Set the number of initial PCM runs (>= 2):
@@ -23,5 +23,5 @@
 nPCM=2
 
-# Set the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on supercomputer:
+# Set the launching mode (0 = "processing scripts"; any other values = "submitting jobs"). The former option is usually used to process the script on a local machine while the latter is used to submit jobs on a supercomputer:
 mode=1
 ########################################################################
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/lib_launchPEM.sh	(revision 3498)
@@ -52,10 +52,10 @@
 
     if [ -v n_mars_years ] && [ ! -z "$n_mars_years" ]; then
-        if [ $n_mars_years -lt 1 ]; then
+        if [ $n_mars_years -lt 0. ]; then
             echo "Error: the value of 'n_mars_years' must be >0!"
             errlaunch
         fi
     elif [ -v n_earth_years ] && [ ! -z "$n_earth_years" ]; then
-        if [ $n_earth_years -lt 1 ]; then
+        if [ $n_earth_years -lt 0. ]; then
             echo "Error: the value of 'n_earth_years' must be >0!"
             errlaunch
@@ -132,5 +132,5 @@
         echo "Number of years to be simulated: $n_myear Martian years."
     elif [ -v n_earth_years ]; then
-        n_myear=$(echo "($n_earth_years/$convert_years + 0.999999)/1" | bc) # Ceiling of n_earth_years/convert_years
+        n_myear=$(echo "$n_earth_years/$convert_years" | bc -l)
         echo "Number of years to be simulated: $n_earth_years Earth years = $n_myear Martian years."
     fi
@@ -141,5 +141,5 @@
     echo "This is a chained simulation for PEM and PCM runs in $dir on $machine by $user."
     convertyears
-    i_myear=0
+    i_myear=0.
     iPEM=1
     iPCM=1
@@ -176,5 +176,5 @@
                 errlaunch
             fi
-        else # Mode: launching jobs
+        else # Mode: submitting jobs
             cp PCMrun.job PCMrun${iPCM}.job
             sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPCM}/" PCMrun${iPCM}.job
@@ -187,5 +187,5 @@
         fi
         ((iPCM++))
-        ((i_myear++))
+        i_myear = $(echo "$i_myear + 1." | bc - l)
         ((ii++))
     else
@@ -201,5 +201,5 @@
                     errlaunch
                 fi
-            else # Mode: launching jobs
+            else # Mode: submitting jobs
                 cp PCMrun.job PCMrun${iPCM}.job
                 sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPCM}/" PCMrun${iPCM}.job
@@ -209,5 +209,5 @@
             fi
             ((iPCM++))
-            ((i_myear++))
+            i_myear = $(echo "$i_myear + 1." | bc - l)
         else
             endlaunch
@@ -226,5 +226,5 @@
                 errlaunch
             fi
-        else # Mode: launching jobs
+        else # Mode: submitting jobs
             sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
             jobID=$(eval "$submit_job PEMrun.job")
@@ -255,5 +255,5 @@
                 errlaunch
             fi
-        else # Mode: launching jobs
+        else # Mode: submitting jobs
             sed -i -E "s/($name_job[^0-9]*[0-9]*[^0-9]*)[0-9]+$/\1${iPEM}/" PEMrun.job
             jobID=$(eval "$submit_dependjob=afterok:${jobID} PEMrun.job")
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3498)
@@ -39,6 +39,6 @@
 # ps_criterion = 0.15
 
-# Time step length of the PEM in Martian years? Default = 1
-# dt_pem=1
+# Time step length of the PEM in Martian years? Default = 1.
+# dt=1
 
 #---------- Subsurface parameters ----------#
Index: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3498)
@@ -7,7 +7,7 @@
 !=======================================================================
 
-SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
+SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
 
-use time_evol_mod, only: dt_pem
+use time_evol_mod, only: dt
 
 implicit none
@@ -25,5 +25,5 @@
 !   OUTPUT
 real, dimension(ngrid,nslope), intent(inout) :: co2_ice      ! Previous and actual density of CO2 ice
-real, dimension(ngrid,nslope), intent(inout) :: tend_co2_ice ! Evolution of perennial ice over one year
+real, dimension(ngrid,nslope), intent(inout) :: d_co2ice ! Evolution of perennial ice over one year
 
 !   local:
@@ -35,8 +35,8 @@
 
 co2_ice_old = co2_ice
-co2_ice = co2_ice + tend_co2_ice*dt_pem
+co2_ice = co2_ice + d_co2ice*dt
 where (co2_ice < 0.)
     co2_ice = 0.
-    tend_co2_ice = -co2_ice_old/dt_pem
+    d_co2ice = -co2_ice_old/dt
 end where
 
@@ -45,7 +45,7 @@
 !=======================================================================
 
-SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
+SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
 
-use time_evol_mod,          only: dt_pem
+use time_evol_mod,          only: dt
 use comslope_mod,           only: subslope_dist, def_slope_mean
 
@@ -73,5 +73,5 @@
 !   OUTPUT
 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice      ! Previous and actual density of h2o ice (kg/m^2)
-real, dimension(ngrid,nslope), intent(inout) :: tend_h2o_ice ! Evolution of perennial ice over one year (kg/m^2/year)
+real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year)
 integer,                       intent(inout) :: stopPEM      ! Stopping criterion code
 
@@ -80,5 +80,5 @@
 integer                       :: i, islope                                           ! Loop variables
 real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
-real, dimension(ngrid,nslope) :: new_tendencies                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
+real, dimension(ngrid,nslope) :: new_tend                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
 !=======================================================================
 if (ngrid /= 1) then ! Not in 1D
@@ -99,8 +99,8 @@
         do islope = 1,nslope
             if (h2o_ice(i,islope) > 0.) then
-                if (tend_h2o_ice(i,islope) > 0.) then
-                    pos_tend = pos_tend + tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
+                if (d_h2oice(i,islope) > 0.) then
+                    pos_tend = pos_tend + d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                 else
-                    neg_tend = neg_tend - tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
+                    neg_tend = neg_tend - d_h2oice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                 endif
             endif
@@ -108,5 +108,5 @@
     enddo
 
-    new_tendencies = 0.
+    new_tend = 0.
     if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
         write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
@@ -118,14 +118,14 @@
         ! We adapt the tendencies to conserve h2o and do only exchange between grid points
         if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
-            where (tend_h2o_ice < 0.) ! We lower the sublimating rate by a coefficient
-                new_tendencies = tend_h2o_ice*pos_tend/neg_tend
+            where (d_h2oice < 0.) ! We lower the sublimating rate by a coefficient
+                new_tend = d_h2oice*pos_tend/neg_tend
             elsewhere ! We don't change the condensing rate
-                new_tendencies = tend_h2o_ice
+                new_tend = d_h2oice
             end where
         else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
-            where (tend_h2o_ice < 0.) ! We don't change the sublimating rate
-                new_tendencies = tend_h2o_ice
+            where (d_h2oice < 0.) ! We don't change the sublimating rate
+                new_tend = d_h2oice
             elsewhere ! We lower the condensing rate by a coefficient
-                new_tendencies = tend_h2o_ice*neg_tend/pos_tend
+                new_tend = d_h2oice*neg_tend/pos_tend
             end where
         endif
@@ -133,5 +133,5 @@
 
     ! Evolution of the h2o ice for each physical point
-    h2o_ice = h2o_ice + new_tendencies*dt_pem
+    h2o_ice = h2o_ice + new_tend*dt
 
     ! We compute the amount of h2o that is sublimated in excess
@@ -142,5 +142,5 @@
                 negative_part = negative_part - h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
                 h2o_ice(i,islope) = 0.
-                tend_h2o_ice(i,islope) = 0.
+                d_h2oice(i,islope) = 0.
             endif
         enddo
@@ -157,9 +157,9 @@
     do islope = 1,nslope
         do i = 1,ngrid
-             if (new_tendencies(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
+             if (new_tend(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tend(i,islope)*real_coefficient*dt*cos(def_slope_mean(islope)*pi/180.)
         enddo
      enddo
 else ! ngrid == 1, i.e. in 1D
-    h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
+    h2o_ice = h2o_ice + d_h2oice*dt
 endif
 
Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3498)
@@ -10,9 +10,9 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-logical, save                           :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
-logical, save                           :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
-real, allocatable, dimension(:,:)       :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
-real, allocatable, dimension(:,:)       :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
-real, allocatable, dimension(:,:,:)     :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
+logical                             :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
+logical                             :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
+real, allocatable, dimension(:,:)   :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
+real, allocatable, dimension(:,:)   :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
+real, allocatable, dimension(:,:,:) :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
 
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3498)
@@ -9,5 +9,5 @@
 !=======================================================================
 
-SUBROUTINE info_PEM(year_iter,stopPEM,i_myear,n_myear)
+SUBROUTINE info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
 
 !=======================================================================
@@ -24,6 +24,7 @@
 
 !----- Arguments
-integer, intent(in) :: year_iter, stopPEM ! # of year and reason to stop
-integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
+integer, intent(in) :: stopPEM          ! Reason to stop
+real,    intent(in) :: i_myear_leg      ! # of years
+real,    intent(in) :: i_myear, n_myear ! Current simulated Martian year and maximum number of Martian years to be simulated
 
 !----- Local variables
@@ -35,6 +36,6 @@
 inquire(file = 'info_PEM.txt',exist = ok)
 if (ok) then
-    write(ich1,'(i0)') i_myear
-    write(ich2,'(i0)') n_myear
+    write(ich1,'(f0.4)') i_myear
+    write(ich2,'(f0.4)') n_myear
     write(fch,'(f0.4)') convert_years ! 4 digits to the right of the decimal point to respect the precision of Martian year in "launch_pem.sh"
     write(ich3,'(i0)') iPCM
@@ -51,5 +52,5 @@
     ! Martian date, Number of Martians years done by the PEM run, Number of Martians years done by the chainded simulation, Code of the stopping criterion
     ! The conversion ratio from Planetary years to Earth years is given in the header of the file
-    write(1,*) year_bp_ini + i_myear, year_iter, i_myear, stopPEM
+    write(1,*) year_bp_ini + i_myear, i_myear_leg, i_myear, stopPEM
     close(1)
 else
Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3498)
@@ -18,5 +18,5 @@
 
 ! Physical parameters
-real, parameter :: tend_dust = 5.78e-2          ! Tendency of dust [kg.m-2.y-1]
+real, parameter :: d_dust = 5.78e-2             ! Tendency of dust [kg.m-2.y-1]
 real, parameter :: dry_lag_porosity = 0.5       ! Porosity of dust lag
 real, parameter :: dry_regolith_porosity = 0.45 ! Porosity of regolith
@@ -452,5 +452,5 @@
 !=======================================================================
 ! Layering algorithm
-SUBROUTINE make_layering(this,tend_co2ice,tend_h2oice,tend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
+SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
 
 implicit none
@@ -461,8 +461,8 @@
 logical,                intent(inout) :: new_str, new_lag1, new_lag2
 integer,                intent(inout) :: stopPEM
-real,     intent(in) :: tend_co2ice, tend_h2oice, tend_dust
-
-!---- Local variables
-real                   :: htend_co2ice, htend_h2oice, htend_dust
+real,     intent(in) :: d_co2ice, d_h2oice, d_dust
+
+!---- Local variables
+real                   :: dh_co2ice, dh_h2oice, dh_dust
 real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
 real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
@@ -471,12 +471,12 @@
 
 !---- Code
-htend_co2ice = tend_co2ice/rho_co2ice
-htend_h2oice = tend_h2oice/rho_h2oice
-htend_dust = tend_dust/rho_dust
-
-if (htend_dust < 0.) then ! Dust lifting only
-    if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
+dh_co2ice = d_co2ice/rho_co2ice
+dh_h2oice = d_h2oice/rho_h2oice
+dh_dust = d_dust/rho_dust
+
+if (dh_dust < 0.) then ! Dust lifting only
+    if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
     write(*,'(a)') ' Stratification -> Dust lifting'
-    h2lift_tot = abs(htend_dust)
+    h2lift_tot = abs(dh_dust)
     do while (h2lift_tot > 0. .and. associated(current1))
         ! Is the considered stratum a dust lag?
@@ -508,11 +508,11 @@
 
 !------ Dust sedimentation only
-    if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
+    if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
         write(*,'(a)') ' Stratification -> Dust sedimentation'
         ! New stratum at the layering top by sedimentation of dust
-        thickness = htend_dust/(1. - dry_regolith_porosity)
+        thickness = dh_dust/(1. - dry_regolith_porosity)
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
                  new_str = .false.
              else
@@ -523,11 +523,11 @@
 
 !------ Condensation of CO2 ice + H2O ice
-    else if ((htend_co2ice >= 0. .and. htend_h2oice > 0.) .or. (htend_co2ice > 0. .and. htend_h2oice >= 0.)) then
+    else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then
         write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
         ! New stratum at the layering top by condensation of CO2 and H2O ice
-        thickness = htend_co2ice/(1. - co2ice_porosity) + htend_h2oice/(1. - h2oice_porosity) + htend_dust
+        thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
                  new_str = .false.
              else
@@ -538,8 +538,8 @@
 
 !------ Sublimation of CO2 ice + Condensation of H2O ice
-    else if (htend_co2ice < 0. .and. htend_h2oice > 0.) then
+    else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
         write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
-        h2subl_tot = abs(htend_co2ice)
+        h2subl_tot = abs(dh_co2ice)
         do while (h2subl_tot > 0. .and. associated(current1))
             h_co2ice_old = current1%co2ice_volfrac*current1%thickness
@@ -585,8 +585,8 @@
         endif
         ! New stratum at the layering top by condensation of H2O ice
-        thickness = htend_h2oice/(1. - h2oice_porosity) + htend_dust
+        thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust
         if (thickness > 0.) then ! Only if the stratum is building up
              if (new_str) then
-                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
+                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
                  new_str = .false.
              else
@@ -597,8 +597,8 @@
 
 !------ Sublimation of CO2 ice + H2O ice
-    else if ((htend_co2ice <= 0. .and. htend_h2oice < 0.) .or. (htend_co2ice < 0. .and. htend_h2oice <= 0.)) then
+    else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then
         write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
-        h2subl_tot = abs(htend_co2ice)
+        h2subl_tot = abs(dh_co2ice)
         do while (h2subl_tot > 0. .and. associated(current1))
             h_co2ice_old = current1%co2ice_volfrac*current1%thickness
@@ -644,5 +644,5 @@
         endif
         ! H2O ice sublimation in the considered stratum + New stratum for dust lag
-        h2subl_tot = abs(htend_h2oice)
+        h2subl_tot = abs(dh_h2oice)
         do while (h2subl_tot > 0. .and. associated(current2))
             h_co2ice_old = current2%co2ice_volfrac*current2%thickness
@@ -689,5 +689,5 @@
 
 !------ Condensation of CO2 ice + Sublimation of H2O ice
-    else if (htend_co2ice > 0. .and. htend_h2oice < 0.) then
+    else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then
         error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
     endif
Index: trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3498)
@@ -15,5 +15,5 @@
 !=======================================================================
 
-SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
+SUBROUTINE orbit_param_criterion(i_myear,n_myear_leg)
 
 #ifdef CPP_IOIPSL
@@ -38,20 +38,20 @@
 ! Input Variables
 !--------------------------------------------------------
-integer, intent(in) :: i_myear ! Martian year of the simulation
+real, intent(in) :: i_myear ! Martian year of the simulation
 
 !--------------------------------------------------------
 ! Output Variables
 !--------------------------------------------------------
-integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
+real, intent(out) :: n_myear_leg ! Maximum number of iteration of the PEM before orb changes too much
 
 !--------------------------------------------------------
 ! Local variables
 !--------------------------------------------------------
-integer :: Year                                           ! Year of the simulation
+real    :: Year                                           ! Year of the simulation
 real    :: Lsp                                            ! Ls perihelion
 real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
 real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
 real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
-integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
+real    :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
 integer :: ilask                                          ! Loop variable
 real    :: xa, xb, ya, yb                                 ! Variables for interpolation
@@ -113,7 +113,7 @@
 ! If we do not want some orb parameter to change, they should not be a stopping criterion,
 ! So the number of iteration corresponding is set to maximum
-max_obl_iter = 1000000
-max_ecc_iter = 1000000
-max_lsp_iter = 1000000
+max_obl_iter = 1000000.
+max_ecc_iter = 1000000.
+max_lsp_iter = 1000000.
 
 ! The maximum reachable year is found by interpolation according to the maximum admissible change of orbital parameters
@@ -136,8 +136,8 @@
         yb = obllask(ilask - 1)
         if (yb < min_obl) then ! The minimum accepted is overtaken
-            max_obl_iter = floor((min_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_obl_iter = (min_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_obl = .true.
         else if (max_obl < yb) then ! The maximum accepted is overtaken
-            max_obl_iter = floor((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_obl_iter = (max_obl - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_obl = .true.
         endif
@@ -149,8 +149,8 @@
         yb = ecclask(ilask - 1)
         if (yb < min_ecc) then ! The minimum accepted is overtaken
-            max_ecc_iter = floor((min_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_ecc_iter = (min_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_ecc = .true.
         else if (max_ecc < yb) then ! The maximum accepted is overtaken
-            max_ecc_iter = floor((max_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_ecc_iter = (max_ecc - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_ecc = .true.
         endif
@@ -187,8 +187,8 @@
         endif
         if (yb < min_lsp) then ! The minimum accepted is overtaken
-            max_lsp_iter = floor((min_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_lsp_iter = (min_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_lsp = .true.
         else if (max_lsp < yb) then ! The maximum accepted is overtaken
-            max_lsp_iter = floor((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            max_lsp_iter = (max_lsp - ya)*(xb - xa)/(yb - ya) + xa - Year
             found_lsp = .true.
         endif
@@ -206,7 +206,7 @@
 write(*,*) 'Max. number of iterations for the Lsp parameter  =', max_lsp_iter
 
-year_iter_max = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
-year_iter_max = max(year_iter_max,1)
-write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
+n_myear_leg = min(max_obl_iter,max_ecc_iter,max_lsp_iter)
+if (n_myear_leg < 1.) n_myear_leg = 1. ! n_myear_leg should be at least equal to 1
+write(*,*) 'So the max. number of iterations for the orbital criteria =', n_myear_leg
 
 END SUBROUTINE orbit_param_criterion
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3498)
@@ -51,5 +51,5 @@
                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
-use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
+use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
 use orbit_param_criterion_mod,  only: orbit_param_criterion
 use recomp_orb_param_mod,       only: recomp_orb_param
@@ -70,5 +70,5 @@
 use writediagpem_mod,           only: writediagpem, writediagsoilpem
 use co2condens_mod,             only: CO2cond_ps
-use layering_mod,               only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
+use layering_mod,               only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
 
 #ifndef CPP_STD
@@ -195,11 +195,11 @@
 
 ! Variables for slopes
-real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
-real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
-real, dimension(:,:),   allocatable :: tend_h2o_ice       ! physical point x slope field: Tendency of evolution of perennial h2o ice
-real, dimension(:,:),   allocatable :: flag_co2flow       ! (ngrid,nslope): Flag where there is a CO2 glacier flow
-real, dimension(:),     allocatable :: flag_co2flow_mesh  ! (ngrid)       : Flag where there is a CO2 glacier flow
-real, dimension(:,:),   allocatable :: flag_h2oflow       ! (ngrid,nslope): Flag where there is a H2O glacier flow
-real, dimension(:),     allocatable :: flag_h2oflow_mesh  ! (ngrid)       : Flag where there is a H2O glacier flow
+real, dimension(:,:),   allocatable :: d_co2ice          ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
+real, dimension(:,:),   allocatable :: d_co2ice_ini      ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
+real, dimension(:,:),   allocatable :: d_h2oice          ! physical point x slope field: Tendency of evolution of perennial h2o ice
+real, dimension(:,:),   allocatable :: flag_co2flow      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
+real, dimension(:),     allocatable :: flag_co2flow_mesh ! (ngrid)       : Flag where there is a CO2 glacier flow
+real, dimension(:,:),   allocatable :: flag_h2oflow      ! (ngrid,nslope): Flag where there is a H2O glacier flow
+real, dimension(:),     allocatable :: flag_h2oflow_mesh ! (ngrid)       : Flag where there is a H2O glacier flow
 
 ! Variables for surface and soil
@@ -230,8 +230,8 @@
 ! Some variables for the PEM run
 real, parameter :: year_step = 1   ! Timestep for the pem
-integer         :: year_iter       ! Number of iteration
-integer         :: year_iter_max   ! Maximum number of iterations before stopping
-integer         :: i_myear         ! Global number of Martian years of the chained simulations
-integer         :: n_myear         ! Maximum number of Martian years of the chained simulations
+real            :: i_myear_leg       ! Number of iteration
+real            :: n_myear_leg     ! Maximum number of iterations before stopping
+real            :: i_myear         ! Global number of Martian years of the chained simulations
+real            :: n_myear         ! Maximum number of Martian years of the chained simulations
 real            :: timestep        ! Timestep [s]
 character(20)   :: job_id          ! Job id provided as argument passed on the command line when the program was invoked
@@ -604,11 +604,11 @@
 !    I_f Compute tendencies
 !------------------------
-allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
-allocate(tend_co2_ice_ini(ngrid,nslope))
+allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
+allocate(d_co2ice_ini(ngrid,nslope))
 
 ! Compute the tendencies of the evolution of ice over the years
-call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice)
-call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice)
-tend_co2_ice_ini = tend_co2_ice
+call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
+call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
+d_co2ice_ini = d_co2ice
 deallocate(min_co2_ice,min_h2o_ice)
 
@@ -646,6 +646,6 @@
 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
               icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
-              ps_timeseries,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,               &
-              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,          &
+              ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
+              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,         &
               delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
 
@@ -662,9 +662,9 @@
     do islope = 1,nslope
         if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
-        if (tend_co2_ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
+        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
             ini_co2ice_sublim(i,islope) = .true.
             co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
         endif
-        if (tend_h2o_ice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
+        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
             ini_h2oice_sublim(i,islope) = .true.
             h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
@@ -712,6 +712,6 @@
 #endif
 
-year_iter_max = Max_iter_pem
-if (evol_orbit_pem) call orbit_param_criterion(i_myear,year_iter_max)
+n_myear_leg = Max_iter_pem
+if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
 
 !-------------------------- END INITIALIZATION -------------------------
@@ -722,5 +722,5 @@
 !    II_a Update pressure, ice and tracers
 !------------------------
-year_iter = 0
+i_myear_leg = 0
 stopPEM = 0
 if (layering_algo) then
@@ -737,10 +737,10 @@
 endif
 
-do while (year_iter < year_iter_max .and. i_myear < n_myear)
+do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
 ! II.a.1. Compute updated global pressure
     write(*,*) "Recomputing the new pressure..."
     do i = 1,ngrid
         do islope = 1,nslope
-            global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
+            global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
         enddo
     enddo
@@ -817,10 +817,10 @@
 !    II_b Evolution of ice
 !------------------------
-    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
-    call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
+    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
+    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
     if (layering_algo) then
         do islope = 1,nslope
             do ig = 1,ngrid
-                call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
+                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
             enddo
         enddo
@@ -832,5 +832,5 @@
 !------------------------
     if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
-                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
+                                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
     if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
 
@@ -878,5 +878,5 @@
                     emis(ig,islope) = emissiv
                 endif
-            else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
+            else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
                 ave = 0.
                 do t = 1,timelen
@@ -945,9 +945,9 @@
 
 ! II_d.4 Update the soil thermal properties
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
 
 ! II_d.5 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
-            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,                   &
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                   &
                                      h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
                                      h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
@@ -986,6 +986,6 @@
         call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
         call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
-        call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
-        call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
+        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
+        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
         call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
         call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
@@ -1012,5 +1012,5 @@
 !    II_f Update the tendencies
 !------------------------
-    call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
+    call recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
                                global_avg_press_PCM,global_avg_press_new)
 
@@ -1023,7 +1023,7 @@
     call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
     call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
-    year_iter = year_iter + dt_pem
-    i_myear = i_myear + dt_pem
-    if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5
+    i_myear_leg = i_myear_leg + dt
+    i_myear = i_myear + dt
+    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
     if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
     call system_clock(c2)
@@ -1053,5 +1053,5 @@
     else
         write(*,*) "The PEM can continue!"
-        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
+        write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
         write(*,*) "Number of simulated Martian years: i_myear =", i_myear
     endif
@@ -1151,5 +1151,5 @@
 enddo
 
-if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
+if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
 
 !------------------------
@@ -1213,10 +1213,10 @@
 write(*,*) "restartpem.nc has been written"
 
-call info_PEM(year_iter,stopPEM,i_myear,n_myear)
-
-write(*,*) "The PEM has run for", year_iter, "Martian years."
+call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
+
+write(*,*) "The PEM has run for", i_myear_leg, "Martian years."
 write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
 write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
-write(*,*) "LL & RV & JBC: so far, so good!"
+write(*,*) "PEM: so far, so good!"
 
 if (layering_algo) then
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3498)
@@ -8,5 +8,5 @@
 
 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness,      &
-                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, &
+                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice, &
                     global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                         &
                     m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
@@ -49,6 +49,6 @@
 real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
 real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
-real, dimension(ngrid,nslope),  intent(in) :: tend_h2o_ice        ! tendencies on h2o ice
-real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! tendencies on co2 ice
+real, dimension(ngrid,nslope),  intent(in) :: d_h2oice        ! tendencies on h2o ice
+real, dimension(ngrid,nslope),  intent(in) :: d_co2ice        ! tendencies on co2 ice
 real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! outputs
@@ -105,5 +105,5 @@
 inquire(file = filename,exist = startpem_file)
 
-write(*,*)'Is start PEM?',startpem_file
+write(*,*)'Is there a "startpem" file?',startpem_file
 
 !0.2 Set to default values
@@ -315,5 +315,5 @@
                 write(*,*)'will reconstruct the values of the ice table given the current state'
                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -327,5 +327,5 @@
                 write(*,*)'will reconstruct the values of the ice table given the current state'
                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -360,5 +360,5 @@
             enddo
 
-            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
+            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
                                         m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
 
@@ -504,5 +504,5 @@
         if (icetable_equilibrium) then
             call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness)
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -512,5 +512,5 @@
             ice_porefilling = 0.
             ice_table_depth = 0.
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -524,5 +524,5 @@
             m_co2_regolith_phys = 0.
             m_h2o_regolith_phys = 0.
-            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
+            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
                                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
             deltam_co2_regolith_phys = 0.
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3498)
@@ -73,5 +73,6 @@
 
 character(*),                            intent(in) :: filename
-integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope, i_myear
+integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope
+real,                                    intent(in) :: i_myear
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3498)
@@ -13,5 +13,5 @@
 !=======================================================================
 
-SUBROUTINE recomp_orb_param(i_myear,year_iter)
+SUBROUTINE recomp_orb_param(i_myear,i_myear_leg)
 
 use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp
@@ -30,6 +30,6 @@
 ! Input Variables
 !--------------------------------------------------------
-integer, intent(in) :: i_myear   ! Number of simulated Martian years
-integer, intent(in) :: year_iter ! Number of iterations of the PEM
+real, intent(in) :: i_myear   ! Number of simulated Martian years
+real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
 
 !--------------------------------------------------------
@@ -40,5 +40,5 @@
 ! Local variables
 !--------------------------------------------------------
-integer :: Year           ! Year of the simulation
+real    :: Year           ! Year of the simulation
 real    :: Lsp            ! Ls perihelion
 integer :: ilask          ! Loop variable
@@ -61,5 +61,5 @@
 Lsp = lsperi*180./pi         ! We convert in degree
 
-write(*,*) 'recomp_orb_param, Old year =', Year - year_iter
+write(*,*) 'recomp_orb_param, Old year =', Year - i_myear_leg
 write(*,*) 'recomp_orb_param, Old obl. =', obliquit
 write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3498)
@@ -7,5 +7,5 @@
 !=======================================================================
 
-SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
+SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_phys_ini,co2ice_slope,emissivity_slope, &
                                  vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
 
@@ -29,9 +29,9 @@
 real,                           intent(in) :: global_avg_press_PCM        ! global averaged pressure at previous timestep
 real,                           intent(in) :: global_avg_press_new        ! global averaged pressure at current timestep
-real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field: Evolution of perennial ice over one year
+real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_phys_ini ! physical point field: Evolution of perennial ice over one year
 real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope (kg/m^2)
 real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
 !   OUTPUT
-real, dimension(ngrid,nslope), intent(inout) ::  tendencies_co2_ice_phys ! physical point field: Evolution of perennial ice over one year
+real, dimension(ngrid,nslope), intent(inout) ::  d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
 
 !   local:
@@ -47,5 +47,5 @@
         coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
         ave = 0.
-        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
+        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
             do t=1,timelen
                 ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
@@ -53,5 +53,5 @@
             enddo
             if (ave < 1e-4) ave = 0.
-            tendencies_co2_ice_phys(i,islope) = tendencies_co2_ice_phys_ini(i,islope) - coef*ave/timelen
+            d_co2ice_phys(i,islope) = d_co2ice_phys_ini(i,islope) - coef*ave/timelen
         endif
     enddo
Index: trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3498)
@@ -3,11 +3,11 @@
 implicit none
 
-integer :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
-integer :: dt_pem         ! Time step used by the PEM in Planetary years
+real    :: year_bp_ini    ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
+real    :: dt             ! Time step used by the PEM in Planetary years
 real    :: convert_years  ! Conversion ratio from Planetary years to Earth years
 real    :: h2o_ice_crit   ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
 real    :: co2_ice_crit   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
 real    :: ps_criterion   ! Percentage of change of averaged surface pressure before stopping the PEM
-integer :: Max_iter_pem   ! Maximal number of iteration when converging to a steady state, read in evol.def
+integer :: Max_iter_pem   ! Maximal number of iterations when converging to a steady state, read in evol.def
 logical :: evol_orbit_pem ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
 logical :: var_obl        ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
Index: trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3497)
+++ trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90	(revision 3498)
@@ -51,5 +51,5 @@
 use mod_phys_lmdz_para, only: is_parallel, is_mpi_root, is_master, gather
 use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
-use time_evol_mod,      only: ecritpem, dt_pem
+use time_evol_mod,      only: ecritpem, dt
 
 implicit none
@@ -255,5 +255,5 @@
 ! to writediagpem
 !------------------------------------------------------------------------
-if (nom == firstnom) zitau = zitau + dt_pem
+if (nom == firstnom) zitau = zitau + 1
 
 !--------------------------------------------------------
@@ -600,5 +600,5 @@
 use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat
 use mod_grid_phy_lmdz,  only: grid_type, unstructured
-use time_evol_mod,      only: ecritpem, dt_pem
+use time_evol_mod,      only: ecritpem, dt
 use iniwritesoil_mod,   only: iniwritesoil
 
@@ -747,5 +747,5 @@
 if (name.eq.firstname) then
   ! if we run across 'firstname', then it is a new time step
-  zitau = zitau + dt_pem
+  zitau = zitau + dt
 endif
 
