Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2991)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2994)
@@ -54,4 +54,5 @@
       use aeropacity_mod, only: iddist, topdustref
       USE mod_phys_lmdz_transfert_para, ONLY: bcast
+      USE paleoclimate_mod,ONLY: paleoclimate
       IMPLICIT NONE
       include "callkeys.h"
@@ -346,4 +347,11 @@
          call getin_p("rayleigh",rayleigh)
          write(*,*)" rayleigh = ",rayleigh
+
+! PALEOCLIMATE
+
+         write(*,*)"Is it paleoclimate run?"
+         paleoclimate=.false. ! default value
+         call getin_p("paleoclimate",paleoclimate)
+         write(*,*)" paleoclimate = ",paleoclimate
 
 
Index: trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 2994)
+++ trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 2994)
@@ -0,0 +1,38 @@
+MODULE paleoclimate_mod
+!=======================================================================
+!   subject: Module dedicated to paleoclimates studies
+!   --------
+!
+!   author: LL, 06/2023
+!   ------
+!   
+!=======================================================================
+
+    IMPLICIT NONE
+
+    LOGICAL :: paleoclimate                        ! False by default, is activate for paleoclimates specific processes (e.g., lag layer)
+    real, save, allocatable :: lag_h2o_ice(:,:)  ! Thickness of the lag before H2O ice [m]
+    real, save, allocatable :: lag_co2_ice(:,:)  ! Thickness of the lag before CO2 ice [m]
+
+    CONTAINS
+
+
+  subroutine ini_paleoclimate_h(ngrid,nslope)
+
+  implicit none
+  integer,intent(in) :: ngrid  ! number of atmospheric columns
+  integer,intent(in) :: nslope ! number of slope within a mesh 
+
+    allocate(lag_h2o_ice(ngrid,nslope))
+    allocate(lag_co2_ice(ngrid,nslope))
+  end subroutine ini_paleoclimate_h
+
+  subroutine end_paleoclimate_h
+
+  implicit none
+    if (allocated(lag_h2o_ice)) deallocate(lag_h2o_ice)
+    if (allocated(lag_co2_ice)) deallocate(lag_co2_ice)
+  end subroutine end_paleoclimate_h
+
+
+    END MODULE paleoclimate_mod
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2991)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2994)
@@ -27,4 +27,5 @@
   use comsoil_h, only: flux_geo
   USE comslope_mod, ONLY: nslope, major_slope
+  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice
   implicit none
   
@@ -752,5 +753,33 @@
 endif !startphy_file
 
-!
+
+if (paleoclimate) then
+  if (startphy_file) then
+   call get_field("lag_h2o_ice",lag_h2o_ice,found,indextime)
+   write(*,*) 'paleo found start?',found
+   
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <lag_h2o_ice> : ", &
+                          "<lag_h2o_ice> is set as -1 (no subsurface ice)"
+     lag_h2o_ice(:,:) = -1.
+   endif
+
+   call get_field("lag_co2_ice",lag_co2_ice,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <lag_co2_ice> : ", &
+                          "<lag_co2_ice> is set as -1 (no subsurface ice)"
+     lag_co2_ice(:,:) = -1.
+   endif
+  else
+     lag_h2o_ice(:,:) = -1.
+     lag_co2_ice(:,:) = -1.
+  endif !startphy_file
+else
+   write(*,*) 'paleo found? nostart',found
+   
+  lag_h2o_ice(:,:) = -1.
+  lag_co2_ice(:,:) = -1.
+endif !paleoclimate
+
 ! close file:
 !
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2991)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2994)
@@ -25,4 +25,5 @@
   use time_phylmdz_mod, only: daysec
   use comslope_mod, ONLY: nslope
+  USE paleoclimate_mod, ONLY: paleoclimate, lag_h2o_ice, lag_co2_ice  
   implicit none
  
@@ -46,5 +47,5 @@
   integer :: ig
   real :: watercaptag_tmp(ngrid)
-  
+
   ! Create physics start file
   call open_restartphy(filename)
@@ -153,4 +154,9 @@
   call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
 
+  ! Paleoclimate outputs
+  if(paleoclimate) then
+    call put_field("lag_h2o_ice","Depth of the H2O lags",lag_h2o_ice)
+    call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
+  endif
   ! Close file
   call close_restartphy
Index: trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2991)
+++ trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2994)
@@ -62,4 +62,5 @@
                                      end_dust_rad_adjust_mod
       use comslope_mod, ONLY: nslope,end_comslope_h,ini_comslope_h
+      use paleoclimate_mod, ONLY: end_paleoclimate_h,ini_paleoclimate_h
       use netcdf
       USE mod_phys_lmdz_para, ONLY: is_master,bcast
@@ -186,4 +187,8 @@
       call ini_comslope_h(ngrid,nslope)
 
+      ! allocate arrays in "paleoclimate_mod"
+      call end_paleoclimate_h
+      call ini_paleoclimate_h(ngrid,nslope)
+
       END SUBROUTINE phys_state_var_init