Index: unk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90	(revision 3038)
+++ 	(revision )
@@ -1,51 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,&
-     min_ice_Y2,tendencies_ice)
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!  Compute the tendencies of the evolution of water ice over the years
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-
-     INTEGER, intent(in) :: ngrid, nslope                           ! # of grid points along longitude/latitude/ total
-     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
-     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
-
-!   OUTPUT
-     REAL, intent(out) , dimension(ngrid,nslope)   :: tendencies_ice            ! physical point field : difference between the minima = evolution of perenial ice
-
-!   local:
-!   ------
-     INTEGER :: ig,islope                                                           ! loop variable
-
-!=======================================================================
-
-!  We compute the difference
-
-  DO ig=1,ngrid
-    DO islope = 1, nslope
-      tendencies_ice(ig,islope)=min_ice_Y2(ig,islope)-min_ice_Y1(ig,islope)
-    enddo
-  ENDDO
-
-!  If the difference is too small; there is no evolution
-  DO ig=1,ngrid
-    DO islope = 1, nslope
-      if(abs(tendencies_ice(ig,islope)).LT.1.0E-10) then
-        tendencies_ice(ig,islope)=0.
-      endif
-    enddo
-  ENDDO
-
-END SUBROUTINE compute_tendencies_slope
-
Index: /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope.F90	(revision 3039)
+++ /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope.F90	(revision 3039)
@@ -0,0 +1,51 @@
+!
+! $Id $
+!
+SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,&
+     min_ice_Y2,tendencies_ice)
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+!  Compute the tendencies of the evolution of water ice over the years
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+
+!   INPUT
+
+     INTEGER, intent(in) :: ngrid, nslope                           ! # of grid points along longitude/latitude/ total
+     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
+     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
+
+!   OUTPUT
+     REAL, intent(out) , dimension(ngrid,nslope)   :: tendencies_ice            ! physical point field : difference between the minima = evolution of perenial ice
+
+!   local:
+!   ------
+     INTEGER :: ig,islope                                                           ! loop variable
+
+!=======================================================================
+
+!  We compute the difference
+
+  DO ig=1,ngrid
+    DO islope = 1, nslope
+      tendencies_ice(ig,islope)=min_ice_Y2(ig,islope)-min_ice_Y1(ig,islope)
+    enddo
+  ENDDO
+
+!  If the difference is too small; there is no evolution
+  DO ig=1,ngrid
+    DO islope = 1, nslope
+      if(abs(tendencies_ice(ig,islope)).LT.1.0E-10) then
+        tendencies_ice(ig,islope)=0.
+      endif
+    enddo
+  ENDDO
+
+END SUBROUTINE compute_tendencies_slope
+
Index: /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3039)
@@ -5,5 +5,5 @@
 ! Purpose: Read the parameter for a PEM run in run_pem.def
 !
-! Author: RV 
+! Author: RV, JBC
 !=======================================================================
 
@@ -12,5 +12,5 @@
 CONTAINS
 
-  SUBROUTINE conf_pem
+  SUBROUTINE conf_pem(i_myear,n_myears)
 
 #ifdef CPP_IOIPSL
@@ -21,119 +21,131 @@
 #endif
   
-  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, &
-                Max_iter_pem, evol_orbit_pem, var_obl, var_ex, var_lsp
-  USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom,depth_breccia,depth_bedrock,reg_thprop_dependp 
-  USE adsorption_mod,only: adsorption_pem
-  USE glaciers_mod, only: co2glaciersflow,h2oglaciersflow
-  use ice_table_mod, only: icetable_equilibrium, icetable_dynamic 
+  use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
+                            evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
+  use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp 
+  use adsorption_mod, only: adsorption_pem
+  use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
+  use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic 
 
-  CHARACTER(len=20),parameter :: modname ='conf_pem' 
+  integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
+
+  character(len=20), parameter :: modname ='conf_pem'
+  integer                      :: ierr
+  integer                      :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
 
 !PEM parameters
 
-! #---------- ORBITAL parameters --------------#
+! #---------- Martian years parameters from launching script
+open(100,file = 'tmp_PEMyears.txt',status = 'old',form = 'formatted',iostat = ierr)
+if (ierr /= 0) then
+    write(*,*) 'Cannot find required file "tmp_PEMyears.txt"!'
+    write(*,*) 'It should be created by the launching script...'
+    stop
+else
+    read(100,*) i_myear, n_myears, convert_years
+endif
+close(100)
 
+! #---------- Orbital parameters
   evol_orbit_pem=.false.
-  CALL getin('evol_orbit_pem', evol_orbit_pem)
+  call getin('evol_orbit_pem',evol_orbit_pem)
 
-  year_bp_ini=0.
-  CALL getin('year_bp_ini', year_bp_ini)
+  year_bp_ini = 0.
+  call getin('year_earth_bp_ini',year_earth_bp_ini)
+  year_bp_ini = year_earth_bp_ini/convert_years
 
   var_obl = .true.
-  CALL getin('var_obl',var_obl)
-  print*,'Does obliquity vary ?',var_obl
+  call getin('var_obl',var_obl)
+  write(*,*) 'Does obliquity vary ?',var_obl
 
-  var_ex = .true.
-  CALL getin('var_ex',var_ex)
-  print*,'Does excentricity vary ?',var_ex
+  var_ecc = .true.
+  call getin('var_ecc',var_ecc)
+  write(*,*) 'Does excentricity vary ?',var_ecc
 
   var_lsp = .true.
-  CALL getin('var_lsp',var_lsp)
-  print*,'Does Ls peri vary ?',var_lsp
+  call getin('var_lsp',var_lsp)
+  write(*,*) 'Does Ls peri vary ?',var_lsp
 
-! #---------- Stopping criterion parameters --------------#
-
-  Max_iter_pem=99999999
-  CALL getin('Max_iter_pem', Max_iter_pem)
+! #---------- Stopping criteria parameters
+  Max_iter_pem=100000000
+  call getin('Max_iter_pem', Max_iter_pem)
 
   water_ice_criterion=0.2
-  CALL getin('water_ice_criterion', water_ice_criterion)
+  call getin('water_ice_criterion', water_ice_criterion)
 
   co2_ice_criterion=0.2
-  CALL getin('co2_ice_criterion', co2_ice_criterion)
+  call getin('co2_ice_criterion', co2_ice_criterion)
 
   ps_criterion = 0.15
-  CALL getin('ps_criterion',ps_criterion)
+  call getin('ps_criterion',ps_criterion)
 
   dt_pem=1
-  CALL getin('dt_pem', dt_pem)
+  call getin('dt_pem', dt_pem)
 
-! #---------- Subsurface parameters --------------#
-
+! #---------- Subsurface parameters
   soil_pem=.true.
-  CALL getin('soil_pem', soil_pem)
+  call getin('soil_pem', soil_pem)
 
   adsorption_pem = .true.
-  CALL getin('adsorption_pem',adsorption_pem)
+  call getin('adsorption_pem',adsorption_pem)
 
   co2glaciersflow = .true.
-  CALL getin('co2glaciersflow', co2glaciersflow)
+  call getin('co2glaciersflow', co2glaciersflow)
 
   h2oglaciersflow = .true.
-  CALL getin('h2oglaciersflow', h2oglaciersflow)
+  call getin('h2oglaciersflow', h2oglaciersflow)
 
   reg_thprop_dependp = .false.
-  CALL getin('reg_thprop_dependp',reg_thprop_dependp)
-  print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
+  call getin('reg_thprop_dependp',reg_thprop_dependp)
+  write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
 
-! #---------- Layering parameters --------------#
-
+! #---------- Layering parameters
   fluxgeo = 0.
-  CALL getin('fluxgeo',fluxgeo)
-  print*,'Flux Geothermal is set to',fluxgeo
+  call getin('fluxgeo',fluxgeo)
+  write(*,*) 'Flux Geothermal is set to',fluxgeo
    
   depth_breccia   = 10.
-  CALL getin('depth_breccia',depth_breccia)
-  print*,'Depth of breccia is set to',depth_breccia
+  call getin('depth_breccia',depth_breccia)
+  write(*,*) 'Depth of breccia is set to',depth_breccia
 
   depth_bedrock   = 1000.
-  CALL getin('depth_bedrock',depth_bedrock)
-  print*,'Depth of bedrock is set to',depth_bedrock
+  call getin('depth_bedrock',depth_bedrock)
+  write(*,*) 'Depth of bedrock is set to',depth_bedrock
 
    icetable_equilibrium = .true.
-   CALL getin('icetable_equilibrium',icetable_equilibrium)
-   print*, 'Do we compute the ice table at equilibrium?', icetable_equilibrium
+   call getin('icetable_equilibrium',icetable_equilibrium)
+   write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
 
    icetable_dynamic = .false.
-   CALL getin('icetable_dynamic',icetable_dynamic)
-   print*, 'Do we compute the ice table with the dynamic method?', icetable_dynamic
+   call getin('icetable_dynamic',icetable_dynamic)
+   write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
   if ((.not.soil_pem).and.((icetable_equilibrium).or.(icetable_dynamic))) then
-       print*,'Ice table  must be used when soil_pem = T'
+       write(*,*) 'Ice table  must be used when soil_pem = T'
        call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
   endif
 
   if ((.not.soil_pem).and.adsorption_pem) then
-       print*,'Adsorption must be used when soil_pem = T'
+       write(*,*) 'Adsorption must be used when soil_pem = T'
        call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
   endif
   
   if ((.not.soil_pem).and.(fluxgeo.gt.0.)) then
-       print*,'Soil is not activated but Flux Geo > 0.'
+       write(*,*) 'Soil is not activated but Flux Geo > 0.'
        call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
   endif
  
   if ((.not.soil_pem).and.reg_thprop_dependp) then
-     print*,'Regolith properties vary with Ps only when soil is set to true'
+     write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
      call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
   endif
 
-  if (evol_orbit_pem.and.year_bp_ini.eq.0.) then
-     print*,'You want to follow the file ob_ex_lsp.asc for changing orb parameters,'
-     print*,'but you did not specify from which year to start.'
+  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,'
+     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)
   endif 
 
   water_reservoir_nom = 1e4
-  CALL getin('water_reservoir_nom',water_reservoir_nom)
+  call getin('water_reservoir_nom',water_reservoir_nom)
 
   END SUBROUTINE conf_pem
Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 3039)
@@ -15,6 +15,6 @@
 SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
 
-  USE temps_mod_evol, ONLY: water_ice_criterion
-  use comslope_mod, ONLY: subslope_dist,nslope
+  use time_evol_mod, only: water_ice_criterion
+  use comslope_mod,  only: subslope_dist,nslope
 
       IMPLICIT NONE
@@ -79,6 +79,6 @@
 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)
 
-  USE temps_mod_evol, ONLY: co2_ice_criterion,ps_criterion
-  use comslope_mod, ONLY: subslope_dist
+  use time_evol_mod, only: co2_ice_criterion,ps_criterion
+  use comslope_mod,  only: subslope_dist
 
       IMPLICIT NONE
Index: /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90	(revision 3039)
@@ -8,5 +8,5 @@
                              iim_input,jjm_input,ngrid,cell_area,nslope)
 
-  USE temps_mod_evol, ONLY: dt_pem
+  use time_evol_mod, only: dt_pem
 
       IMPLICIT NONE
Index: /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90	(revision 3039)
@@ -7,8 +7,8 @@
 SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING)
 
-  use temps_mod_evol, only: dt_pem
-  use comslope_mod, only: subslope_dist,def_slope_mean
+  use time_evol_mod,          only: dt_pem
+  use comslope_mod,           only: subslope_dist, def_slope_mean
   use criterion_pem_stop_mod, only: criterion_waterice_stop
-  use comconst_mod,only: pi
+  use comconst_mod,           only: pi
 
   IMPLICIT NONE
Index: /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90	(revision 3039)
@@ -1,25 +1,33 @@
-SUBROUTINE info_run_PEM(year_iter, criterion_stop)
+SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
 
 !=======================================================================
 !
 ! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
+!          Update the file tmp_PEMyears.txt to count the number of simulated Martian years
 !
-! Author: RV 
+! Author: RV, JBC
 !=======================================================================
 
-   IMPLICIT NONE
+  use time_evol_mod, only: convert_years
 
-  INTEGER, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
+  IMPLICIT NONE
 
-  logical :: exist
+  integer, intent(in) :: year_iter, criterion_stop ! # 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
 
-  inquire(file='info_run_PEM.txt', exist=exist)
-  if (exist) then
-    open(12, file='info_run_PEM.txt', status="old", position="append", action="write")
-  else
-    open(12, file='info_run_PEM.txt', status="new", action="write")
-  end if
-  write(12, *) year_iter, criterion_stop
-  close(12)
+  logical :: ok
+
+inquire(file = 'info_run_PEM.txt', exist = ok)
+if (ok) then
+    open(12,file = 'info_run_PEM.txt',status = "old",position = "append",action = "write")
+else
+    open(12,file = 'info_run_PEM.txt',status = "new",action = "write")
+endif
+write(12,*) year_iter, i_myear, criterion_stop
+close(12)
+
+open(100,file = 'tmp_PEMyears.txt',status = 'replace')
+write(100,*) i_myear, n_myear, convert_years
+close(100)
 
 END SUBROUTINE info_run_PEM
Index: /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3039)
@@ -4,5 +4,5 @@
 !!! Purpose: Define parameters from Laskar et al., 2004 evolution
 !!! 
-!!! Author: RV
+!!! Author: RV, JBC
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -10,20 +10,30 @@
 implicit none
 
-  real,save,allocatable :: yearlask(:)          ! year before present from Laskar et al. Tab
-  real,save,allocatable :: oblask(:)            ! obliquity    [deg]
-  real,save,allocatable :: exlask(:)            ! excentricity [deg]
-  real,save,allocatable :: lsplask(:)           ! ls perihelie [deg]
-  integer, save :: last_ilask                   ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
+  real,save,allocatable :: yearlask(:) ! year before present from Laskar et al. Tab
+  real,save,allocatable :: obllask(:)  ! obliquity    [deg]
+  real,save,allocatable :: ecclask(:)  ! excentricity [deg]
+  real,save,allocatable :: lsplask(:)  ! ls perihelie [deg]
+  integer, save         :: last_ilask  ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
 
 contains
 
-  subroutine ini_lask_param_mod(nlask)
+  subroutine ini_lask_param_mod
   
   implicit none
-  integer,intent(in) :: nlask ! number of line
-  
+
+  integer :: nlask ! number of lines in Laskar files
+  integer :: ierr
+
+  nlask = 0
+  open(1,file = 'obl_ecc_lsp.asc')
+  do
+      read(1,*,iostat = ierr)
+      if (ierr /= 0) exit
+      nlask = nlask + 1
+  enddo
+  close(1)
   allocate(yearlask(nlask))
-  allocate(oblask(nlask))
-  allocate(exlask(nlask))
+  allocate(obllask(nlask))
+  allocate(ecclask(nlask))
   allocate(lsplask(nlask))
   
@@ -36,6 +46,6 @@
 
   if (allocated(yearlask)) deallocate(yearlask)
-  if (allocated(oblask)) deallocate(oblask)
-  if (allocated(exlask)) deallocate(exlask)
+  if (allocated(obllask)) deallocate(obllask)
+  if (allocated(ecclask)) deallocate(ecclask)
   if (allocated(lsplask)) deallocate(lsplask)
 
Index: /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3039)
@@ -1,3 +1,3 @@
-      MODULE orbit_param_criterion_mod
+MODULE orbit_param_criterion_mod
 
 !=======================================================================
@@ -9,226 +9,214 @@
 !=======================================================================
 
-      IMPLICIT NONE
-
-      CONTAINS
-
-      SUBROUTINE orbit_param_criterion(year_iter_max)
+IMPLICIT NONE
+
+CONTAINS
+
+SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
+
 #ifdef CPP_IOIPSL
-  use IOIPSL, only: getin
+    use IOIPSL,          only: getin
 #else
-  ! if not using IOIPSL, we still need to use (a local version of) getin
-  use ioipsl_getincom, only: getin
+    ! if not using IOIPSL, we still need to use (a local version of) getin
+    use ioipsl_getincom, only: getin
 #endif
-
-      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
 #ifndef CPP_STD
-      USE planete_h, ONLY: e_elips, obliquit, timeperi
+    use planete_h,   only: e_elips, obliquit, timeperi
 #else
-      use planete_mod, only: e_elips, obliquit, timeperi
+    use planete_mod, only: e_elips, obliquit, timeperi
 #endif
-      USE comconst_mod, ONLY: pi
-      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
-   				ini_lask_param_mod,last_ilask
-
-      IMPLICIT NONE
+  use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
+  use comconst_mod,   only: pi
+  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
+
+  IMPLICIT NONE
 
 !--------------------------------------------------------
 ! Input Variables
 !--------------------------------------------------------
+  integer, 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
+  integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
 
 !--------------------------------------------------------
 ! Local variables 
 !--------------------------------------------------------
-
-      real :: Year	                ! Year of the simulation
-      real :: Lsp	                ! Year of the simulation
-      integer nlask, ilask, iylask      ! Loop variables
-      parameter (nlask = 20001)         ! Number of line in Laskar file
-
-      real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible
-
-      real max_obl,max_ex,max_lsp       ! Maximum value of orbit param given the acceptable percentage 
-      real min_obl,min_ex,min_lsp       ! Maximum value of orbit param given the acceptable percentage
-      real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value
-      real xa,xb,ya,yb,yc
-
-      ! **********************************************************************
-      ! 0. Initializations
-      ! **********************************************************************
-
-          Year = year_bp_ini + year_PEM      ! We compute the current year
-          Lsp = 360. - timeperi*360./(2.*pi) ! We convert in degree
-
-	  call ini_lask_param_mod(nlask) ! Allocation of variables
-
-          print*, "orbit_param_criterion, Year in pem.def=", year_bp_ini
-          print*, "orbit_param_criterion, Year in the startpem.nc =", year_PEM
-          print*, "orbit_param_criterion, Current year=", Year
-          print*, "orbit_param_criterion, Current obl=", obliquit
-          print*, "orbit_param_criterion, Current ex=", e_elips
-          print*, "orbit_param_criterion, Current lsp=", Lsp
+  integer :: 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
+  integer :: ilask, iylask                                  ! Loop variables
+  real    :: xa, xb, ya, yb, yc_obl, yc_ecc, yc_lsp         ! Variables for interpolation
+
+! **********************************************************************
+! 0. Initializations
+! **********************************************************************
+Year = year_bp_ini + i_myear        ! We compute the current year
+Lsp = modulo(-timeperi*180./pi,360.) ! We convert in degree
+
+call ini_lask_param_mod ! Allocation of variables
+
+write(*,*) 'orbit_param_criterion, Current year =', Year
+write(*,*) 'orbit_param_criterion, Current obl. =', obliquit
+write(*,*) 'orbit_param_criterion, Current ecc. =', e_elips
+write(*,*) 'orbit_param_criterion, Current Lsp  =', Lsp
 
 ! We read the file
-
-          open(73,file='ob_ex_lsp.asc')
-          do ilask = 1,nlask
-            read(73,*) yearlask(ilask), oblask(ilask), exlask(ilask), lsplask(ilask)
-            yearlask(ilask) = yearlask(ilask)*1000.
-            if (yearlask(ilask) > Year) last_ilask = ilask + 1
-          enddo
-          close(73)
-
-       print*, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask
-
-!Constant max change case 
-
-        max_change_obl = 0.5
-        max_change_ex = 0.1
-        max_change_lsp = 40.
-
-        CALL getin('max_change_obl', max_change_obl)
-	max_obl = obliquit + max_change_obl
-	min_obl = obliquit - max_change_obl
-        print*, "Maximum obliquity accepted=", max_obl
-        print*, "Minimum obliquity accepted=", min_obl
-
-        CALL getin('max_change_ex', max_change_ex)
-	max_ex = e_elips + max_change_ex
-	min_ex = e_elips - max_change_ex
-        print*, "Maximum excentricity accepted=", max_ex
-        print*, "Minimum excentricity accepted=", min_ex
-
-        CALL getin('max_change_lsp', max_change_lsp)
-	max_lsp = Lsp + max_change_lsp
-	min_lsp = Lsp - max_change_lsp
-        print*, "Maximum lsp accepted=", max_lsp
-        print*, "Minimum lsp accepted=", min_lsp
-
-!End Constant max change case 
+open(73,file = 'obl_ecc_lsp.asc')
+do ilask = 1,size(yearlask,1)
+    read(73,*) yearlask(ilask), obllask(ilask), ecclask(ilask), lsplask(ilask)
+    yearlask(ilask) = yearlask(ilask)*1000./convert_years ! Conversion from Laskar's kilo Earth years to PEM Martian years
+    if (yearlask(ilask) > Year) last_ilask = ilask + 1
+enddo
+close(73)
+write(*,*) 'Corresponding line in the obl_ecc_lsp.asc file =', last_ilask
+
+! Constant max change case 
+max_change_obl = 0.6
+max_change_ecc = 2.8e-3
+max_change_lsp = 18.
+
+call getin('max_change_obl', max_change_obl)
+max_obl = obliquit + max_change_obl
+min_obl = obliquit - max_change_obl
+write(*,*) 'Maximum obliquity accepted =', max_obl
+write(*,*) 'Minimum obliquity accepted =', min_obl
+
+call getin('max_change_ecc', max_change_ecc)
+max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ex < 1.
+min_ecc = max(e_elips - max_change_ecc,0.) ! ex >= 0.
+write(*,*) 'Maximum eccentricity accepted =', max_ecc
+write(*,*) 'Minimum eccentricity accepted =', min_ecc
+
+call getin('max_change_lsp', max_change_lsp)
+max_lsp = Lsp + max_change_lsp
+min_lsp = Lsp - max_change_lsp
+write(*,*) 'Maximum Lsp accepted =', max_lsp
+write(*,*) 'Minimum Lsp accepted =', min_lsp
+! End Constant max change case 
 
 ! 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 = 999999
-        max_ex_iter = 999999
-        max_lsp_iter = 999999
+max_obl_iter = 1000000
+max_ecc_iter = 1000000
+max_lsp_iter = 1000000
 
         ! Tendency of the orbital parameter for the considered year gives the limitation between min and max
-        do ilask = last_ilask,2,-1
-          xa = yearlask(ilask)
-          xb = yearlask(ilask - 1)
-          if (xa <= Year .and. Year < xb) then
-            ! Obliquity
-            if (var_obl) then
-              ya = oblask(ilask)
-              yb = oblask(ilask - 1)
-              if (ya < yb) then ! Increasing -> max is the limitation
-                yc = max_obl
-              else ! Decreasing -> min is the limitation
-                yc = min_obl
-              endif
-            endif
-
-            ! Excentricity
-            if (var_ex) then
-              ya = exlask(ilask)
-              yb = exlask(ilask - 1)
-              if (ya < yb) then ! Increasing -> max is the limitation
-                yc = max_ex
-              else ! Decreasing -> min is the limitation
-                yc = min_ex
-              endif
-            endif
-
-            ! Lsp
-            if (var_lsp) then
-              ya = lsplask(ilask)
-              yb = lsplask(ilask - 1)
-              if (ya < yb) then ! Increasing -> max is the limitation
-                if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
-                  !yb = yb - 360.
-                  yc = min_lsp
-                else
-                  yc = max_lsp
-                endif
-              else ! Decreasing -> min is the limitation
-                if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
-                  !yb = yb + 360.
-                  yc = max_lsp
-                else
-                  yc = min_lsp
-                endif
-              endif
-            endif
-            iylask = ilask
-            exit ! The loop is left as soon as the right interval is found
-          endif
-        enddo
-        if (ilask == 1) then
-          write(*,*) 'The year does not match with Laskar data in the file ob_ex_lsp.asc.'
-          stop
-        endif
-
-        ! Linear interpolation gives the maximum reachable year according to the limitation
+do ilask = last_ilask,2,-1
+    xa = yearlask(ilask)
+    xb = yearlask(ilask - 1)
+    if (xa <= Year .and. Year < xb) then
         ! Obliquity
         if (var_obl) then
-          do ilask = iylask,2,-1
-            ya = oblask(ilask)
-            yb = oblask(ilask - 1)
-            if (ya <= yc .and. yc < yb) then
-              xa = yearlask(ilask)
-              xb = yearlask(ilask - 1)
-              max_obl_iter = int((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
-              exit ! The loop is left as soon as the right interval is found
+            ya = obllask(ilask)
+            yb = obllask(ilask - 1)
+            if (ya < yb) then ! Increasing -> max is the limitation
+                yc_obl = max_obl
+            else ! Decreasing -> min is the limitation
+                yc_obl = min_obl
             endif
-          enddo
-	endif
-        ! Excentricity
-        if (var_ex) then
-          do ilask = iylask,2,-1
-            ya = exlask(ilask)
-            yb = exlask(ilask - 1)
-            if (ya <= yc .and. yc < yb) then
-              xa = yearlask(ilask)
-              xb = yearlask(ilask - 1)
-              max_ex_iter = int((max_ex - ya)*(xb - xa)/(yb - ya) + xa) - Year
-              exit ! The loop is left as soon as the right interval is found
-            endif            
-          enddo
-	endif
+        endif
+
+        ! Eccentricity
+        if (var_ecc) then
+            ya = ecclask(ilask)
+            yb = ecclask(ilask - 1)
+            if (ya < yb) then ! Increasing -> max is the limitation
+                yc_ecc = max_ecc
+            else ! Decreasing -> min is the limitation
+                yc_ecc = min_ecc
+            endif
+        endif
+
         ! Lsp
         if (var_lsp) then
-          do ilask = iylask,2,-1
             ya = lsplask(ilask)
             yb = lsplask(ilask - 1)
-            if (ya <= yc .and. yc < yb) then
-              xa = yearlask(ilask)
-              xb = yearlask(ilask - 1)
-              max_lsp_iter = int((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
-              exit ! The loop is left as soon as the right interval is found
+             if (ya < yb) then ! Increasing -> max is the limitation
+                if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around
+                    yc_lsp = min_lsp
+                else
+                    yc_lsp = max_lsp
+                endif
+            else ! Decreasing -> min is the limitation
+                if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around
+                    yc_lsp = max_lsp
+                else
+                    yc_lsp = min_lsp
+                endif
             endif
-          enddo
-	endif
-
-      print*, "Maximum number of iteration for the obl. parameter=", max_obl_iter
-      print*, "Maximum number of iteration for the ex. parameter=", max_ex_iter
-      print*, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter
-
-      year_iter_max = min(max_obl_iter,max_ex_iter,max_lsp_iter)
-      if (year_iter_max > 0) then
-        print*, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max
-      else
-        write(*,*) 'The max. number of iteration (year) is not compatible with Laskar data in the file ob_ex_lsp.asc.'
-        stop
-      endif
-
-        END SUBROUTINE orbit_param_criterion
+        endif
+        iylask = ilask
+        exit ! The loop is left as soon as the right interval is found
+    endif
+enddo
+
+if (ilask == 1) then
+    write(*,*) 'The year does not match with Laskar data in the file obl_ecc_lsp.asc.'
+    stop
+endif
+
+! Linear interpolation gives the maximum reachable year according to the limitation
+! Obliquity
+if (var_obl) then
+    do ilask = iylask,2,-1
+        ya = obllask(ilask)
+        yb = obllask(ilask - 1)
+        if (min(ya,yb) <= yc_obl .and. yc_obl < max(ya,yb)) then
+            xa = yearlask(ilask)
+            xb = yearlask(ilask - 1)
+            max_obl_iter = floor((yc_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            exit ! The loop is left as soon as the right interval is found
+        endif
+    enddo
+endif
+! Eccentricity
+if (var_ecc) then
+    do ilask = iylask,2,-1
+        ya = ecclask(ilask)
+        yb = ecclask(ilask - 1)
+        if (min(ya,yb) <= yc_ecc .and. yc_ecc < max(ya,yb)) then
+            xa = yearlask(ilask)
+            xb = yearlask(ilask - 1)
+            max_ecc_iter = floor((yc_ecc - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            exit ! The loop is left as soon as the right interval is found
+        endif
+    enddo
+endif
+! Lsp
+if (var_lsp) then
+    do ilask = iylask,2,-1
+        ya = lsplask(ilask)
+        yb = lsplask(ilask - 1)
+        if (min(ya,yb) <= yc_lsp .and. yc_lsp < max(ya,yb)) then
+            xa = yearlask(ilask)
+            xb = yearlask(ilask - 1)
+            max_lsp_iter = floor((yc_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year
+            exit ! The loop is left as soon as the right interval is found
+        endif
+    enddo
+endif
+
+write(*,*) 'Max. number of iterations for the obl. parameter =', max_obl_iter
+write(*,*) 'Max. number of iterations for the ecc. parameter =', max_ecc_iter
+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)
+if (year_iter_max > 0) then
+    write(*,*) 'So the max. number of iterations for the orbital criteria =', year_iter_max
+else
+    write(*,*) 'The max. number of iterations is not compatible with Laskar data in the file obl_ecc_lsp.asc.'
+    stop
+endif
+
+END SUBROUTINE orbit_param_criterion
 
 !********************************************************************************   
-      
-      END MODULE orbit_param_criterion_mod
+
+END MODULE orbit_param_criterion_mod
+
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3039)
@@ -53,5 +53,5 @@
                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
-use temps_mod_evol,             only: dt_pem, evol_orbit_pem, Max_iter_pem
+use time_evol_mod,              only: dt_pem, 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
@@ -59,4 +59,5 @@
                                       ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
+use time_phylmdz_mod,           only: daysec, dtphys, day_end
 
 #ifndef CPP_STD
@@ -79,16 +80,16 @@
                                   PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
                                   ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
-    use aerosol_mod,  only: iniaerosol
-    use planete_mod,  only: apoastr, periastr, year_day, peri_day, obliquit
-    use comcstfi_mod, only: r, mugaz
+    use aerosol_mod,        only: iniaerosol
+    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
+    use comcstfi_mod,       only: r, mugaz
 #endif
 
 #ifndef CPP_1D
-    use iniphysiq_mod,            only: iniphysiq
-    use control_mod,              only: iphysiq, day_step, nsplit_phys
+    use iniphysiq_mod,    only: iniphysiq
+    use control_mod,      only: iphysiq, day_step, nsplit_phys
 #else
-    use time_phylmdz_mod,         only:  iphysiq, day_step
+    use time_phylmdz_mod, only: iphysiq, day_step
 #endif
-use time_phylmdz_mod, only: daysec,dtphys
+
 #ifdef CPP_1D
     use regular_lonlat_mod,       only: init_regular_lonlat
@@ -104,4 +105,5 @@
 include "comgeom.h"
 include "iniprint.h"
+include "callkeys.h"
 
 integer ngridmx
@@ -114,4 +116,5 @@
 integer :: day_ini   ! First day of the simulation
 real    :: pday      ! Physical day
+real    :: ptime     ! Physical time
 real    :: time_phys ! Same as GCM
 real    :: ptimestep ! Same as GCM
@@ -154,7 +157,7 @@
 real                                :: global_ave_press_old  ! constant: Global average pressure of initial/previous time step
 real                                :: global_ave_press_new  ! constant: Global average pressure of current time step
-real, dimension(:,:),   allocatable ::  zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
-real, dimension(:,:),   allocatable ::  zplev_gcm            ! same but retrieved from the gcm [kg/m^2]
-real, dimension(:,:,:), allocatable ::  zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
+real, dimension(:,:),   allocatable :: zplev_new             ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
+real, dimension(:,:),   allocatable :: zplev_gcm             ! same but retrieved from the gcm [kg/m^2]
+real, dimension(:,:,:), allocatable :: zplev_new_timeseries  ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
 real, dimension(:,:,:), allocatable :: zplev_old_timeseries  ! same but with the time series, for oldest time step
 logical                             :: STOPPING_water        ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
@@ -185,8 +188,8 @@
 real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
 real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical pointx slope  field : Tendency of evolution of perenial 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 :: 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
@@ -194,6 +197,6 @@
 real, allocatable :: tsoil_ave(:,:,:)            ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
 real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
-real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
+real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
+real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
 real, allocatable :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
 real, allocatable :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
@@ -216,4 +219,6 @@
 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            :: timestep      ! timestep [s]
 real            :: watercap_sum  ! total mass of water cap [kg/m^2] 
@@ -247,5 +252,5 @@
 
 ! Loop variables
-integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil,icap
+integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap
 
 ! Parallel variables
@@ -318,5 +323,5 @@
 #endif
 
-call conf_pem
+call conf_pem(i_myear,n_myear)
 
 !------------------------
@@ -449,5 +454,5 @@
     if (noms(nnq) == "co2") igcm_co2 = nnq
 enddo 
-r = 8.314511E+0*1000.E+0/mugaz
+r = 8.314511*1000./mugaz
 
 !------------------------
@@ -661,9 +666,9 @@
      call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
 #else
-     call iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
+     call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 #endif
 
 if (evol_orbit_pem) then
-    call orbit_param_criterion(year_iter_max)
+    call orbit_param_criterion(i_myear,year_iter_max)
 else
     year_iter_max = Max_iter_pem
@@ -677,6 +682,6 @@
 !------------------------
 year_iter = 0
-do while (year_iter < year_iter_max)
-
+
+do while (year_iter < year_iter_max .and. i_myear < n_myear)
 ! II.a.1. Compute updated global pressure
     write(*,*) "Recomputing the new pressure..."
@@ -911,5 +916,5 @@
 
         call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, &
-                                          delta_h2o_icetablesublim)        ! Mass of H2O exchange between the ssi and the atmosphere.
+                                          delta_h2o_icetablesublim)        ! Mass of H2O exchange between the ssi and the atmosphere
 
         write(*,*) "Update soil propreties"
@@ -962,4 +967,5 @@
 
     year_iter = year_iter + dt_pem
+    i_myear = i_myear + dt_pem
 
     write(*,*) "Checking all the stopping criterion."
@@ -984,4 +990,8 @@
         criterion_stop = 4
     endif
+    if (i_myear >= n_myear) then
+        write(*,*) "STOPPING because maximum number of Martian years to be simulated reached"
+        criterion_stop = 5
+    endif
 
     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
@@ -989,5 +999,6 @@
     else
         write(*,*) "We continue!"
-        write(*,*) "Number of iteration of the PEM : year_iter=", year_iter
+        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
+        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
     endif
 
@@ -1118,5 +1129,5 @@
 enddo
 
-if (evol_orbit_pem) call recomp_orb_param(year_iter)
+if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
 
 !------------------------
@@ -1126,6 +1137,8 @@
 ! III_b.1 Write restart_evol.nc
 ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys
+print*, ptimestep,dtphys,iphysiq*daysec/real(day_step),nsplit_phys
 pday = day_ini
-ztime_fin = 0.
+ptime = 0.
+ztime_fin = pday - day_end + ptime + ptimestep/(float(iphysiq)*daysec)
 
 allocate(p(ip1jmp1,nlayer + 1))
@@ -1133,7 +1146,6 @@
     call pression (ip1jmp1,ap,bp,ps,p)
     call massdair(p,masse)
-    call dynredem0("restart_evol.nc", day_ini, phis)
-    call dynredem1("restart_evol.nc",   &
-    time_0,vcov,ucov,teta,q,masse,ps)
+    call dynredem0("restart_evol.nc",day_ini,phis)
+    call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
     write(*,*) "restart_evol.nc has been written"
 #else
@@ -1173,13 +1185,15 @@
 
 
-call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope ,           &
-             tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
+call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
+             TI_PEM, porefillingice_depth,porefillingice_thickness,         &
              co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
-
-call info_run_PEM(year_iter,criterion_stop)
-
 write(*,*) "restartfi_PEM.nc has been written"
-write(*,*) "The PEM had run for ", year_iter, " years."
-write(*,*) "LL & RV : So far so good"
+
+call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
+
+write(*,*) "The PEM has run for", year_iter, "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!"
 
 deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3039)
@@ -1,31 +1,27 @@
-   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
+SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
                       tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
                       global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
                       m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    
-   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
-   use comsoil_h, only:  volcapa,inertiedat 
-   use adsorption_mod, only : regolith_adsorption,adsorption_pem
-   USE temps_mod_evol, ONLY: year_PEM
-   USE ice_table_mod, only: computeice_table_equilibrium
-   USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
-                                   TI_breccia,TI_bedrock
-   use soil_thermalproperties_mod, only: update_soil_thermalproperties
-   use tracer_mod, only: mmol,igcm_h2o_vap ! tracer names and molar masses
+  use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
+  use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+  use comsoil_h,                  only: volcapa,inertiedat 
+  use adsorption_mod,             only: regolith_adsorption, adsorption_pem
+  use ice_table_mod,              only: computeice_table_equilibrium
+  use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
+  use soil_thermalproperties_mod, only: update_soil_thermalproperties
+  use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
 
 #ifndef CPP_STD
-      USE comcstfi_h, only: r, mugaz
+    use comcstfi_h,   only: r, mugaz
+    use surfdat_h,    only: watercaptag
 #else
-      USE comcstfi_mod, only: r, mugaz
+    use comcstfi_mod, only: r, mugaz
 #endif
 
-
-#ifndef CPP_STD   
-   use surfdat_h, only: watercaptag
-#endif
-
   implicit none
- 
+
+  include "callkeys.h"
+
   character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
   integer,intent(in) :: ngrid                                 ! # of physical grid points
@@ -73,5 +69,4 @@
    real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
    real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
-   real :: year_PEM_read                                              ! Year of the PEM previous run
    LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
 #ifdef CPP_STD   
@@ -97,5 +92,4 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-
 !0. Check if the start_PEM exist.
 
@@ -110,13 +104,5 @@
    call open_startphy(filename)
 
-   call get_var("Time",year_PEM_read,found)
-   year_PEM=INT(year_PEM_read)
-   if(.not.found) then
-      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
-   else
-      write(*,*)'year_PEM of startpem=', year_PEM
-   endif
-
-    if(soil_pem) then
+    if (soil_pem) then
  
 !1. Thermal Inertia
@@ -161,5 +147,5 @@
 ENDDO ! islope
 
-print *,'PEMETAT0: THERMAL INERTIA DONE'
+write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
 
 ! b. Special case for inertiedat, inertiedat_PEM
@@ -204,7 +190,5 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
 !2. Soil Temperature
-
 DO islope=1,nslope
   write(num,fmt='(i2.2)') islope
@@ -260,15 +244,10 @@
         enddo
       enddo
-
-    
 ENDDO ! islope
 
-print *,'PEMETAT0: SOIL TEMP  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !3. Ice Table
-
-
   call get_field("ice_table",ice_table,found)
    if(.not.found) then
@@ -282,10 +261,8 @@
    endif
 
-print *,'PEMETAT0: ICE TABLE  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+write(*,*) 'PEMETAT0: ICE TABLE  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !4. CO2 & H2O Adsorption
-
  if(adsorption_pem) then
   DO islope=1,nslope
@@ -318,13 +295,10 @@
        deltam_h2o_regolith_phys(:) = 0.  
     endif
-print *,'PEMETAT0: CO2 & H2O adsorption done  '
+write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
     endif
 endif ! soil_pem
 
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. 5 water reservoir
-
 #ifndef CPP_STD   
    call get_field("water_reservoir",water_reservoir,found)
@@ -342,12 +316,9 @@
 #endif
 
-
 call close_startphy
 
 else !No startfi, let's build all by hand
 
-    year_PEM=0
-
-    if(soil_pem) then
+    if (soil_pem) then
 
 !a) Thermal inertia
@@ -407,5 +378,4 @@
 enddo
 
-
 !!! zone de transition
 delta = depth_bedrock
@@ -424,9 +394,7 @@
  enddo
 
-print *,'PEMETAT0: TI  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
+write(*,*) 'PEMETAT0: TI DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !b) Soil temperature 
     do islope = 1,nslope
@@ -458,31 +426,24 @@
         enddo
      enddo
- 
-      do isoil = nsoil_GCM+1,nsoil_PEM
-        do ig = 1,ngrid
-        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
+
+        do isoil = nsoil_GCM+1,nsoil_PEM
+          do ig = 1,ngrid
+            watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
+          enddo
         enddo
-      enddo
 enddo !islope
-print *,'PEMETAT0: TSOIL DONE  '
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!c) Ice table    
-     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
-     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
-     do islope = 1,nslope
-     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-     enddo
-    
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-print *,'PEMETAT0: Ice table DONE  '
-
-
+write(*,*) 'PEMETAT0: TSOIL DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!c) Ice table
+       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
+       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
+       do islope = 1,nslope
+         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+       enddo
+       write(*,*) 'PEMETAT0: Ice table DONE'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed
-
  if(adsorption_pem) then
     m_co2_regolith_phys(:,:,:) = 0.
@@ -496,18 +457,10 @@
  endif
 
-print *,'PEMETAT0: CO2 adsorption done  '
-
+write(*,*) 'PEMETAT0: CO2 adsorption DONE'
 endif !soil_pem
 
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. e) water reservoir
-
 #ifndef CPP_STD   
-
       do ig=1,ngrid
         if(watercaptag(ig)) then
@@ -517,5 +470,4 @@
         endif
       enddo
-
 #endif
 
@@ -527,7 +479,5 @@
       DO islope = 1,nslope
         DO iloop = 1,nsoil_PEM
-           if(isnan(tsoil_PEM(ig,iloop,islope))) then 
-         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1) 
-           endif
+          if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1) 
         ENDDO
       ENDDO
@@ -535,5 +485,3 @@
 endif!soil_pem
 
-
-
-   END SUBROUTINE
+END SUBROUTINE
Index: /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3039)
@@ -13,35 +13,32 @@
 contains
 
-subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid, &
-                         day_ini,time,nslope,def_slope,subslope_dist)
+subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid,day_ini,time,nslope,def_slope,subslope_dist)
 ! create physics restart file and write time-independent variables
-  use comsoil_h_PEM, only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
-  use iostart_PEM, only : open_restartphy, close_restartphy, & 
-                      put_var, put_field, length
-  use mod_grid_phy_lmdz, only : klon_glo
+  use comsoil_h_PEM,     only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
+  use iostart_PEM,       only: open_restartphy, close_restartphy, put_var, put_field, length
+  use mod_grid_phy_lmdz, only: klon_glo
+  use time_phylmdz_mod,  only: daysec
 #ifndef CPP_STD
-  use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
-                       peri_day, periheli, year_day
-  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
+  use planete_h,    only: aphelie, emin_turb, lmixmin, obliquit, peri_day, periheli, year_day
+  use comcstfi_h,   only: g, mugaz, omeg, rad, rcp
 #else
-  use planete_mod, only: apoastr, emin_turb, lmixmin, obliquit, &
-                       peri_day, periastr, year_day
-  USE comconst_mod, ONLY: g, omeg, rad
+  use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
+  use comconst_mod, only: g, omeg, rad
   use comcstfi_mod, only: mugaz, rcp
 #endif
-  use time_phylmdz_mod, only: daysec
+
   implicit none
  
   character(len=*), intent(in) :: filename
-  real,intent(in) :: lonfi(ngrid)
-  real,intent(in) :: latfi(ngrid)
-  integer,intent(in) :: nsoil_PEM
-  integer,intent(in) :: ngrid
-  real,intent(in) :: day_ini
-  real,intent(in) :: time
-  real, intent(in) :: cell_area(ngrid) !boundaries for bining of the slopes
-  integer,intent(in) :: nslope
-  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
-  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
+  real, intent(in)             :: lonfi(ngrid)
+  real, intent(in)             :: latfi(ngrid)
+  integer, intent(in)          :: nsoil_PEM
+  integer, intent(in)          :: ngrid
+  real, intent(in)             :: day_ini
+  real, intent(in)             :: time
+  real, intent(in)             :: cell_area(ngrid) !boundaries for bining of the slopes
+  integer, intent(in)          :: nslope
+  real, intent(in)             :: def_slope(nslope+1) !boundaries for bining of the slopes
+  real, intent(in)             :: subslope_dist(ngrid,nslope) !undermesh statistics
   
   ! Create physics start file
@@ -70,5 +67,5 @@
 end subroutine pemdem0
 
-subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
+subroutine pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope, &
                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table_depth,ice_table_thickness, &
                     m_co2_regolith,m_h2o_regolith,water_reservoir)
@@ -76,9 +73,7 @@
 
   ! write time-dependent variable to restart file
-  use iostart_PEM, only : open_restartphy, close_restartphy, & 
-                      put_var, put_field
-
-  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM,soil_pem
-  USE temps_mod_evol, ONLY: year_PEM
+  use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
+  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo, inertiedat_PEM,soil_pem
+  use time_evol_mod, only: year_bp_ini, convert_years
                 
   implicit none
@@ -88,27 +83,19 @@
 #endif
   
-  character(len=*),intent(in) :: filename
-  integer,intent(in) :: nsoil_PEM
-  integer,intent(in) :: ngrid
-  real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
-  real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
-  real,intent(IN) :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
-  real,intent(IN) :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
-  integer,intent(IN) :: nslope
-  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
-  real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
-  real,intent(in) :: water_reservoir(ngrid)
-  integer :: islope
-  CHARACTER*2 :: num  
-  integer, intent(IN) :: year_iter
+  character(len=*), intent(in) :: filename
+  integer, intent(in)          :: nsoil_PEM, ngrid, nslope, i_myear
+  real, intent(in)             :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
+  real, intent(in)             :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
+  real, intent(in)             :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
+  real, intent(in)             :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
+  real, intent(in)             :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
+  real, intent(in)             :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
+  real, intent(in)             :: water_reservoir(ngrid)
 
-  real :: year_tot
- 
-  integer :: iq
-  character(len=30) :: txt ! to store some text
-  ! indexes of water vapour & water ice tracers (if any):
-
-  year_tot=real(year_iter+year_PEM)
-  print *, "Year of simulation=", year_tot  
+  integer           :: islope
+  character*2       :: num
+  integer           :: iq
+  character(len=30) :: txt  ! To store some text
+  real              :: Year ! Year of the simulation
 
   ! Open file
@@ -117,32 +104,29 @@
   ! First variable to write must be "Time", in order to correctly
   ! set time counter in file
-  call put_var("Time","Year of simulation",year_tot)
-  call put_field('water_reservoir','water_reservoir', &
-           water_reservoir,year_tot)
+  Year = (year_bp_ini + i_myear)*convert_years
+  call put_var("Time","Year of simulation",Year)
+  call put_field('water_reservoir','water_reservoir',water_reservoir,Year)
+
     if(soil_pem) then
   
   ! Multidimensionnal variables (undermesh slope statistics)
-DO islope=1,nslope
+do islope = 1,nslope
   write(num,fmt='(i2.2)') islope
   call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
-                 tsoil_slope_PEM(:,:,islope),year_tot)
+                 tsoil_slope_PEM(:,:,islope),Year)
   call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
-                 inertiesoil_slope_PEM(:,:,islope),year_tot)
+                 inertiesoil_slope_PEM(:,:,islope),Year)
 
- call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
-                    m_co2_regolith(:,:,islope), year_tot)
+  call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
+                m_co2_regolith(:,:,islope),Year)
   call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith", &
-                    m_h2o_regolith(:,:,islope), year_tot)
+                 m_h2o_regolith(:,:,islope),Year)
+enddo
 
-ENDDO
+  call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
 
-  call put_field("ice_table_depth","Depth of ice table", &
-                 ice_table_depth,year_tot)
+  call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
 
-  call put_field("ice_table_thickness","Depth of ice table", &
-                 ice_table_thickness,year_tot)
-
-  call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
-                 inertiedat_PEM,year_tot)
+  call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
 
    endif !soil_pem
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3038)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3039)
@@ -1,3 +1,3 @@
-      MODULE recomp_orb_param_mod
+MODULE recomp_orb_param_mod
 !=======================================================================
 !
@@ -6,30 +6,28 @@
 ! Author: RV, JBC
 !=======================================================================
-      IMPLICIT NONE
+IMPLICIT NONE
 
-      CONTAINS
+CONTAINS
 
-      SUBROUTINE recomp_orb_param(final_iter)
+SUBROUTINE recomp_orb_param(i_myear,year_iter)
 
-      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
+  use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
 #ifndef CPP_STD      
-      USE comconst_mod, ONLY: pi
-      USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day
+    use comconst_mod, only: pi
+    use planete_h,    only: e_elips, obliquit, timeperi, periheli, aphelie, p_elips, peri_day, year_day
 #else
-      use planete_mod, only: e_elips, obliquit, timeperi,periastr,apoastr,p_elips,peri_day,year_day
-      USE comcstfi_mod, only: pi
-
+    use planete_mod,  only: e_elips, obliquit, timeperi, periastr, apoastr, p_elips, peri_day, year_day
+    use comcstfi_mod, only: pi
 #endif
 
-      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
-   				end_lask_param_mod,last_ilask
+  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
 
-      IMPLICIT NONE
+  IMPLICIT NONE
 
 !--------------------------------------------------------
 ! Input Variables
 !--------------------------------------------------------
-
-      integer,intent(in) :: final_iter ! Number of year iteration of the PEM 
+  integer, intent(in) :: i_myear   ! Number of simulated Martian years
+  integer, intent(in) :: year_iter ! Number of iterations of the PEM
 
 !--------------------------------------------------------
@@ -40,89 +38,85 @@
 ! Local variables 
 !--------------------------------------------------------
+  integer :: Year           ! Year of the simulation
+  real    :: Lsp            ! Ls perihelion
+  integer :: ilask          ! Loop variable
+  real    :: halfaxe        ! Million of km
+  real    :: unitastr
+  real    :: xa, xb, ya, yb ! Variables for interpolation
 
-      real :: Year     ! Year of the simulation
-      real :: Lsp      ! time of peri in ls
-      integer ilask    ! Loop variable
-      real :: halfaxe  ! Million of km
-      real :: unitastr
+! **********************************************************************
+! 0. Initializations
+! **********************************************************************
+#ifdef CPP_STD 
+    real aphelie
+    real periheli
 
-      real xa,xb,ya,yb
-
-
-      ! **********************************************************************
-      ! 0. Initializations
-      ! **********************************************************************
-#ifdef CPP_STD 
-      real aphelie
-      real periheli
-
-      aphelie=apoastr
-      periheli=periastr
+    aphelie=apoastr
+    periheli=periastr
 #endif
 
-          Year = year_bp_ini + year_PEM + final_iter ! We compute the current year
-          Lsp = 360. - timeperi*360./(2.*pi)         ! We convert in degree
+Year = year_bp_ini + i_myear        ! We compute the current year
+Lsp = modulo(-timeperi*180./pi,360.) ! We convert in degree
 
-!          Year=-953200
-          print*, "recomp_orb_param, Old year=", year_bp_ini+year_PEM
-          print*, "recomp_orb_param, New year=", Year
-          print*, "recomp_orb_param, Old obl=", obliquit
-          print*, "recomp_orb_param, Old ex=", e_elips
-          print*, "recomp_orb_param, Old lsp=", Lsp
-          print*, "recomp_orb_param, Old timeperi=", timeperi
+write(*,*) 'recomp_orb_param, Old year =', Year - year_iter
+write(*,*) 'recomp_orb_param, Old obl. =', obliquit
+write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
+write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp
 
-        do ilask = last_ilask,2,-1
-          xa = yearlask(ilask)
-          xb = yearlask(ilask - 1)
-          if (xa <= Year .and. Year < xb) then
-            ! Obliquity
-            if (var_obl) then
-              ya = oblask(ilask)
-              yb = oblask(ilask - 1)
-              obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
+do ilask = last_ilask,2,-1
+    xa = yearlask(ilask)
+    xb = yearlask(ilask - 1)
+    if (xa <= Year .and. Year < xb) then
+        ! Obliquity
+        if (var_obl) then
+            ya = obllask(ilask)
+            yb = obllask(ilask - 1)
+            obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya
+        endif
+
+        ! Eccentricity
+        if (var_ecc) then
+            ya = ecclask(ilask)
+            yb = ecclask(ilask - 1)
+            e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
+        endif
+
+        ! Lsp
+        if (var_lsp) then
+            ya = lsplask(ilask)
+            yb = lsplask(ilask - 1)
+            if (yb - ya > 300.) then ! If modulo is "crossed"
+                if (ya < yb) then ! Increasing
+                    yb = yb - 360.
+                else ! Decreasing
+                    yb = yb + 360.
+                endif
             endif
+            Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya
+            Lsp = modulo(Lsp,360.)
+        endif
+        exit ! The loop is left as soon as the right interval is found
+    endif
+enddo
+halfaxe = 227.94
+timeperi = Lsp*pi/180.
+periheli = halfaxe*(1 - e_elips)
+aphelie =  halfaxe*(1 + e_elips)
+unitastr = 149.597927
+p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
 
-            ! Excentricity
-            if (var_ex) then
-              ya = exlask(ilask)
-              yb = exlask(ilask - 1)
-              e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya
-            endif
-
-            ! Lsp
-            if (var_lsp) then
-              ya = lsplask(ilask)
-              yb = lsplask(ilask - 1)
-              if (yb - ya > 300.) then ! If modulo is "crossed"
-                if (ya < yb) then ! Increasing
-                  yb = yb - 360.
-                else ! Decreasing
-                  yb = yb + 360.
-                endif
-              endif
-              Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya
-              Lsp = modulo(Lsp,360.)
-            endif
-            exit ! The loop is left as soon as the right interval is found
-          endif
-        enddo
-        halfaxe = 227.94
-        timeperi = Lsp*2.*pi/360.
-        periheli = halfaxe*(1 - e_elips)
-        aphelie =  halfaxe*(1 + e_elips)
-        unitastr = 149.597927
-        p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
 #ifndef CPP_STD  
-        call call_dayperi(timeperi,e_elips,peri_day,year_day)
+    call call_dayperi(timeperi,e_elips,peri_day,year_day)
 #endif
 
-       print*, "recomp_orb_param, Final year of the PEM run:", Year
-       print*, "recomp_orb_param, New obl=", obliquit
-       print*, "recomp_orb_param, New ex=", e_elips
-       print*, "recomp_orb_param, New timeperi=", timeperi
+write(*,*) 'recomp_orb_param, New year =', Year
+write(*,*) 'recomp_orb_param, New obl. =', obliquit
+write(*,*) 'recomp_orb_param, New ecc. =', e_elips
+write(*,*) 'recomp_orb_param, New Lsp  =', Lsp
 
-      END SUBROUTINE recomp_orb_param
+END SUBROUTINE recomp_orb_param
 
 !********************************************************************************   
       
-      END MODULE recomp_orb_param_mod
+END MODULE recomp_orb_param_mod
+
Index: unk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 3038)
+++ 	(revision )
@@ -1,17 +1,0 @@
-MODULE temps_mod_evol
-
-IMPLICIT NONE  
-
-  INTEGER   year_bp_ini         !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
-  INTEGER   dt_pem              !     dt_pem  : in Planetary years, the time step used by the PEM    
-  REAL      water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM 
-  REAL      co2_ice_criterion   !     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   year_PEM            !     year written in startfiPEM.nc
-  INTEGER   Max_iter_pem        !     Maximal number of iteration when converging to a steady state, read in evol.def
-  LOGICAL   evol_orbit_pem      !     True if we want to follow the orbital parameters of ob_ex_lsp.asc, read in evol.def
-  LOGICAL   var_obl             !     True if we want the PEM to follow ob_ex_lsp.asc parameters for obliquity 
-  LOGICAL   var_ex              !     True if we want the PEM to follow ob_ex_lsp.asc parameters for excenticity
-  LOGICAL   var_lsp             !     True if we want the PEM to follow ob_ex_lsp.asc parameters for ls perihelie
-
-END MODULE temps_mod_evol
Index: /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3039)
+++ /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3039)
@@ -0,0 +1,17 @@
+MODULE time_evol_mod
+
+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    :: convert_years       !     Conversion ratio from Planetary years to Earth years
+  real    :: water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM 
+  real    :: co2_ice_criterion   !     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
+  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 
+  logical :: var_ecc             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
+  logical :: var_lsp             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
+
+END MODULE time_evol_mod
Index: /trunk/LMDZ.MARS/changelog.txt
===================================================================
--- /trunk/LMDZ.MARS/changelog.txt	(revision 3038)
+++ /trunk/LMDZ.MARS/changelog.txt	(revision 3039)
@@ -4185,2 +4185,10 @@
 * Some changes in run_PEM.def and xml files.
 JBC
+
+
+== 11/09/2023 == JBC
+* New management of time in PEM. While definitions in "run.def"/"launching script" and Laskar's data are in Earth years, the PEM works in Martian years. It follows several changes and the new variable 'convert_years' to make the conversion;
+* New parameter for the years to be simulated. The user does not ask anymore for a number of PEM iterations to do but for a number of Earth years to simulate. The GCM years are now counted! To do so, a temporary file "tmp_PEMyears.txt" gives few basic information between the runs;
+* New values for the maximal admissible change of orbital parameters. They are now coherent and allows the PEM to run for ~1000 Martian years;
+* Laskar's data interpolation has been cleaned further;
+* Some cleaning and renaming of PEM subroutines in the course of the previous modifications.
