Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3149)
@@ -154,2 +154,16 @@
 - Addition of a script in LMDZ.MARS/deftank/pem/ to launch a chained simulation of 1D PCM runs which follow, year by year, the orbital parameters (obliquity, eccentricity, Ls perihelion) given in a specified file.
 - Small changes to other files of the deftank directory (check and cosmetic).
+
+== 06/12/2023 == JBC
+- Simplification of the algorithm managing the stopping criteria;
+- Complete rework of the ice management in the PEM (H2O & CO2);
+      > Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
+      > Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
+      > Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
+      > Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
+      > Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
+      > Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
+      > New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').
+- Renaming of variables/subroutines for clarity;
+- Some cleanings throughout the code;
+- Small updates in files of the deftank.
Index: trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90	(revision 3149)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_tend_mod.F90	(revision 3149)
@@ -0,0 +1,37 @@
+MODULE compute_tend_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE compute_tend(ngrid,nslope,min_ice,tendencies_ice)
+
+implicit none
+
+!=======================================================================
+!
+! Compute the tendencies of the evolution of water ice over the years
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,                         intent(in) :: ngrid   ! # of grid points
+integer,                         intent(in) :: nslope  ! # of subslopes
+real, dimension(ngrid,nslope,2), intent(in) :: min_ice ! Minima of ice at each point for the PCM years
+
+!   OUTPUT
+real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! Difference between the minima = evolution of perennial ice
+!=======================================================================
+! We compute the difference
+tendencies_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.
+
+END SUBROUTINE compute_tend
+
+END MODULE compute_tend_mod
Index: trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90	(revision 3145)
+++ 	(revision )
@@ -1,39 +1,0 @@
-MODULE compute_tendencies_slope_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-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, dimension(ngrid,nslope), intent(in) :: min_ice_Y1    ! LON x LAT field : minimum of water ice at each point for the first year
-real, dimension(ngrid,nslope), intent(in) :: min_ice_Y2    ! LON x LAT field : minimum of water ice at each point for the second year
-
-!   OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perennial ice
-!=======================================================================
-
-!  We compute the difference
-tendencies_ice = min_ice_Y2 - min_ice_Y1
-
-!  If the difference is too small; there is no evolution
-where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
-
-END SUBROUTINE compute_tendencies_slope
-
-END MODULE compute_tendencies_slope_mod
-
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 3149)
@@ -3,27 +3,26 @@
 implicit none
 
-integer, parameter                        :: nsoilmx_PEM = 68    ! number of layers in the PEM
-real, allocatable, dimension(:),     save :: layer_PEM           ! soil layer depths [m]
-real, allocatable, dimension(:),     save :: mlayer_PEM          ! soil mid-layer depths [m]
-real, allocatable, dimension(:,:,:), save :: TI_PEM              ! soil thermal inertia [SI]
-real, allocatable, dimension(:,:),   save :: inertiedat_PEM      ! soil thermal inertia saved as reference for current climate [SI]
+integer, parameter                        :: nsoilmx_PEM = 68     ! number of layers in the PEM
+real, allocatable, dimension(:),     save :: layer_PEM            ! soil layer depths [m]
+real, allocatable, dimension(:),     save :: mlayer_PEM           ! soil mid-layer depths [m]
+real, allocatable, dimension(:,:,:), save :: TI_PEM               ! soil thermal inertia [SI]
+real, allocatable, dimension(:,:),   save :: inertiedat_PEM       ! soil thermal inertia saved as reference for current climate [SI]
 ! Variables (FC: built in firstcall in soil.F)
-real, allocatable, dimension(:,:,:), save :: tsoil_PEM           ! sub-surface temperatures [K]
-real, allocatable, dimension(:,:),   save :: mthermdiff_PEM      ! (FC) mid-layer thermal diffusivity [SI]
-real, allocatable, dimension(:,:),   save :: thermdiff_PEM       ! (FC) inter-layer thermal diffusivity [SI]
-real, allocatable, dimension(:),     save :: coefq_PEM           ! (FC) q_{k+1/2} coefficients [SI]
-real, allocatable, dimension(:,:),   save :: coefd_PEM           ! (FC) d_k coefficients [SI]
-real, allocatable, dimension(:,:),   save :: alph_PEM            ! (FC) alpha_k coefficients [SI]
-real, allocatable, dimension(:,:),   save :: beta_PEM            ! beta_k coefficients [SI]
-real,                                save :: mu_PEM              ! mu coefficient [SI]
-real,                                save :: fluxgeo             ! Geothermal flux [W/m^2]
-real,                                save :: depth_breccia       ! Depth at which we have breccia [m]
-real,                                save :: depth_bedrock       ! Depth at which we have bedrock [m]
-integer,                             save :: index_breccia       ! last index of the depth grid before having breccia
-integer,                             save :: index_bedrock       ! last index of the depth grid before having bedrock
-logical                                   :: soil_pem            ! True by default, to run with the subsurface physic. Read in pem.def
-real, allocatable, dimension(:),     save :: water_reservoir     ! Large reserve of water [kg/m^2]
-real,                                save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2]
-logical,                             save :: reg_thprop_dependp  ! thermal properites of the regolith vary with the surface pressure
+real, allocatable, dimension(:,:,:), save :: tsoil_PEM            ! sub-surface temperatures [K]
+real, allocatable, dimension(:,:),   save :: mthermdiff_PEM       ! (FC) mid-layer thermal diffusivity [SI]
+real, allocatable, dimension(:,:),   save :: thermdiff_PEM        ! (FC) inter-layer thermal diffusivity [SI]
+real, allocatable, dimension(:),     save :: coefq_PEM            ! (FC) q_{k+1/2} coefficients [SI]
+real, allocatable, dimension(:,:),   save :: coefd_PEM            ! (FC) d_k coefficients [SI]
+real, allocatable, dimension(:,:),   save :: alph_PEM             ! (FC) alpha_k coefficients [SI]
+real, allocatable, dimension(:,:),   save :: beta_PEM             ! beta_k coefficients [SI]
+real,                                save :: mu_PEM               ! mu coefficient [SI]
+real,                                save :: fluxgeo              ! Geothermal flux [W/m^2]
+real,                                save :: depth_breccia        ! Depth at which we have breccia [m]
+real,                                save :: depth_bedrock        ! Depth at which we have bedrock [m]
+integer,                             save :: index_breccia        ! last index of the depth grid before having breccia
+integer,                             save :: index_bedrock        ! last index of the depth grid before having bedrock
+logical                                   :: soil_pem             ! True by default, to run with the subsurface physic. Read in pem.def
+real,                                save :: ini_h2o_bigreservoir ! Initial value for the large reservoir of water [kg/m^2]
+logical,                             save :: reg_thprop_dependp   ! thermal properites of the regolith vary with the surface pressure
 
 !=======================================================================
@@ -39,15 +38,14 @@
 
 allocate(layer_PEM(nsoilmx_PEM))
-allocate(mlayer_PEM(0:nsoilmx_PEM-1))
+allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
 allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
 allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
-allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
-allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
-allocate(coefq_PEM(0:nsoilmx_PEM-1))
-allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
-allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
-allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
+allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1))
+allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1))
+allocate(coefq_PEM(0:nsoilmx_PEM - 1))
+allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1))
+allocate(alph_PEM(ngrid,nsoilmx_PEM - 1))
+allocate(beta_PEM(ngrid,nsoilmx_PEM - 1))
 allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
-allocate(water_reservoir(ngrid))
 
 END SUBROUTINE ini_comsoil_h_PEM
@@ -70,5 +68,4 @@
 if (allocated(beta_PEM)) deallocate(beta_PEM)
 if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
-if (allocated(water_reservoir)) deallocate(water_reservoir)
 
 END SUBROUTINE end_comsoil_h_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3149)
@@ -25,7 +25,7 @@
 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 comsoil_h_pem,    only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
 use adsorption_mod,   only: adsorption_pem
-use glaciers_mod,     only: co2glaciersflow, h2oglaciersflow
+use glaciers_mod,     only: co2_ice_flow, h2o_ice_flow
 use ice_table_mod,    only: icetable_equilibrium, icetable_dynamic
 use time_phylmdz_mod, only: ecritphy
@@ -71,5 +71,5 @@
 var_ecc = .true.
 call getin('var_ecc',var_ecc)
-write(*,*) 'Does excentricity vary ?',var_ecc
+write(*,*) 'Does eccentricity vary ?',var_ecc
 
 var_lsp = .true.
@@ -100,9 +100,9 @@
 call getin('adsorption_pem',adsorption_pem)
 
-co2glaciersflow = .true.
-call getin('co2glaciersflow',co2glaciersflow)
+co2_ice_flow = .true.
+call getin('co2_ice_flow',co2_ice_flow)
 
-h2oglaciersflow = .true.
-call getin('h2oglaciersflow',h2oglaciersflow)
+h2o_ice_flow = .true.
+call getin('h2o_ice_flow',h2o_ice_flow)
 
 reg_thprop_dependp = .false.
@@ -156,6 +156,6 @@
 endif
 
-water_reservoir_nom = 1.e4
-call getin('water_reservoir_nom',water_reservoir_nom)
+ini_h2o_bigreservoir = 1.e4
+call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
 
 END SUBROUTINE conf_pem
Index: trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3149)
@@ -12,32 +12,31 @@
 real, parameter :: m_h2o = 18.01528e-3 ! Molecular weight of h2o (kg/mol)
 
-!     Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
+! Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
 real, parameter :: alpha_clap_co2 = 23.3494 ! Uniteless, James et al. 1992
 real, parameter :: beta_clap_co2 = 3182.48  ! Kelvin, James et al. 1992
 
-!     Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
+! Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
 real, parameter :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphy and Koop 2005
 real, parameter :: beta_clap_h2o = -6143.7  ! Kelvin, Murphy and Koop 2005
 
-!     Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
+! Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
 real, parameter :: rho_regolith = 2000. ! kg/m^3
 
-!     Average  Thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
+! Average thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
 real, parameter :: TI_regolith_avg = 250. ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
 real, parameter :: TI_breccia = 750.      ! Thermal inertia of Breccia following Wood 2009 [SI]
 real, parameter :: TI_bedrock = 2300.     ! Thermal inertia of Bedrock following Wood 2009 [SI]
 
-!     Porosity of the soil
+! Porosity of the soil
 real, parameter :: porosity = 0.4 ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
 
-!     Stefan Boltzmann constant
+! Stefan Boltzmann constant
 real, parameter :: sigmaB = 5.678e-8 
 
-!     Latent heat of CO2
+! Latent heat of CO2
 real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
 
-! Conversion H2O/CO2 frost to perennial frost and vice versa
-real, parameter :: threshold_h2o_frost2perennial = 1000.  !~ 1 m
-real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m
+! Threshold to consider the amount of H2O ice as an infinite reservoir
+real, parameter :: inf_h2oice_threshold = 2000.  !~ 1 m
 
 END MODULE constants_marspem_mod
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 3145)
+++ 	(revision )
@@ -1,163 +1,0 @@
-MODULE criterion_pem_stop_mod
-
-implicit none
-
-! ----------------------------------------------------------------------
-contains
-! ----------------------------------------------------------------------
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!
-!!! Purpose: Criterions to check if the PEM needs to call the PCM
-!!! Author: RV & LL, 02/2023
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
-
-use time_evol_mod, only: water_ice_criterion
-use comslope_mod,  only: subslope_dist, nslope
-
-implicit none
-
-!=======================================================================
-!
-!  Routine that checks if the water ice criterion to stop the PEM is reached
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-!   INPUT
-integer,                       intent(in) :: ngrid           ! # of physical grid points
-real, dimension(ngrid),        intent(in) :: cell_area       ! Area of the cells
-real, dimension(ngrid,nslope), intent(in) :: qsurf           ! Actual density of water ice
-real,                          intent(in) :: ini_surf        ! Initial surface of h2o ice that was sublimating
-real, dimension(ngrid,nslope), intent(in) :: initial_h2o_ice ! Grid point that initialy were covered by h2o_ice
-!   OUTPUT
-logical, intent(out) :: STOPPING ! Is the criterion reached?
-!   local:
-!   ------
-integer :: i, islope    ! Loop
-real    :: present_surf ! Initial/Actual surface of water ice
-
-!=======================================================================
-! Initialisation to false
-STOPPING = .false.
-
-! Computation of the present surface of h2o ice
-present_surf = 0.
-do i = 1,ngrid
-    do islope = 1,nslope
-        if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
-    enddo
-enddo
-
-! Check of the criterion
-if (present_surf < ini_surf*(1. - water_ice_criterion)) then
-    STOPPING = .true.
-    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
-    write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion)
-    write(*,*) "Current surface of water ice sublimating =", present_surf
-    write(*,*) "Initial surface of water ice sublimating =", ini_surf
-    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
-endif
-if (present_surf > ini_surf*(1. + water_ice_criterion)) then
-    STOPPING = .true.
-    write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
-    write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion)
-    write(*,*) "Current surface of water ice sublimating =", present_surf
-    write(*,*) "Initial surface of water ice sublimating =", ini_surf
-    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
-endif
-
-if (ini_surf < 1.e-5 .and. ini_surf > -1.e-5) STOPPING = .false.
-
-END SUBROUTINE criterion_waterice_stop
-
-! ----------------------------------------------------------------------
-
-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 time_evol_mod, only: co2_ice_criterion, ps_criterion
-use comslope_mod,  only: subslope_dist
-
-implicit none
-
-!=======================================================================
-!
-!  Routine that checks if the co2 and pressure criteria to stop the PEM are reached
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-!   INPUT
-integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
-real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
-real, dimension(ngrid,nslope), intent(in) :: qsurf                ! Actual density of water ice
-real,                          intent(in) :: ini_surf             ! Initial surface of co2 ice that was sublimating
-real, dimension(ngrid,nslope), intent(in) :: initial_co2_ice      ! Grid point that initialy were covered by co2_ice
-real,                          intent(in) :: global_ave_press_GCM ! Planet average pressure from the PCM start files
-real,                          intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations
-!   OUTPUT
-logical, intent(out) :: STOPPING_ice ! Is the criterion for co2 ice reached?
-logical, intent(out) :: STOPPING_ps  ! Is the criterion for pressure reached?
-!   local:
-!   ------
-integer :: i, islope    ! Loop
-real    :: present_surf ! Initial/Actual surface of water ice
-
-!=======================================================================
-! Initialisation to false
-STOPPING_ice = .false.
-STOPPING_ps = .false.
-
-! Computation of the present surface of co2 ice
-present_surf = 0.
-do i = 1,ngrid
-    do islope = 1,nslope
-        if (initial_co2_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
-    enddo
-enddo
-
-! Check of the criterion
-if (present_surf < ini_surf*(1. - co2_ice_criterion)) then
-    STOPPING_ice = .true.
-    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold"
-    write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion)
-    write(*,*) "Current surface of co2 ice sublimating =", present_surf
-    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
-    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
-endif
-if (present_surf > ini_surf*(1. + co2_ice_criterion)) then
-    STOPPING_ice = .true.
-    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
-    write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion)
-    write(*,*) "Current surface of co2 ice sublimating =", present_surf
-    write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
-    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
-endif
-
-if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.
-
-if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then
-    STOPPING_ps = .true.
-    write(*,*) "Reason of stopping: the global pressure reach the threshold"
-    write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)
-    write(*,*) "Current global pressure =", global_ave_press_new
-    write(*,*) "PCM global pressure =", global_ave_press_GCM
-    write(*,*) "Percentage of change accepted =", ps_criterion*100.
-endif
-if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
-    STOPPING_ps = .true.
-    write(*,*) "Reason of stopping: the global pressure reach the threshold"
-    write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)
-    write(*,*) "Current global pressure =", global_ave_press_new
-    write(*,*) "PCM global pressure =", global_ave_press_GCM
-    write(*,*) "Percentage of change accepted =", ps_criterion*100.
-endif
-
-END SUBROUTINE criterion_co2_stop
-
-END MODULE
Index: trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod.F90	(revision 3145)
+++ 	(revision )
@@ -1,47 +1,0 @@
-MODULE evol_co2_ice_s_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE evol_co2_ice_s(qsurf,tendencies_co2_ice_phys,iim_input,jjm_input,ngrid,cell_area,nslope)
-
-use time_evol_mod, only: dt_pem
-
-implicit none
-
-!=======================================================================
-!
-!  Routine that compute the evolution of the CO2 ice
-!
-!=======================================================================
-!   arguments:
-!   ----------
-!   INPUT
-integer,                intent(in) :: iim_input, jjm_input, ngrid, nslope ! # of grid points along longitude/latitude/ total
-REAL, dimension(ngrid), intent(in) :: cell_area
-
-!   OUTPUT
-real, dimension(ngrid,nslope), intent(inout) :: qsurf                   ! physical point field: Previous and actual density of water ice
-real, dimension(ngrid,nslope), intent(inout) :: tendencies_co2_ice_phys ! physical point field: Evolution of perennial ice over one year
-
-!   local:
-!   ----
-integer :: i, j, ig0, islope ! loop variable
-!=======================================================================
-! Evolution of the CO2 ice for each physical point
-do i=1,ngrid
-    do islope=1,nslope
-        qsurf(i,islope) = qsurf(i,islope) + tendencies_co2_ice_phys(i,islope)*dt_pem
-        if (qsurf(i,islope) < 0) then
-            qsurf(i,islope) = 0.
-            tendencies_co2_ice_phys(i,islope) = 0.
-        endif
-    enddo
-enddo
-
-END SUBROUTINE evol_co2_ice_s
-
-END MODULE evol_co2_ice_s_mod
Index: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod.F90	(revision 3145)
+++ 	(revision )
@@ -1,147 +1,0 @@
-MODULE evol_h2o_ice_s_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf,tendencies_h2o_ice_phys,STOPPING)
-
-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
-
-#ifndef CPP_STD
-    use comcstfi_h,   only: pi
-#else
-    use comcstfi_mod, only: pi
-#endif
-
-implicit none
-
-!=======================================================================
-!
-!  Routine that compute the evolution of the water ice
-!
-!=======================================================================
-!   arguments:
-!   ----------
-!   INPUT
-integer,                intent(in) :: ngrid                    ! # of grid points along longitude/latitude grid;
-integer,                intent(in) :: nslope                   ! # of subslope
-real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
-real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded      ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2)
-real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
-
-!   OUTPUT
-real, dimension(ngrid,nslope), intent(inout) :: qsurf                   ! physical point field: Previous and actual density of water ice (kg/m^2)
-real, dimension(ngrid,nslope), intent(inout) :: tendencies_h2o_ice_phys ! physical point field: Evolution of perennial ice over one year (kg/m^2/year)
-logical,                       intent(inout) :: STOPPING                ! Stopping criterion
-
-!   local:
-!   ------
-integer :: i,j,islope                                          ! loop variable
-real    :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve water
-real    ::  new_tendencies(ngrid,nslope)                       ! Tendencies computed in order to conserve water ice on the surface, only exchange between surface are done
-!=======================================================================
-
-  STOPPING=.false.
-
-  pos_tend=0.
-  neg_tend=0.
-if (ngrid.NE.1) then ! to make sure we are not in 1D
-  ! We compute the amount of water accumulating and sublimating
-  do i=1,ngrid
-     if(delta_h2o_adsorbded(i).GT.0) then
-        pos_tend=pos_tend+delta_h2o_adsorbded(i)*cell_area(i)
-     else
-        neg_tend=neg_tend+delta_h2o_adsorbded(i)*cell_area(i)
-     endif
-     if(delta_h2o_icetablesublim(i).GT.0) then
-        pos_tend=pos_tend+delta_h2o_icetablesublim(i)*cell_area(i)
-     else
-        neg_tend=neg_tend+delta_h2o_icetablesublim(i)*cell_area(i)
-     endif
-     do islope=1,nslope
-     if (qsurf(i,islope).GT.0) then
-         if (tendencies_h2o_ice_phys(i,islope).GT.0) then
-            pos_tend=pos_tend+tendencies_h2o_ice_phys(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
-         else
-            neg_tend=neg_tend-tendencies_h2o_ice_phys(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
-         endif
-     endif
-     enddo
-  enddo
-  ! We adapt the tendencies to conserve water and do only exchange between grid points
-   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then ! We are sublimating more in the planet than condensing
-     do i=1,ngrid
-       do islope=1,nslope
-       if(tendencies_h2o_ice_phys(i,islope).LT.0) then ! We lower the sublimating rate by a coefficient
-          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
-       else                                            ! We dont't change the accumulating rate
-          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)
-       endif
-       enddo
-     enddo
-   elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then ! We are condensing more in the planet than sublimating
-     do i=1,ngrid
-       do islope=1,nslope
-       if(tendencies_h2o_ice_phys(i,islope).LT.0) then ! We dont't change the sublimating rate
-          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)
-       else                                            ! We lower the condensing rate by a coefficient
-          new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(neg_tend/pos_tend)
-       endif
-       enddo
-     enddo
-   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then
-    write(*,*) "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
-    write(*,*) "Tendencies on ice sublimating=", neg_tend
-    write(*,*) "Tendencies on ice increasing=", pos_tend
-    write(*,*) "This can be due to the absence of water ice in the PCM run!!"
-    call criterion_waterice_stop(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,qsurf(:,:)*0.)
-    do i=1,ngrid
-       do islope=1,nslope
-         new_tendencies(i,islope)=0
-       enddo
-    enddo
-   endif
-  negative_part = 0.
-
-! Evolution of the water ice for each physical point
-  do i=1,ngrid
-    do islope=1, nslope
-      qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
-      ! We compute the amount of water that is sublimated in excess
-      if (qsurf(i,islope).lt.0) then
-        negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
-        qsurf(i,islope)=0.
-        tendencies_h2o_ice_phys(i,islope)=0.
-      endif
-    enddo
-  enddo
-
-
-  if(pos_tend.eq.0) then
-   real_coefficient = 0.
-  else
-   real_coefficient = negative_part/pos_tend ! We compute a coefficient by which we should remove the ice that has been added
-                                             ! to places even if this ice was contributing to an unphysical negative amount
-                                             ! of ice at other places
-  endif
-  do i=1,ngrid
-    do islope=1, nslope
-     if(new_tendencies(i,islope).GT.0) then  ! In the place of accumulation of ice, we remove a bit of ice in order to conserve water
-         qsurf(i,islope)=qsurf(i,islope)-new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
-     endif
-    enddo
-  enddo
-else ! ngrid==1;
-  do islope=1, nslope
-    qsurf(1,islope)=qsurf(1,islope)+tendencies_h2o_ice_phys(1,islope)*dt_pem
-  enddo
-endif
-
-END SUBROUTINE evol_h2o_ice_s
-
-END MODULE evol_h2o_ice_s_mod
Index: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3149)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3149)
@@ -0,0 +1,160 @@
+MODULE evol_ice_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
+
+use time_evol_mod, only: dt_pem
+
+implicit none
+
+!=======================================================================
+!
+! Routine to compute the evolution of CO2 ice
+!
+!=======================================================================
+!   arguments:
+!   ----------
+!   INPUT
+integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
+
+!   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
+
+!   local:
+!   ----
+
+!=======================================================================
+! Evolution of CO2 ice for each physical point
+write(*,*) 'Evolution of co2 ice'
+co2_ice = co2_ice + tend_co2_ice*dt_pem
+where (co2_ice < 0.)
+    co2_ice = 0.
+    tend_co2_ice = 0.
+end where
+
+END SUBROUTINE evol_co2_ice
+
+!=======================================================================
+
+SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
+
+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
+
+#ifndef CPP_STD
+    use comcstfi_h,   only: pi
+#else
+    use comcstfi_mod, only: pi
+#endif
+
+implicit none
+
+!=======================================================================
+!
+! Routine that compute the evolution of the h2o ice
+!
+!=======================================================================
+!   arguments:
+!   ----------
+!   INPUT
+integer,                intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
+real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
+real, dimension(ngrid), intent(in) :: delta_h2o_adsorbded      ! Mass of H2O adsorbded/desorbded in the soil (kg/m^2)
+real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
+
+!   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)
+integer,                       intent(inout) :: stopPEM      ! Stopping criterion code
+
+!   local:
+!   ------
+integer                       :: i, j, 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
+!=======================================================================
+if (ngrid /= 1) then ! Not in 1D
+    ! We compute the amount of condensing and sublimating h2o ice
+    pos_tend = 0.
+    neg_tend = 0.
+    do i = 1,ngrid
+        if (delta_h2o_adsorbded(i) > 0.) then
+            pos_tend = pos_tend + delta_h2o_adsorbded(i)*cell_area(i)
+        else
+            neg_tend = neg_tend + delta_h2o_adsorbded(i)*cell_area(i)
+        endif
+        if (delta_h2o_icetablesublim(i) > 0.) then
+            pos_tend = pos_tend + delta_h2o_icetablesublim(i)*cell_area(i)
+        else
+            neg_tend = neg_tend + delta_h2o_icetablesublim(i)*cell_area(i)
+        endif
+        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.)
+                else
+                    neg_tend = neg_tend - tend_h2o_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180.)
+                endif
+            endif
+        enddo
+    enddo
+    ! 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
+        elsewhere ! We don't change the condensing rate
+            new_tendencies = tend_h2o_ice
+        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
+        elsewhere ! We lower the condensing rate by a coefficient
+            new_tendencies = tend_h2o_ice*neg_tend/pos_tend
+        end where
+    else 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 h20 ice!"
+        write(*,*) "Tendencies on ice sublimating =", neg_tend
+        write(*,*) "Tendencies on ice increasing =", pos_tend
+        write(*,*) "This can be due to the absence of h2o ice in the PCM run!"
+        stopPEM = 2
+        new_tendencies = 0.
+    endif
+
+    ! Evolution of the h2o ice for each physical point
+    h2o_ice = h2o_ice + new_tendencies*dt_pem
+
+    ! We compute the amount of h2o that is sublimated in excess
+    negative_part = 0.
+    do i = 1,ngrid
+        do islope = 1,nslope
+            if (h2o_ice(i,islope) < 0.) then
+                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.
+            endif
+        enddo
+    enddo
+
+    ! We compute a coefficient by which we should remove the ice that has been added to places
+    ! even if this ice was contributing to an unphysical negative amount of ice at other places
+    if (abs(pos_tend) < 1.e-10) then
+        real_coefficient = 0.
+    else
+        real_coefficient = negative_part/pos_tend
+    endif
+    ! In the place of accumulation of ice, we remove a bit of ice in order to conserve h2o
+    where (new_tendencies > 0.) h2o_ice = h2o_ice - new_tendencies*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
+else ! ngrid == 1, i.e. in 1D
+    h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
+endif
+
+END SUBROUTINE evol_h2o_ice
+
+END MODULE evol_ice_mod
Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3149)
@@ -1,18 +1,14 @@
-module glaciers_mod
-       
- implicit none 
- LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
- LOGICAL  h2oglaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!
-!!! Purpose: Compute CO2 glacier flows
-!!!
-!!! Author: LL
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+MODULE glaciers_mod
+
+implicit none
+
+logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
+logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
+
+!=======================================================================
 contains
-
-subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
+!=======================================================================
+
+SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -26,5 +22,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-IMPLICIT NONE
+implicit none
 
 ! arguments
@@ -32,35 +28,33 @@
 
 ! Inputs:
-      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
-      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
-      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
-      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
-      REAL,INTENT(IN) :: global_ave_ps_GCM          ! Global averaged surface pressure from the GCM [Pa]
-      REAL,INTENT(IN) :: global_ave_ps_PEM          ! global averaged surface pressure during the PEM iteration [Pa]
+      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
+      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
+      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
+      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
+      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
+      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
+      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
      
 ! Ouputs:
-      REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
-      REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
-      REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid)  ! same but within the mesh
-
+      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
+      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
 
 ! Local 
-      REAL :: Tcond(ngrid,nslope) !  Physical field: CO2 condensation temperature [K]
-      REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
+      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
 
 !-----------------------------
-      call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
-
-      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
-
-      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
-   RETURN   
-end subroutine
-
-
-
-
-
-subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
+write(*,*) "Flow of CO2 glacier"
+    
+call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
+call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
+call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
+
+END SUBROUTINE flow_co2glaciers
+
+!=======================================================================
+
+SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -74,6 +68,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-
-IMPLICIT NONE
+implicit none
 
 ! arguments
@@ -81,24 +74,25 @@
 
 ! Inputs:
-      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
-      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
-      REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K]
+      integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
+      real,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
+      real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
 ! Ouputs:
-      REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
-      REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
-      REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
+      real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+      real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
+      real,intent(inout) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
 ! Local 
-      REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
+      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow 
 
 !-----------------------------
-
-      call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
-      call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
-
-   RETURN
-end subroutine
-
-
-subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
+write(*,*) "Flow of H2O glaciers"
+
+call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
+call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
+
+END SUBROUTINE flow_h2oglaciers
+
+!=======================================================================
+
+SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
 
 use abort_pem_mod, only: abort_pem
@@ -118,5 +112,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-   IMPLICIT NONE
+implicit none
 
 ! arguments
@@ -124,16 +118,16 @@
 
 ! Inputs
-      INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
-      INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
-      REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
-      REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
-      character(len=3), INTENT(IN) :: name_ice ! Nature of the ice
+      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
+      integer,intent(in) :: iflat        ! index of the flat subslope
+      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
+      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
+      character(len=3), intent(in) :: name_ice ! Nature of the ice
 ! Outputs
-      REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
+      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
 ! Local
       DOUBLE PRECISION :: tau_d    ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
-      REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
-      INTEGER :: ig,islope ! loop variables
-      REAL :: slo_angle
+      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
+      integer :: ig,islope ! loop variables
+      real :: slo_angle
 
 ! 1. Compute rho
@@ -167,12 +161,9 @@
        ENDDO
     ENDDO
-RETURN
-
-end subroutine
-
-
-
-
-subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
+END SUBROUTINE compute_hmaxglaciers
+
+!=======================================================================
+
+SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -196,20 +187,20 @@
 
 ! Inputs
-      INTEGER, INTENT(IN) :: ngrid,nslope               !# of physical points and subslope
-      INTEGER, INTENT(IN) :: iflat                      ! index of the flat subslope
-      REAL, INTENT(IN) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
-      REAL, INTENT(IN) :: def_slope_mean(nslope)        ! values of the subgrid slopes
-      REAL, INTENT(IN) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
-      REAL, INTENT(IN) :: Tice(ngrid,nslope)            ! Ice temperature[K]
-      character(len=3), INTENT(IN) :: name_ice              ! Nature of the ice
+      integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
+      integer, intent(in) :: iflat                      ! index of the flat subslope
+      real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
+      real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
+      real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
+      real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
+      character(len=3), intent(in) :: name_ice              ! Nature of the ice
 
 ! Outputs
-      REAL, INTENT(INOUT) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
-      REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
-      REAL, INTENT(INOUT) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
+      real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
+      real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
+      real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
 ! Local
-      INTEGER ig,islope ! loop
-      REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
-      INTEGER iaval ! ice will be transfered here
+      integer ig,islope ! loop
+      real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
+      integer iaval ! ice will be transfered here
 
 ! 0. Compute rho
@@ -268,8 +259,9 @@
         ENDDO !islope 
        ENDDO !ig
-RETURN
-end subroutine
-
-     subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
+END SUBROUTINE
+
+!=======================================================================
+
+SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -288,32 +280,27 @@
 
 ! INPUT
-      INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
-      REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
-      REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
-      REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
-      REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
+integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
+real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
+real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
+real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
+
 ! OUTPUT
-      REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged 
+real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged 
 
 ! LOCAL
-
-      INTEGER :: ig,it,islope ! for loop
-      REAL :: ave ! intermediate to compute average
+integer :: ig, it ! For loop
+real    :: ave    ! Intermediate to compute average
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-      DO ig = 1,ngrid
-        ave = 0
-        DO it = 1,timelen
-           ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100))
-        ENDDO
-        DO islope = 1,nslope
-          Tcond(ig,islope) = ave/timelen
-        ENDDO
-      ENDDO
-RETURN
-
-
-end subroutine
-end module
+do ig = 1,ngrid
+    ave = 0
+    do it = 1,timelen
+        ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100))
+    enddo
+    Tcond(ig,:) = ave/timelen
+enddo
+
+END SUBROUTINE computeTcondCO2
+
+END MODULE glaciers_mod
Index: trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3149)
@@ -7,5 +7,5 @@
 !=======================================================================
 
-SUBROUTINE info_PEM(year_iter,criterion_stop,i_myear,n_myear)
+SUBROUTINE info_PEM(year_iter,stopPEM,i_myear,n_myear)
 
 !=======================================================================
@@ -22,5 +22,5 @@
 
 !----- Arguments
-integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
+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
 
@@ -43,5 +43,5 @@
     endif
     open(1,file = 'info_PEM.txt',status = "old",position = "append",action = "write")
-    write(1,*) year_iter, i_myear, criterion_stop
+    write(1,*) year_iter, i_myear, stopPEM
     close(1)
 else
Index: trunk/LMDZ.COMMON/libf/evolution/interpol_TI_PEM2PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/interpol_TI_PEM2PCM_mod.F90	(revision 3149)
+++ trunk/LMDZ.COMMON/libf/evolution/interpol_TI_PEM2PCM_mod.F90	(revision 3149)
@@ -0,0 +1,39 @@
+MODULE interpol_TI_PEM2PCM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE interpol_TI_PEM2PCM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PEM,TI_PCM)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Transfer the thermal inertia from the PEM vertical grid to the PCM vertical grid
+!!!
+!!!
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+
+!======================================================================
+!  arguments
+!  ---------
+!  inputs:
+integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
+integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
+integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
+integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the GCM
+real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
+
+real, dimension(ngrid,nsoil_PCM,nslope), intent(inout) :: TI_PCM ! Thermal inertia in the PCM vertical grid [J/m^2/K/s^{1/2}]
+
+!----- Code
+TI_PCM = TI_PEM(:,:nsoil_PCM,:)
+
+END SUBROUTINE interpol_TI_PEM2PCM
+
+END MODULE interpol_TI_PEM2PCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM_mod.F90	(revision 3145)
+++ 	(revision )
@@ -1,39 +1,0 @@
-MODULE interpolate_TIPEM_TIGCM_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!
-!!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
-!!!
-!!!
-!!! Author: LL
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-implicit none
-
-!======================================================================
-!  arguments
-!  ---------
-!  inputs:
-integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
-integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
-integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
-integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
-real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
-
-real, dimension(ngrid,nsoil_GCM,nslope), intent(inout) :: TI_GCM ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
-
-!----- Code
-TI_GCM(:,:,:) = TI_PEM(:,:nsoil_GCM,:)
-
-END SUBROUTINE interpolate_TIPEM_TIGCM
-
-END MODULE interpolate_TIPEM_TIGCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 3149)
@@ -5,5 +5,5 @@
 !!! 
 !!! 
-!!! Author: LL, inspired by iostart from the GCM
+!!! Author: LL, inspired by iostart from the PCM
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3149)
@@ -13,5 +13,5 @@
 real, dimension(:), allocatable, save :: yearlask   ! year before present from Laskar et al. Tab
 real, dimension(:), allocatable, save :: obllask    ! obliquity    [deg]
-real, dimension(:), allocatable, save :: ecclask    ! excentricity [deg]
+real, dimension(:), allocatable, save :: ecclask    ! eccentricity [deg]
 real, dimension(:), allocatable, save :: 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
Index: trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90	(revision 3149)
@@ -19,5 +19,5 @@
 !=======================================================================
 !
-! Purpose: Read in the data_PCM_Yr*.nc the number of time step
+! Purpose: Read in the data_PCM_Y*.nc the number of time steps
 !
 ! Author: RV
Index: trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3149)
@@ -99,6 +99,6 @@
 
 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.
+max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1.
+min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0.
 write(*,*) 'Maximum eccentricity accepted =', max_ecc
 write(*,*) 'Minimum eccentricity accepted =', min_ecc
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3149)
@@ -13,6 +13,6 @@
 ! II  Run
 !    II_a Update pressure, ice and tracers
-!    II_b Evolution of the ice
-!    II_c CO2 & H2O glaciers flows
+!    II_b Evolution of ice
+!    II_c Flow of glaciers
 !    II_d Update surface and soil temperatures
 !    II_e Outputs
@@ -28,47 +28,45 @@
 PROGRAM pem
 
-use phyetat0_mod,                 only: phyetat0
-use phyredem,                     only: physdem0, physdem1
-use netcdf,                       only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
-use turb_mod,                     only: q2, wstar
-use comslope_mod,                 only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
-use logic_mod,                    only: iflag_phys
-use mod_const_mpi,                only: COMM_LMDZ
+use phyetat0_mod,                only: phyetat0
+use phyredem,                    only: physdem0, physdem1
+use netcdf,                      only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
+use turb_mod,                    only: q2, wstar
+use comslope_mod,                only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
+use logic_mod,                   only: iflag_phys
+use mod_const_mpi,               only: COMM_LMDZ
 use infotrac
-use geometry_mod,                 only: latitude_deg
-use conf_pem_mod,                 only: conf_pem
-use pemredem,                     only: pemdem0, pemdem1
-use glaciers_mod,                 only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
-use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
-use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                        m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial
-use evol_co2_ice_s_mod,           only: evol_co2_ice_s
-use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
-use comsoil_h_PEM,                only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
-                                        TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
-                                        tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
-                                        fluxgeo,                          & ! Geothermal flux for the PEM and PCM
-                                        water_reservoir                     ! Water ressources
-use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
-                                        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 orbit_param_criterion_mod,    only: orbit_param_criterion
-use recomp_orb_param_mod,         only: recomp_orb_param
-use ice_table_mod,                only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                        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
-use abort_pem_mod,                only: abort_pem
-use soil_settings_PEM_mod,        only: soil_settings_PEM
-use compute_tendencies_slope_mod, only: compute_tendencies_slope
-use info_PEM_mod,                 only: info_PEM
-use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
-use nb_time_step_PCM_mod,         only: nb_time_step_PCM
-use pemetat0_mod,                 only: pemetat0
-use read_data_PCM_mod,            only: read_data_PCM
-use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
-use soil_pem_compute_mod,         only: soil_pem_compute
-use writediagpem_mod,             only: writediagpem
+use geometry_mod,                only: latitude_deg
+use conf_pem_mod,                only: conf_pem
+use pemredem,                    only: pemdem0, pemdem1
+use glaciers_mod,                only: flow_co2glaciers, flow_h2oglaciers, co2_ice_flow, h2o_ice_flow
+use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
+use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
+                                       m_noco2, inf_h2oice_threshold
+use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
+use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
+                                       TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
+                                       tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
+                                       fluxgeo                             ! Geothermal flux for the PEM and PCM
+use adsorption_mod,              only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
+                                       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 orbit_param_criterion_mod,   only: orbit_param_criterion
+use recomp_orb_param_mod,        only: recomp_orb_param
+use ice_table_mod,               only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
+                                       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
+use abort_pem_mod,               only: abort_pem
+use soil_settings_PEM_mod,       only: soil_settings_PEM
+use compute_tend_mod,            only: compute_tend
+use info_PEM_mod,                only: info_PEM
+use interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM
+use nb_time_step_PCM_mod,        only: nb_time_step_PCM
+use pemetat0_mod,                only: pemetat0
+use read_data_PCM_mod,           only: read_data_PCM
+use recomp_tend_co2_slope_mod,   only: recomp_tend_co2_slope
+use soil_pem_compute_mod,        only: soil_pem_compute
+use writediagpem_mod,            only: writediagpem
 
 #ifndef CPP_STD
@@ -77,12 +75,11 @@
                                   albedodat, zmea, zstd, zsig, zgam, zthe,     &
                                   hmons, summit, base,albedo_h2o_frost,        &
-                                  frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, &
-                                  emisice, albedice
+                                  frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, emisice, albedice
     use dimradmars_mod,     only: totcloudfrac, albedo
     use dust_param_mod,     only: tauscaling
-    use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
+    use tracer_mod,         only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
     use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
     use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
-    use paleoclimate_mod,   only: albedo_perennialco2
+    use paleoclimate_mod,   only: albedo_perennialco2, paleoclimate
     use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
     use surfini_mod,        only: surfini
@@ -127,7 +124,7 @@
 integer            :: day_ini      ! First day of the simulation
 real               :: pday         ! Physical day
-real               :: time_phys    ! Same as PCM
-real               :: ptimestep    ! Same as PCM
-real               :: ztime_fin    ! Same as PCM
+real               :: time_phys    ! Same as in PCM
+real               :: ptimestep    ! Same as in PCM
+real               :: ztime_fin    ! Same as in PCM
 
 ! Variables to read start.nc
@@ -140,6 +137,6 @@
 real, dimension(:,:,:), allocatable :: q             ! champs advectes
 real, dimension(ip1jmp1)            :: ps            ! pression au sol
-real, dimension(:),     allocatable :: ps_start_GCM  ! (ngrid) pression  au sol
-real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées
+real, dimension(:),     allocatable :: ps_start_PCM  ! (ngrid) surface pressure
+real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
 real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
 real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
@@ -161,16 +158,16 @@
 
 ! Variables for h2o_ice evolution
-real                                :: ini_surf_h2o         ! Initial surface of sublimating h2o ice
-real                                :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
-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 PCM [kg/m^2]
+real, dimension(:,:),   allocatable :: h2o_ice              ! h2o ice in the PEM
+real, dimension(:,:),   allocatable :: h2o_ice_ini          ! Initial amount of h2o ice in the PEM
+real, dimension(:,:,:), allocatable :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
+real                                :: h2o_surf_ini         ! Initial surface of sublimating h2o ice
+real                                :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
+real                                :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
+real                                :: global_avg_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_PCM            ! same but retrieved from the PCM [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?
-logical                             :: STOPPING_1_water     ! Logical: is there still water ice to sublimate?
-logical                             :: STOPPING_pressure    ! Logical: is the criterion (% of change in the surface pressure) reached?
-integer                             :: criterion_stop       ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
+integer                             :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
 real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
 real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
@@ -181,44 +178,39 @@
 
 ! Variables for co2_ice evolution
-real                              :: ini_surf_co2         ! Initial surface of sublimating co2 ice
-logical                           :: STOPPING_co2         ! Logical: is the criterion (% of change in the surface of sublimating co2 ice) reached?
-real, dimension(:,:), allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
-real, dimension(:,:), allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
-real, dimension(:,:), allocatable :: q_co2_PEM_phys       ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
-real, dimension(:,:), allocatable :: co2frost_ini         ! Initial amount of co2 frost (at the start of the PEM run)
-real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
+real, dimension(:,:),   allocatable :: co2_ice          ! co2 ice in the PEM
+real, dimension(:,:),   allocatable :: co2_ice_ini      ! Initial amount of co2 ice in the PEM
+real, dimension(:,:,:), allocatable :: min_co2_ice      ! Minimum of co2 ice at each point for the first year [kg/m^2]
+real                                :: co2_surf_ini     ! Initial surface of sublimating co2 ice
+real, dimension(:,:),   allocatable :: vmr_co2_PCM      ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real, dimension(:,:),   allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
+real, dimension(:,:),   allocatable :: q_co2_PEM_phys   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
 
 ! Variables for slopes
-real, dimension(:,:),   allocatable :: min_co2_ice_1          ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2]
-real, dimension(:,:),   allocatable :: min_co2_ice_2          ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2]
-real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field: minimum of water ice at each point for the first year [kg/m^2]
-real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field: minimum of water ice at each point for the second year [kg/m^2]
-real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field: co2 ice given by the PCM  [kg/m^2]
-real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
-real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field: Logical array indicating if there is water ice at initial state
-real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field: Logical array indicating if there is co2 ice at initial state
-real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
-real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
-real, dimension(:,:),   allocatable :: tendencies_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 :: co2_ice_PCM        ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
+real, dimension(:,:),   allocatable :: co2_ice_ini_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
+real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is water ice at initial state
+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 x nslope: 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
 
 ! Variables for surface and soil
-real, dimension(:,:),     allocatable :: tsurf_ave                          ! Physic x SLOPE field : Averaged Surface Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsoil_ave                          ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsurf_GCM_timeseries               ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries               ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-real, dimension(:,:),     allocatable :: tsurf_ave_yr1                      ! Physic x SLOPE field : Averaged Surface Temperature of first call of the PCM [K]
+real, dimension(:,:),     allocatable :: tsurf_ave                          ! Physic x SLOPE field: Averaged Surface Temperature [K]
+real, dimension(:,:,:),   allocatable :: tsoil_ave                          ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
+real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries               ! ngrid x SLOPE XTULES field: Surface Temperature in timeseries [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_PCM_timeseries               ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
+real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
 real, dimension(:,:),     allocatable :: TI_locslope                        ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
 real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: intermediate when computing Tsoil [K]
-real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
+real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
 real, dimension(:,:),     allocatable :: watersurf_density_ave              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
 real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_ave          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
-real, dimension(:,:),     allocatable :: Tsurfave_before_saved              ! Surface temperature saved from previous time step [K]
+real, dimension(:,:),     allocatable :: Tsurfavg_before_saved              ! Surface temperature saved from previous time step [K]
 real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
 real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
@@ -236,6 +228,4 @@
 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]
-real            :: water_sum     ! total mass of water in the mesh [kg/m^2]
 
 #ifdef CPP_STD
@@ -341,9 +331,9 @@
 !------------------------
 ! I_b.1 Read start_evol.nc
-allocate(ps_start_GCM(ngrid))
+allocate(ps_start_PCM(ngrid))
 #ifndef CPP_1D
     call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
 
-    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
+    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
 
     call iniconst !new
@@ -359,8 +349,7 @@
     status = nf90_close(ncid)
 
-    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
-                   rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
 #else
-    ps_start_GCM(1) = ps(1)
+    ps_start_PCM(1) = ps(1)
 #endif
 
@@ -387,5 +376,5 @@
 #ifndef CPP_STD
     call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
-                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
+                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
                   watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
 
@@ -414,6 +403,6 @@
     allocate(inertiesoil(ngrid,nsoilmx,1))
     call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
-                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,            &
-                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,     &
+                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
+                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
                   rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
@@ -422,13 +411,13 @@
     call ini_comslope_h(ngrid,1)
 
-    qsurf(:,:,1) = qsurf_read_generic(:,:)
-    tsurf(:,1) = tsurf_read_generic(:)
-    tsoil(:,:,1) = tsoil_read_generic(:,:)
-    emis(:,1) = emis_read_generic(:)
+    qsurf(:,:,1) = qsurf_read_generic
+    tsurf(:,1) = tsurf_read_generic
+    tsoil(:,:,1) = tsoil_read_generic
+    emis(:,1) = emis_read_generic
     watercap(:,1) = 0.
     watercaptag(:) = .false.
     albedo(:,1,1) = albedo_read_generic(:,1)
     albedo(:,2,1) = albedo_read_generic(:,2)
-    inertiesoil(:,:,1) = inertiedat(:,:)
+    inertiesoil(:,:,1) = inertiedat
 
     if (nslope == 1) then
@@ -440,5 +429,5 @@
 
     ! Remove unphysical values of surface tracer
-    qsurf(:,:,1) = qsurf_read_generic(:,:)
+    qsurf(:,:,1) = qsurf_read_generic
     where (qsurf < 0.) qsurf = 0.
 #endif
@@ -472,8 +461,8 @@
 allocate(flag_h2oflow_mesh(ngrid))
 
-flag_co2flow(:,:) = 0
-flag_co2flow_mesh(:) = 0
-flag_h2oflow(:,:) = 0
-flag_h2oflow_mesh(:) = 0
+flag_co2flow = 0
+flag_co2flow_mesh = 0
+flag_h2oflow = 0
+flag_h2oflow_mesh = 0
 
 !------------------------
@@ -486,20 +475,18 @@
 allocate(tsoil_ave(ngrid,nsoilmx,nslope))
 allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
-allocate(vmr_co2_gcm(ngrid,timelen))
+allocate(vmr_co2_PCM(ngrid,timelen))
 allocate(ps_timeseries(ngrid,timelen))
-allocate(min_co2_ice_1(ngrid,nslope))
-allocate(min_h2o_ice_1(ngrid,nslope))
-allocate(min_co2_ice_2(ngrid,nslope))
-allocate(min_h2o_ice_2(ngrid,nslope))
-allocate(tsurf_ave_yr1(ngrid,nslope))
+allocate(min_co2_ice(ngrid,nslope,2))
+allocate(min_h2o_ice(ngrid,nslope,2))
+allocate(tsurf_avg_yr1(ngrid,nslope))
 allocate(tsurf_ave(ngrid,nslope))
-allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
-allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
+allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
+allocate(tsoil_PCM_timeseries(ngrid,nsoilmx,nslope,timelen))
 allocate(q_co2_PEM_phys(ngrid,timelen))
 allocate(q_h2o_PEM_phys(ngrid,timelen))
-allocate(co2_ice_GCM(ngrid,nslope,timelen))
+allocate(co2_ice_PCM(ngrid,nslope,timelen))
 allocate(watersurf_density_ave(ngrid,nslope))
 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
-allocate(Tsurfave_before_saved(ngrid,nslope))
+allocate(Tsurfavg_before_saved(ngrid,nslope))
 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
@@ -508,17 +495,17 @@
 allocate(delta_h2o_icetablesublim(ngrid))
 allocate(delta_h2o_adsorbded(ngrid))
-allocate(vmr_co2_pem_phys(ngrid,timelen))
+allocate(vmr_co2_PEM_phys(ngrid,timelen))
 
 write(*,*) "Downloading data Y1..."
-call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
-                   tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
-                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
+call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
+                   tsurf_avg_yr1,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                      &
+                   co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
 write(*,*) "Downloading data Y1 done"
 
 ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
 write(*,*) "Downloading data Y2"
-call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
-                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
-                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
+call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), &
+                   tsurf_ave,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                         &
+                   co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
 write(*,*) "Downloading data Y2 done"
 
@@ -536,16 +523,14 @@
 if (soil_pem) then
     call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
-    do l = 1,nsoilmx
-        tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
-        tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
-        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
-    enddo
+    tsoil_PEM(:,1:nsoilmx,:) = tsoil_ave
+    tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_PCM_timeseries
+    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
     do l = nsoilmx + 1,nsoilmx_PEM
         tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
         watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
     enddo
-    watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+    watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
 endif !soil_pem
-deallocate(tsoil_ave,tsoil_GCM_timeseries)
+deallocate(tsoil_ave,tsoil_PCM_timeseries)
 
 !------------------------
@@ -553,14 +538,11 @@
 !    I_f Compute tendencies
 !------------------------
-allocate(tendencies_co2_ice(ngrid,nslope))
-allocate(tendencies_co2_ice_ini(ngrid,nslope))
-allocate(tendencies_h2o_ice(ngrid,nslope))
+allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
+allocate(tend_co2_ice_ini(ngrid,nslope))
 
 ! Compute the tendencies of the evolution of ice over the years
-call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
-tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
-call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
-
-deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
+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
 
 !------------------------
@@ -568,53 +550,30 @@
 !    I_g Save initial PCM situation
 !------------------------
-allocate(initial_co2_ice_sublim(ngrid,nslope))
-allocate(initial_co2_ice(ngrid,nslope))
-allocate(initial_h2o_ice(ngrid,nslope))
-
 ! We save the places where water ice is sublimating
 ! We compute the surface of water ice sublimating
-ini_surf_co2 = 0.
-ini_surf_h2o = 0.
+co2_surf_ini = 0.
+h2o_surf_ini = 0.
 Total_surface = 0.
 do i = 1,ngrid
     Total_surface = Total_surface + cell_area(i)
     do islope = 1,nslope
-        if (tendencies_co2_ice(i,islope) < 0) then
-            initial_co2_ice_sublim(i,islope) = 1.
-            ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
-        else
-            initial_co2_ice_sublim(i,islope) = 0.
-        endif
-        if (qsurf(i,igcm_co2,islope) > 0) then
-            initial_co2_ice(i,islope) = 1.
-        else
-            initial_co2_ice(i,islope) = 0.
-        endif
-        if (tendencies_h2o_ice(i,islope) < 0) then
-            initial_h2o_ice(i,islope) = 1.
-            ini_surf_h2o = ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
-        else
-            initial_h2o_ice(i,islope) = 0.
-        endif
+        if (tend_co2_ice(i,islope) < 0.) co2_surf_ini = co2_surf_ini + cell_area(i)*subslope_dist(i,islope)
+        if (tend_h2o_ice(i,islope) < 0.) h2o_surf_ini = h2o_surf_ini + cell_area(i)*subslope_dist(i,islope)
     enddo
 enddo
 
-write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
-write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
+write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2_surf_ini
+write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2o_surf_ini
 write(*,*) "Total surface of the planet =", Total_surface
-allocate(zplev_gcm(ngrid,nlayer + 1))
+allocate(zplev_PCM(ngrid,nlayer + 1))
 
 do ig = 1,ngrid
-    zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
+    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
 enddo
 
-global_ave_press_old = 0.
-do i = 1,ngrid
-    global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
-enddo
-
-global_ave_press_GCM = global_ave_press_old
-global_ave_press_new = global_ave_press_old
-write(*,*) "Initial global average pressure =", global_ave_press_GCM
+global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
+global_avg_press_PCM = global_avg_press_old
+global_avg_press_new = global_avg_press_old
+write(*,*) "Initial global average pressure =", global_avg_press_PCM
 
 !------------------------
@@ -622,23 +581,17 @@
 !    I_h Read the startpem.nc
 !------------------------
+allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
+deallocate(min_co2_ice,min_h2o_ice)
+
 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
-              porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,       &
-              tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                   &
-              qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,          &
-              co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
-
-delta_h2o_icetablesublim(:) = 0.
-
-! We save the initial values for the co2 frost and perennial ice
-allocate(co2frost_ini(ngrid,nslope),perennial_co2ice_ini(ngrid,nslope))
-co2frost_ini = qsurf(:,igcm_co2,:)
-perennial_co2ice_ini = perennial_co2ice
-
-do ig = 1,ngrid
-    do islope = 1,nslope
-        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
-        qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
-    enddo
-enddo
+              porefillingice_thickness,tsurf_avg_yr1,tsurf_ave,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_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded,                &
+              h2o_adsorbded_phys,delta_h2o_adsorbded)
+
+allocate(co2_ice_ini(ngrid,nslope))
+co2_ice_ini = co2_ice
+
+delta_h2o_icetablesublim = 0.
 
 if (adsorption_pem) then
@@ -648,16 +601,15 @@
         do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
-                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
             enddo
         enddo
     enddo
-
     write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
     write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
 endif ! adsorption
-deallocate(tsurf_ave_yr1)
+deallocate(tsurf_avg_yr1)
 
 !------------------------
@@ -684,4 +636,5 @@
 !------------------------
 year_iter = 0
+stopPEM = 0
 
 do while (year_iter < year_iter_max .and. i_myear < n_myear)
@@ -690,5 +643,5 @@
     do i = 1,ngrid
         do islope = 1,nslope
-            global_ave_press_new = global_ave_press_new - g*cell_area(i)*tendencies_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 - g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
         enddo
     enddo
@@ -696,9 +649,9 @@
     if (adsorption_pem) then
         do i = 1,ngrid
-            global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
+            global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
         enddo
     endif
-    write(*,*) 'Global average pressure old time step',global_ave_press_old
-    write(*,*) 'Global average pressure new time step',global_ave_press_new
+    write(*,*) 'Global average pressure old time step',global_avg_press_old
+    write(*,*) 'Global average pressure new time step',global_avg_press_new
 
 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
@@ -714,9 +667,9 @@
     write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
     do ig = 1,ngrid
-        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
+        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
     enddo
 
 ! II.a.4. New pressure levels timeseries
-    allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
+    allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
     write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
     do l = 1,nlayer + 1
@@ -755,5 +708,5 @@
             endif
             mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
-            vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
+            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
         enddo
     enddo
@@ -763,22 +716,16 @@
 !------------------------
 ! II  Run
-!    II_b Evolution of the ice
-!------------------------
-    write(*,*) "Evolution of h2o ice"
-    call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water)
-
-    write(*,*) "Evolution of co2 ice"
-    call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
+!    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)
 
 !------------------------
 ! II  Run
-!    II_c CO2 & H2O glaciers flows
-!------------------------
-    write(*,*) "CO2 glacier flows"
-    if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
-                                               global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
-
-    write(*,*) "H2O glacier flows"
-    if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
+!    II_c Flow of glaciers
+!------------------------
+    if (co2_ice_flow) 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)
+    if (h2o_ice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
 
 !------------------------
@@ -789,12 +736,12 @@
     write(*,*) "Updating the new Tsurf"
     bool_sublim = .false.
-    Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
+    Tsurfavg_before_saved = tsurf_ave
     do ig = 1,ngrid
         do islope = 1,nslope
-            if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice
+            if (co2_ice_ini(ig,islope) > 0.5 .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice
                 if (latitude_deg(ig) > 0) then
                     do ig_loop = ig,ngrid
                         do islope_loop = islope,iflat,-1
-                            if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
+                            if (co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
                                 bool_sublim = .true.
@@ -807,5 +754,5 @@
                     do ig_loop = ig,1,-1
                         do islope_loop = islope,iflat
-                            if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
+                            if(co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
                                 bool_sublim = .true.
@@ -816,6 +763,6 @@
                     enddo
                 endif
-                initial_co2_ice(ig,islope) = 0
-                if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
+                co2_ice_ini(ig,islope) = 0
+                if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
                     albedo(ig,1,islope) = albedo_h2o_frost
                     albedo(ig,2,islope) = albedo_h2o_frost
@@ -825,11 +772,11 @@
                     emis(ig,islope) = emissiv
                 endif
-            else if ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2
+            else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
                 ave = 0.
                 do t = 1,timelen
-                    if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
-                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
+                    if (co2_ice_PCM(ig,islope,t) > 1.e-3) then
+                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
                     else
-                        ave = ave + tsurf_GCM_timeseries(ig,islope,t)
+                        ave = ave + tsurf_PCM_timeseries(ig,islope,t)
                     endif
                 enddo
@@ -849,10 +796,10 @@
 
     do t = 1,timelen
-        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
+        tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_ave - Tsurfavg_before_saved
     enddo
     ! for the start
     do ig = 1,ngrid
         do islope = 1,nslope
-            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
+            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_ave(ig,islope))
         enddo
     enddo
@@ -868,11 +815,11 @@
         ! Soil averaged
         do islope = 1,nslope
-            TI_locslope(:,:) = TI_PEM(:,:,islope)
+            TI_locslope = TI_PEM(:,:,islope)
             do t = 1,timelen
-                Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
-                Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
+                Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t)
+                Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t)
                 call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
                 call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
-                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
+                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope
                 do ig = 1,ngrid
                     do isoil = 1,nsoilmx_PEM
@@ -883,6 +830,6 @@
             enddo
         enddo
-        tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
-        watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+        tsoil_PEM = sum(tsoil_phys_PEM_timeseries,4)/timelen
+        watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
 
         write(*,*) "Update of soil temperature done"
@@ -892,17 +839,16 @@
 ! II_d.3 Update the ice table
         write(*,*) "Compute ice table"
-        porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
+        porefillingice_thickness_prev_iter = porefillingice_thickness
         call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
         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
 
 ! II_d.4 Update the soil thermal properties
-        write(*,*) "Update soil propreties"
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_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,tendencies_h2o_ice,tendencies_co2_ice, &
-                                     qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),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)
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,       &
+                                     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)
 
             totmassco2_adsorbded = 0.
@@ -911,7 +857,7 @@
                 do islope =1, nslope
                     do l = 1,nsoilmx_PEM - 1
-                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
                         subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
                         subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                     enddo
@@ -927,11 +873,11 @@
 !    II_e Outputs
 !------------------------
-    call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
+    call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
     do islope = 1,nslope
         write(str2(1:2),'(i2.2)') islope
-        call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
-        call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
-        call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
-        call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
+        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,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
         call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
@@ -942,7 +888,6 @@
 !    II_f Update the tendencies
 !------------------------
-    write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
-    call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
-                               global_ave_press_GCM,global_ave_press_new)
+    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, &
+                               global_avg_press_PCM,global_avg_press_new)
 
 !------------------------
@@ -950,39 +895,30 @@
 !    II_g Checking the stopping criterion
 !------------------------
-    call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
-
-    call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
-                            initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
+    call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
+    call stopping_crit_co2(cell_area,co2_surf_ini,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
 
-    write(*,*) "Checking all the stopping criteria..."
-    if (STOPPING_water) then
-        write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
-        criterion_stop = 1
-    endif
-    if (STOPPING_1_water) then
-        write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
-        criterion_stop = 1
-    endif
-    if (STOPPING_co2) then
-        write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
-        criterion_stop = 2
-    endif
-    if (STOPPING_pressure) then
-        write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
-        criterion_stop = 3
-    endif
-    if (year_iter >= year_iter_max) then
-        write(*,*) "STOPPING because maximum number of iterations reached"
-        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
+    write(*,*) "Checking the stopping criteria..."
+    if (year_iter >= year_iter_max) stopPEM = 5
+    if (i_myear >= n_myear) stopPEM = 6
+    if (stopPEM > 0) then
+        select case (stopPEM)
+            case(1)
+                write(*,*) "STOPPING because surface of water ice sublimating is too low:", stopPEM, "See message above."
+            case(2)
+                write(*,*) "STOPPING because tendencies on water ice = 0:", stopPEM, "See message above."
+            case(3)
+                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
+            case(4)
+                write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
+            case(5)
+                write(*,*) "STOPPING because maximum number of iterations due to orbital parameters is reached:", stopPEM
+            case(6)
+                write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
+            case default
+                write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
+        end select
         exit
     else
@@ -992,7 +928,7 @@
     endif
 
-    global_ave_press_old = global_ave_press_new
-
-enddo  ! big time iteration loop of the pem
+    global_avg_press_old = global_avg_press_new
+
+enddo ! big time iteration loop of the pem
 !------------------------------ END RUN --------------------------------
 
@@ -1005,84 +941,27 @@
 
 ! H2O ice
+watercap = 0.
 do ig = 1,ngrid
-    if (watercaptag(ig)) then
-        watercap_sum = 0.
-        do islope = 1,nslope
-            ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not.
-            if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then
-                ! 1st case: more ice than initialy. Put back ancien frost
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
-            else
-                ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap
-                watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope)
-                qsurf(ig,igcm_h2o_ice,islope) = 0.
-            endif
-            watercap_sum = watercap_sum + watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-            watercap(ig,islope) = 0.
-        enddo
-        water_reservoir(ig) = water_reservoir(ig) + watercap_sum
+    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
+    ! If there is enough ice to be considered as an infinite reservoir
+        watercaptag(ig) = .true.
+    else
+    ! If there too little ice to be considered as an infinite reservoir
+    ! Ice is transferred to the frost
+        watercaptag(ig) = .false.
+        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
+        h2o_ice(ig,:) = 0.
     endif
 enddo
 
-do ig = 1,ngrid
-    water_sum = 0.
-    do islope = 1,nslope
-        water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-    enddo
-    if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
-        if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
-            watercaptag(ig) = .true.
-            water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
-            do islope = 1,nslope
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
-            enddo
-        endif
-    else ! let's check that the infinite source of water disapear
-        if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then
-            watercaptag(ig) = .false.
-            do islope = 1,nslope
-                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
-            enddo
-            water_reservoir(ig) = 0.
-        endif
-    endif
-enddo
-
-! CO2 ice
-do ig = 1,ngrid
-    do islope = 1,nslope
-        if (qsurf(ig,igcm_co2,islope) == threshold_co2_frost2perennial) then
-        ! If co2 ice is equal to the threshold, then everything is transformed into perennial ice
-            perennial_co2ice(ig,islope) = threshold_co2_frost2perennial
-            qsurf(ig,igcm_co2,islope) = 0.
-        else if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perennial) then
-        ! If co2 ice is superior to the threshold, then co2 frost is equal to the threshold (max possible value)
-        ! and the leftover is transformed into perennial ice
-            perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) - threshold_co2_frost2perennial
-            qsurf(ig,igcm_co2,islope) = threshold_co2_frost2perennial
-        else
-        ! If co2 ice is inferior to the threshold, then we compare with the initial state
-            if (qsurf(ig,igcm_co2,islope) > perennial_co2ice_ini(ig,islope)) then
-            ! If co2 ice higher than the initial perennial ice, then the change is affected only to the frost
-                perennial_co2ice(ig,islope) = perennial_co2ice_ini(ig,islope)
-                qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) - perennial_co2ice_ini(ig,islope)
-            else
-            ! If co2 ice is lower than the initial perennial ice, then there is no frost anymore
-                perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope)
-                qsurf(ig,igcm_co2,islope) = 0.
-            endif
-        endif
-    enddo
-enddo
-
 ! III_a.2  Tsoil update (for startfi)
 if (soil_pem) then
-    call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
-    tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
+    call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
+    tsoil = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
 endif
 
 ! III_a.4 Pressure (for start)
-ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
-ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
+ps = ps*global_avg_press_new/global_avg_press_PCM
+ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
 
 ! III_a.5 Tracer (for start)
@@ -1090,5 +969,5 @@
 
 do l = 1,nlayer + 1
-    zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
+    zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
 enddo
 
@@ -1097,5 +976,5 @@
         do l = 1,llm - 1
             do ig = 1,ngrid
-                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
             enddo
             q(:,llm,nnq) = q(:,llm - 1,nnq)
@@ -1104,6 +983,6 @@
         do l = 1,llm - 1
             do ig = 1,ngrid
-                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
-                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
+                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
             enddo
             q(:,llm,nnq) = q(:,llm - 1,nnq)
@@ -1178,7 +1057,8 @@
 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
              TI_PEM, porefillingice_depth,porefillingice_thickness,      &
-             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
+             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice)
 write(*,*) "restartpem.nc has been written"
-call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
+
+call info_PEM(year_iter,stopPEM,i_myear,n_myear)
 
 write(*,*) "The PEM has run for", year_iter, "Martian years."
@@ -1187,9 +1067,9 @@
 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)
-deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
+deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
+deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
-deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
-deallocate(co2frost_ini,perennial_co2ice_ini)
+deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
+deallocate(co2_ice_ini)
 !----------------------------- END OUTPUT ------------------------------
 
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3149)
@@ -7,12 +7,12 @@
 !=======================================================================
 
-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)
+SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, &
+                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,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)
 
 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 comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_h2o_bigreservoir, 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
@@ -26,5 +26,5 @@
 #ifndef CPP_STD
     use comcstfi_h,   only: r, mugaz
-    use surfdat_h,    only: watercaptag
+    use surfdat_h,    only: watercaptag, perennial_co2ice
 #else
     use comcstfi_mod, only: r, mugaz
@@ -37,23 +37,23 @@
 character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
 integer,                        intent(in) :: ngrid               ! # of physical grid points
-integer,                        intent(in) :: nsoil_GCM           ! # of vertical grid points in the GCM
+integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
 integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
 integer,                        intent(in) :: nslope              ! # of sub-grid slopes
 integer,                        intent(in) :: timelen             ! # time samples
 real,                           intent(in) :: timestep            ! time step [s]
-real,                           intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa]
-real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
-real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM call [K]
+real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
 real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
 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_h2oglaciers    ! tendencies on h2o glaciers
-real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
-real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
-real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
-real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
+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) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! outputs
 real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
 real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
+real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
@@ -63,6 +63,5 @@
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
-real, dimension(ngrid),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! surface water ice density, yearly averaged (kg/m^3)
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! local
 real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
@@ -75,5 +74,5 @@
 character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
 real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr1[K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr 1 [K]
 real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
 real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
@@ -82,6 +81,6 @@
 
 #ifdef CPP_STD
-   logical, dimension(ngrid) :: watercaptag
-   watercaptag(:) = .false.
+    logical, dimension(ngrid) :: watercaptag
+    watercaptag(:) = .false.
 #endif
 
@@ -91,9 +90,9 @@
 !!!
 !!! Order: 0. Previous year of the PEM run
+!!!           Ice initialization
 !!!        1. Thermal Inertia
 !!!        2. Soil Temperature
 !!!        3. Ice table
 !!!        4. Mass of CO2 & H2O adsorbed
-!!!        5. Water reservoir
 !!!
 !!! /!\ This order must be respected !
@@ -113,371 +112,352 @@
     call open_startphy(filename)
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+#ifndef CPP_STD
+    ! h2o ice
+    h2o_ice = 0.
+    call get_field("h2o_ice",h2o_ice,found)
+    if (.not. found) then
+        write(*,*)'Pemetat0: failed loading <h2o_ice>'
+        write(*,*)'will reconstruct the values from watercaptag'
+        write(*,*)'with default value ''ini_h2o_bigreservoir'''
+        do ig = 1,ngrid
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
+        enddo
+    endif
+
+    ! co2 ice
+    co2_ice = perennial_co2ice
+#endif
+
     if (soil_pem) then
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !1. Thermal Inertia
 ! a. General case
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
+            if (.not. found) then
+                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
+                write(*,*)'will reconstruct the values of TI_PEM'
+
+                do ig = 1,ngrid
+                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
+                        !!! transition
+                        delta = depth_breccia
+                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
+                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
+                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
+                        do iloop=index_breccia+2,index_bedrock
+                            TI_PEM(ig,iloop,islope) = TI_breccia
+                        enddo
+                    else ! we keep the high ti values
+                        do iloop = index_breccia + 1,index_bedrock
+                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
+                        enddo
+                    endif ! TI PEM and breccia comparison
+                    !! transition
+                    delta = depth_bedrock
+                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
+                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
+                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
+                    do iloop = index_bedrock + 2,nsoil_PEM
+                        TI_PEM(ig,iloop,islope) = TI_bedrock
+                    enddo
+                enddo
+            else ! found
+                do iloop = nsoil_PCM + 1,nsoil_PEM
+                    TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
+                enddo
+            endif ! not found
+        enddo ! islope
+
+        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
+
+! b. Special case for inertiedat, inertiedat_PEM
+        call get_field("inertiedat_PEM",inertiedat_PEM,found)
+        if (.not.found) then
+            do iloop = 1,nsoil_PCM
+                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
+            enddo
+        !!! zone de transition
+            delta = depth_breccia
+            do ig = 1,ngrid
+                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
+                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
+                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+
+                    do iloop = index_breccia+2,index_bedrock
+                        inertiedat_PEM(ig,iloop) = TI_breccia
+                    enddo
+                else
+                    do iloop=index_breccia+1,index_bedrock
+                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
+                    enddo
+                endif ! comparison ti breccia
+            enddo ! ig
+
+            !!! zone de transition
+            delta = depth_bedrock
+            do ig = 1,ngrid
+                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
+                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
+            enddo
+
+            do iloop = index_bedrock + 2, nsoil_PEM
+                do ig = 1,ngrid
+                    inertiedat_PEM(ig,iloop) = TI_bedrock
+                enddo
+            enddo
+        endif ! not found
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!2. Soil Temperature
+        do islope=1,nslope
+            write(num,fmt='(i2.2)') islope
+            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
+            if (.not.found) then
+                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
+                write(*,*)'will reconstruct the values of Tsoil'
+!                do ig = 1,ngrid
+!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
+!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
+!                    do iloop=index_breccia+2,index_bedrock
+!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
+!                    enddo
+!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
+!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
+!
+!                    do iloop=index_bedrock+2,nsoil_PEM
+!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
+!                    enddo
+!                enddo
+                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+            else
+! predictor corrector: restart from year 1 of the PCM and build the evolution of
+! tsoil at depth
+                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
+                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+                call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
+                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
+
+                do iloop = nsoil_PCM+1,nsoil_PEM
+                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
+                enddo
+            endif !found
+
+            do it = 1,timelen
+                do isoil = nsoil_PCM+1,nsoil_PEM
+                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
+                enddo
+            enddo
+            do isoil = nsoil_PCM+1,nsoil_PEM
+                do ig = 1,ngrid
+                    watersoil_avg(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 ! islope
+        write(*,*) 'PEMETAT0: SOIL TEMP done'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!3. Ice Table
+        call get_field("ice_table",ice_table,found)
         if (.not. found) then
-            write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
-            write(*,*)'will reconstruct the values of TI_PEM'
-
+            write(*,*)'PEM settings: failed loading <ice_table>'
+            write(*,*)'will reconstruct the values of the ice table given the current state'
+            ice_table = -1  ! by default, no ice table
+            ice_table_thickness = -1  ! by default, no ice table
+            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
+            do islope = 1,nslope
+                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+            enddo
+        endif
+
+        write(*,*) 'PEMETAT0: ICE TABLE done'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!4. CO2 & H2O Adsorption
+        if (adsorption_pem) then
+            do islope = 1,nslope
+                write(num,fmt='(i2.2)') islope
+                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
+                if (.not. found) then
+                    m_co2_regolith_phys = 0.
+                    exit
+                endif
+            enddo
+            do islope=1,nslope
+                write(num,fmt='(i2.2)') islope
+                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
+                if (.not. found2) then
+                    m_h2o_regolith_phys = 0.
+                    exit
+                endif
+            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, &
+                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
+
+            if (.not. found) deltam_co2_regolith_phys = 0.
+            if (.not.found2) deltam_h2o_regolith_phys = 0.
+
+            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
+        endif ! adsorption_pem
+    endif ! soil_pem
+
+    call close_startphy
+
+else !No startfi, let's build all by hand
+
+    ! h2o ice
+    h2o_ice = 0.
+    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_h2o_bigreservoir''.'
+        do ig = 1,ngrid
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
+        enddo
+
+    ! co2 ice
+    co2_ice = perennial_co2ice
+
+    if (soil_pem) then
+
+!a) Thermal inertia
+        do islope = 1,nslope
             do ig = 1,ngrid
                 if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
                     !!! transition
                     delta = depth_breccia
-                    TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
-                                                        (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
-                                                        ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
-                    do iloop=index_breccia+2,index_bedrock
+                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
+                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
+                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
+                    do iloop=index_breccia + 2,index_bedrock
                         TI_PEM(ig,iloop,islope) = TI_breccia
                     enddo
                 else ! we keep the high ti values
-                    do iloop = index_breccia + 1,index_bedrock
+                    do iloop=index_breccia + 1,index_bedrock
                         TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
                     enddo
-                endif ! TI PEM and breccia comparison
+                endif
                 !! transition
                 delta = depth_bedrock
-                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
-                                                        (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
-                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
+                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
+                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
+                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
                 do iloop = index_bedrock + 2,nsoil_PEM
                     TI_PEM(ig,iloop,islope) = TI_bedrock
                 enddo
             enddo
-        else ! found
-            do iloop = nsoil_GCM + 1,nsoil_PEM
-                TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
-            enddo
-        endif ! not found
-    enddo ! islope
-
-    write(*,*) 'PEMETAT0: THERMAL INERTIA done'
-
-! b. Special case for inertiedat, inertiedat_PEM
-call get_field("inertiedat_PEM",inertiedat_PEM,found)
-if(.not.found) then
- do iloop = 1,nsoil_GCM
-   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
- enddo
-!!! zone de transition
-delta = depth_breccia
-do ig = 1,ngrid
-if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
-inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
-            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
-                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
-
- do iloop = index_breccia+2,index_bedrock
-       inertiedat_PEM(ig,iloop) = TI_breccia
-  enddo
-
-else
-   do iloop=index_breccia+1,index_bedrock
-      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
-   enddo
-endif ! comparison ti breccia
-enddo!ig
-
-!!! zone de transition
-delta = depth_bedrock
-do ig = 1,ngrid
-inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
-                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
-enddo
-
- do iloop = index_bedrock+2, nsoil_PEM
-   do ig = 1,ngrid
-      inertiedat_PEM(ig,iloop) = TI_bedrock
-   enddo
- enddo
-endif ! not found
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!2. Soil Temperature
-do islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
-  if(.not.found) then
-      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
-      write(*,*)'will reconstruct the values of Tsoil'
-!      do ig = 1,ngrid
-!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
-!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
-!       do iloop=index_breccia+2,index_bedrock
-!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
-!      enddo
-!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
-!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
+        enddo
+
+        do iloop = 1,nsoil_PCM
+            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
+        enddo
+        !!! zone de transition
+        delta = depth_breccia
+        do ig = 1,ngrid
+            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
+                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
+                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
+                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
+                do iloop = index_breccia + 2,index_bedrock
+                    inertiedat_PEM(ig,iloop) = TI_breccia
+                enddo
+            else
+                do iloop = index_breccia + 1,index_bedrock
+                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
+                enddo
+            endif
+        enddo
+
+        !!! zone de transition
+        delta = depth_bedrock
+        do ig = 1,ngrid
+            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
+                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
+                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
+        enddo
+
+        do iloop = index_bedrock + 2,nsoil_PEM
+            do ig = 1,ngrid
+                inertiedat_PEM(ig,iloop) = TI_bedrock
+            enddo
+        enddo
+
+        write(*,*) 'PEMETAT0: TI done'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!b) Soil temperature
+        do islope = 1,nslope
+!            do ig = 1,ngrid
+!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
+!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
 !
-!      do iloop=index_bedrock+2,nsoil_PEM
-!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
-!      enddo
-!      enddo
-
-
-     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-
-
-   else
-! predictor corrector: restart from year 1 of the GCM and build the evolution of
-! tsoil at depth
-
-    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-    call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
-    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
-
-
-     do iloop = nsoil_GCM+1,nsoil_PEM
-       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
-     enddo
-   endif !found
-
-    do it = 1,timelen
-        do isoil = nsoil_GCM+1,nsoil_PEM
-        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
-        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)
-        enddo
-      enddo
-enddo ! islope
-
-write(*,*) 'PEMETAT0: SOIL TEMP done'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!3. Ice Table
-  call get_field("ice_table",ice_table,found)
-   if(.not.found) then
-      write(*,*)'PEM settings: failed loading <ice_table>'
-      write(*,*)'will reconstruct the values of the ice table given the current state'
-      ice_table(:,:) = -1  ! by default, no ice table  
-      ice_table_thickness(:,:) = -1  ! by default, no 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
-   endif
-
-write(*,*) 'PEMETAT0: ICE TABLE done'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!4. CO2 & H2O Adsorption
- if(adsorption_pem) then
-  do islope=1,nslope
-   write(num,fmt='(i2.2)') islope
-   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
-    if((.not.found)) then
-       m_co2_regolith_phys(:,:,:) = 0.
-       exit
-    endif
-
-  enddo
-
-  do islope=1,nslope
-   write(num,fmt='(i2.2)') islope
-   call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
-    if((.not.found2)) then
-       m_h2o_regolith_phys(:,:,:) = 0.
-      exit
-    endif
-
-  enddo
-
-    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,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)
-
-    if((.not.found)) then
-      deltam_co2_regolith_phys(:) = 0.
-    endif
-    if((.not.found2)) then
-       deltam_h2o_regolith_phys(:) = 0.
-    endif
-write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
-    endif
-endif ! soil_pem
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!. 5 water reservoir
-#ifndef CPP_STD
-   water_reservoir = 0.
-   call get_field("water_reservoir",water_reservoir,found)
-   if (.not. found) then
-       write(*,*)'Pemetat0: failed loading <water_reservoir>'
-       write(*,*)'will reconstruct the values from watercaptag'
-       where (watercaptag) water_reservoir = water_reservoir_nom
-   endif
-#endif
-
-call close_startphy
-
-else !No startfi, let's build all by hand
-
-    if (soil_pem) then
-
-!a) Thermal inertia
-   do islope = 1,nslope
-      do ig = 1,ngrid
-
-          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
-!!! transition
-             delta = depth_breccia
-             TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
-            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
-                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
-
-          do iloop=index_breccia+2,index_bedrock
-            TI_PEM(ig,iloop,islope) = TI_breccia
-         enddo
-         else ! we keep the high ti values
-           do iloop=index_breccia+1,index_bedrock
-                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
-           enddo
-         endif
-
-!! transition
-             delta = depth_bedrock
-             TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
-                      ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
-          do iloop=index_bedrock+2,nsoil_PEM
-            TI_PEM(ig,iloop,islope) = TI_bedrock
-         enddo
-      enddo
-enddo
-
- do iloop = 1,nsoil_GCM
-   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
- enddo
-!!! zone de transition
-delta = depth_breccia
-do ig = 1,ngrid
-      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
-inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
-            (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
-                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
-
-
- do iloop = index_breccia+2,index_bedrock
-
-       inertiedat_PEM(ig,iloop) = TI_breccia
-
- enddo
-else
-   do iloop = index_breccia+1,index_bedrock
-       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
-    enddo
-
-endif
-enddo
-
-!!! zone de transition
-delta = depth_bedrock
-do ig = 1,ngrid
-inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-            (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
-                      ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
-enddo
-
-
-
- do iloop = index_bedrock+2, nsoil_PEM
-   do ig = 1,ngrid
-      inertiedat_PEM(ig,iloop) = TI_bedrock
-   enddo
- enddo
-
-write(*,*) 'PEMETAT0: TI done'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!b) Soil temperature
-    do islope = 1,nslope
-!     do ig = 1,ngrid
-!        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
-!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
+!                do iloop=index_breccia+2,index_bedrock
+!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
+!                enddo
+!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
+!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
 !
-!       do iloop=index_breccia+2,index_bedrock
-!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
-!      enddo
-!        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
-!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
-
-!       do iloop=index_bedrock+2,nsoil_PEM
-!            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
-!      enddo
-
-!       enddo
-
-      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-  
-! First raw initialization 
-      do it = 1,timelen
-        do isoil = nsoil_GCM+1,nsoil_PEM
-          tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
-        enddo
-      enddo
-
-      do it = 1,timelen
-        do isoil = nsoil_GCM+1,nsoil_PEM
-          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
-        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)
-         enddo
-       enddo
-enddo !islope
-write(*,*) 'PEMETAT0: TSOIL done'
+!                do iloop=index_bedrock+2,nsoil_PEM
+!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
+!                enddo
+!            enddo
+            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+            call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+
+! First raw initialization
+            do it = 1,timelen
+                do isoil = nsoil_PCM + 1,nsoil_PEM
+                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
+                enddo
+            enddo
+
+            do it = 1,timelen
+                do isoil = nsoil_PCM + 1,nsoil_PEM
+                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
+                enddo
+            enddo
+
+            do isoil = nsoil_PCM + 1,nsoil_PEM
+                do ig = 1,ngrid
+                    watersoil_avg(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 !islope
+        write(*,*) 'PEMETAT0: TSOIL done'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table
-       ice_table(:,:) = -1.  ! by default, no ice table
-       ice_table_thickness(:,:) = -1.
-       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'
+        ice_table = -1.  ! by default, no ice table
+        ice_table_thickness = -1.
+        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
+        call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
+        do islope = 1,nslope
+            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+        enddo
+        write(*,*) 'PEMETAT0: Ice table done'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed
-if (adsorption_pem) then
-    m_co2_regolith_phys(:,:,:) = 0.
-    m_h2o_regolith_phys(:,:,:) = 0.
-
-    call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,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.
-    deltam_h2o_regolith_phys(:) = 0.
-endif
-
-write(*,*) 'PEMETAT0: CO2 adsorption done'
-endif !soil_pem
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!. e) water reservoir
-#ifndef CPP_STD
-    water_reservoir = 0.
-    where (watercaptag) water_reservoir = water_reservoir_nom
-#endif
-
+        if (adsorption_pem) then
+            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, &
+                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
+            deltam_co2_regolith_phys = 0.
+            deltam_h2o_regolith_phys = 0.
+        endif
+
+        write(*,*) 'PEMETAT0: CO2 adsorption done'
+    endif !soil_pem
 endif ! of if (startphy_file)
 
@@ -495,3 +475,2 @@
 
 END MODULE pemetat0_mod
-
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 3149)
@@ -1,9 +1,10 @@
-module pemredem
+MODULE pemredem
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
-!!! Purpose: Write specific netcdf restart for the PEM 
-!!! 
-!!! 
-!!! Author: LL, inspired by phyredem  from the GCM
+!!! Purpose: Write specific netcdf restart for the PEM
+!!!
+!!!
+!!! Author: LL, inspired by phyredem from the PCM
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -11,128 +12,116 @@
 implicit none
 
+!=======================================================================
 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 time_phylmdz_mod,  only: daysec
+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 comcstfi_mod, only: g, mugaz, omeg, rad, rcp
+    use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
+    use comcstfi_mod, only: g, mugaz, omeg, rad, rcp
 #endif
 
-  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
-  
-  ! Create physics start file
-  call open_restartphy(filename)
-  
-  ! Write the mid-layer depths
-  call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)
-  
-  ! Write longitudes
-  call put_field("longitude","Longitudes of physics grid",lonfi)
-  
-  ! Write latitudes
-  call put_field("latitude","Latitudes of physics grid",latfi)
-  
-  ! Write mesh areas
-  call put_field("area","Mesh area",cell_area)
-     
-  ! Multidimensionnal variables (nogcm undermesh slope statistics)
-  call put_var("def_slope","slope criterium stages",def_slope)
-  call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
+implicit none
 
+character(*),                  intent(in) :: filename
+integer,                       intent(in) :: nsoil_PEM, ngrid, nslope
+real, dimension(ngrid),        intent(in) :: lonfi, latfi
+real,                          intent(in) :: day_ini, time
+real, dimension(ngrid),        intent(in) :: cell_area     ! boundaries for bining of the slopes
+real, dimension(ngrid + 1),    intent(in) :: def_slope     ! boundaries for bining of the slopes
+real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! undermesh statistics
 
-  ! Close file
-  call close_restartphy
-  
-end subroutine pemdem0
+! Create physics start file
+call open_restartphy(filename)
 
-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)
+! Write the mid-layer depths
+call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)
 
+! Write longitudes
+call put_field("longitude","Longitudes of physics grid",lonfi)
 
-  ! 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 time_evol_mod, only: year_bp_ini, convert_years
-                
-  implicit none
-  
+! Write latitudes
+call put_field("latitude","Latitudes of physics grid",latfi)
+
+! Write mesh areas
+call put_field("area","Mesh area",cell_area)
+
+! Multidimensionnal variables (nopcm undermesh slope statistics)
+call put_var("def_slope","slope criterium stages",def_slope)
+call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
+
+! Close file
+call close_restartphy
+
+END SUBROUTINE pemdem0
+
+!=======================================================================
+
+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,h2o_ice)
+
+! 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 time_evol_mod, only: year_bp_ini, convert_years
+
+implicit none
+
 #ifndef CPP_STD
-  include "callkeys.h"
+    include "callkeys.h"
 #endif
-  
-  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)
 
-  integer           :: islope
-  character*2       :: num
-  integer           :: iq
-  character(len=30) :: txt  ! To store some text
-  real              :: Year ! Year of the simulation
+character(*),                            intent(in) :: filename
+integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope, 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
+real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth       ! under mesh bining according to slope
+real, dimension(ngrid,nslope),           intent(in) :: ice_table_thickness   ! under mesh bining according to slope
+real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
+real, dimension(ngrid,nslope),           intent(in) :: h2o_ice
 
-  ! Open file
-  call open_restartphy(filename)
-  
-  ! First variable to write must be "Time", in order to correctly
-  ! set time counter in file
-  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)
+integer       :: islope
+character(2)  :: num
+integer       :: iq
+character(30) :: txt  ! To store some text
+real          :: Year ! Year of the simulation
 
-    if(soil_pem) then
-  
+! Open file
+call open_restartphy(filename)
+
+! First variable to write must be "Time", in order to correctly
+! set time counter in file
+Year = (year_bp_ini + i_myear)*convert_years
+call put_var("Time","Year of simulation",Year)
+call put_field('h2o_ice','h2o_ice',h2o_ice,Year)
+
+if (soil_pem) then
   ! Multidimensionnal variables (undermesh slope statistics)
-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)
-  call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
-                 inertiesoil_slope_PEM(:,:,islope),Year)
+    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)
+        call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type",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)
+        call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith",m_h2o_regolith(:,:,islope),Year)
+    enddo
+    call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
+    call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
+    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
+endif ! soil_pem
 
-  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)
-enddo
+! Close file
+call close_restartphy
 
-  call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
+END SUBROUTINE pemdem1
 
-  call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
-
-  call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
-
-   endif !soil_pem
-
-  ! Close file
-  call close_restartphy
-  
-end subroutine pemdem1
-
-end module pemredem
+END MODULE pemredem
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3149)
@@ -15,7 +15,7 @@
 !=======================================================================
 
-SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,          &
-                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
-                         watersurf_density_ave,watersoil_density)
+SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM_phys,ps_timeseries,          &
+                         min_co2_ice,min_h2o_ice,tsurf_avg,tsoil_avg,tsurf_PCM,tsoil_PCM,q_co2,q_h2o,co2_ice_slope, &
+                         watersurf_density_avg,watersoil_density)
 use comsoil_h,             only: nsoilmx
 use comsoil_h_PEM,         only: soil_pem
@@ -36,21 +36,21 @@
 ! Arguments:
 
-character(len = *), intent(in) :: filename                            ! File name
-integer,            intent(in) :: timelen                             ! number of times stored in the file
-integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
+character(*), intent(in) :: filename                            ! File name
+integer,      intent(in) :: timelen                             ! number of times stored in the file
+integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
 ! Ouputs
 real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
 real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
-real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_PCM_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
 real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
 real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
 real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
-real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
+real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per slope of the year [kg/m^2]
 !SOIL
-real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
-real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
-real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
+real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg             ! Average surface temperature of the concatenated file [K]
+real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg             ! Average soil temperature of the concatenated file [K]
+real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_PCM             ! Surface temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_PCM             ! Soil temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg ! Water density at the surface [kg/m^3]
 real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
 !=======================================================================
@@ -67,17 +67,17 @@
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm               ! CO2 volume mixing ratio in the first layer [mol/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                    ! Surface Pressure [Pa]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_PCM               ! CO2 volume mixing ratio in the first layer [mol/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_PCM                    ! Surface Pressure [Pa]
 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn             ! Average surface temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn             ! Average soil temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn             ! Surface temperature of the concatenated file, time series [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn             ! Soil temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Average surface temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Average soil temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_PCM_dyn             ! Surface temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_PCM_dyn             ! Soil temperature of the concatenated file, time series [K]
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_ave ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
 real, dimension(ngrid,nslope,timelen)                               :: watersurf_density         ! Water density at the surface, physical grid, time series [kg/m^3]
@@ -88,5 +88,5 @@
 
 ! Dowload the data from the file
-call get_var3("ps",ps_GCM)
+call get_var3("ps",ps_PCM)
 write(*,*) "Data for surface pressure downloaded!"
 
@@ -104,5 +104,5 @@
     write(*,*) "Data for h2o_ice_s downloaded!"
 
-    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
+    call get_var3("tsurf",tsurf_PCM_dyn(:,:,1,:))
     write(*,*) "Data for tsurf downloaded!"
 
@@ -115,5 +115,5 @@
 
     if (soil_pem) then
-        call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
+        call get_var4("soiltemp",tsoil_PCM_dyn(:,:,:,1,:))
         write(*,*) "Data for soiltemp downloaded!"
 
@@ -140,5 +140,5 @@
     do islope = 1,nslope
         write(num,'(i2.2)') islope
-        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
+        call get_var3("tsurf_slope"//num,tsurf_PCM_dyn(:,:,islope,:))
     enddo
     write(*,*) "Data for tsurf downloaded!"
@@ -160,5 +160,5 @@
         do islope = 1,nslope
             write(num,'(i2.2)') islope
-            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
+            call get_var4("soiltemp_slope"//num,tsoil_PCM_dyn(:,:,:,islope,:))
         enddo
         write(*,*) "Data for soiltemp downloaded!"
@@ -182,18 +182,18 @@
 write(*,*) "Computing the min of h2o_ice..."
 where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
-min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap,4)
+min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
 write(*,*) "Computing the min of co2_ice..."
 where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
-min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn + perennial_co2ice,4)
+min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
 
 ! Compute averages over the year for each point
 write(*,*) "Computing the average of tsurf..."
-tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
+tsurf_avg_dyn = sum(tsurf_PCM_dyn,4)/timelen
 
 if (soil_pem) then
     write(*,*) "Computing average of tsoil..."
-    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
+    tsoil_avg_dyn = sum(tsoil_PCM_dyn,5)/timelen
     write(*,*) "Computing average of waterdensity_surface..."
-    watersurf_density_dyn_ave(:,:,:) = sum(watersurf_density_dyn(:,:,:,:),4)/timelen
+    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
 endif
 
@@ -215,5 +215,5 @@
             endif
             mmean = 1/(A*q_co2_dyn(i,j,t) + B)
-            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
+            vmr_co2_PCM(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
         enddo
     enddo
@@ -221,6 +221,6 @@
 
 #ifndef CPP_1D
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_PCM,vmr_co2_PCM_phys)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_PCM,ps_timeseries)
     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
@@ -231,22 +231,22 @@
         if (soil_pem) then
             do l = 1,nsoilmx
-                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
+                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
                 do t = 1,timelen
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_PCM_dyn(:,:,l,islope,t),tsoil_PCM(:,l,islope,t))
                     call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
                 enddo
             enddo
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_ave(:,:,islope),watersurf_density_ave(:,islope))
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
         endif ! soil_pem
         do t = 1,timelen
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_PCM_dyn(:,:,islope,t),tsurf_PCM(:,islope,t))
             call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
         enddo
     enddo
 
-    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
+    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
 #else
-    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
-    ps_timeseries(1,:) = ps_GCM(1,1,:)
+    vmr_co2_PCM_phys(1,:) = vmr_co2_PCM(1,1,:)
+    ps_timeseries(1,:) = ps_PCM(1,1,:)
     q_co2(1,:) = q_co2_dyn(1,1,:)
     q_h2o(1,:) = q_h2o_dyn(1,1,:)
@@ -254,12 +254,12 @@
     min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
     if (soil_pem) then
-        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
-        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
+        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
+        tsoil_PCM(1,:,:,:) = tsoil_PCM_dyn(1,1,:,:,:)
         watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
-        watersurf_density_ave(1,:) = watersurf_density_dyn_ave(1,1,:)
+        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
     endif ! soil_pem
-    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
+    tsurf_PCM(1,:,:) = tsurf_PCM_dyn(1,1,:,:)
     co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
-    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
+    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
 #endif
 
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3149)
@@ -8,7 +8,7 @@
 
 SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
-                                 vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new)
+                                 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
 
-use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
+use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
 
 implicit none
@@ -24,11 +24,11 @@
 !   INPUT
 integer,                        intent(in) :: timelen, ngrid, nslope
-real, dimension(ngrid,timelen), intent(in) :: vmr_co2_gcm                 ! physical point field : Volume mixing ratio of co2 in the first layer
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM                 ! physical point field : Volume mixing ratio of co2 in the first layer
 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_pem                 ! physical point field : Volume mixing ratio of co2 in the first layer
-real, dimension(ngrid,timelen), intent(in) :: ps_GCM_2                    ! physical point field : Surface pressure in the GCM
-real,                           intent(in) :: global_ave_press_GCM        ! global averaged pressure at previous timestep
-real,                           intent(in) :: global_ave_press_new        ! global averaged pressure at current timestep
+real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2                    ! physical point field : Surface pressure in the PCM
+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) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope(kg/m^2)
+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
@@ -40,4 +40,6 @@
 real    :: coef, ave
 
+write(*,*) "Update of the CO2 tendency from the current pressure"
+    
 ! Evolution of the water ice for each physical point
 do i = 1,ngrid
@@ -47,6 +49,6 @@
         if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
             do t=1,timelen
-                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 &
-                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
+                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
+                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**4
             enddo
         endif
Index: trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3149)
@@ -7,5 +7,5 @@
 !=======================================================================
 
-SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_GCM,TI_PEM)
+SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
 
 !      use netcdf
@@ -39,6 +39,6 @@
 integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
 integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
-integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
-real, dimension(ngrid,nsoil_GCM,nslope), intent(in) :: TI_GCM    ! Thermal inertia  in the GCM [SI]
+integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the PCM
+real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM    ! Thermal inertia  in the PCM [SI]
 
 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
@@ -57,10 +57,10 @@
 lay1 = 2.e-4
 alpha = 2
-do iloop = 0,nsoil_GCM - 1
+do iloop = 0,nsoil_PCM - 1
     mlayer_PEM(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
 enddo
 
 do iloop = 0, nsoil_PEM-1
-    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_GCM-1)) then
+    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM-1)) then
         index_powerlaw = iloop
         exit
@@ -68,6 +68,6 @@
 enddo
 
-do iloop = nsoil_GCM, nsoil_PEM-1 
-    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_GCM)-0.5))
+do iloop = nsoil_PCM, nsoil_PEM-1 
+    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_PCM)-0.5))
 enddo
 
@@ -83,10 +83,10 @@
 do ig = 1,ngrid
     do islope = 1,nslope
-        do iloop = 1,nsoil_GCM
-            TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
+        do iloop = 1,nsoil_PCM
+            TI_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
         enddo
-        if (nsoil_PEM > nsoil_GCM) then
-            do iloop = nsoil_GCM + 1,nsoil_PEM
-                TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
+        if (nsoil_PEM > nsoil_PCM) then
+            do iloop = nsoil_PCM + 1,nsoil_PEM
+                TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
             enddo
         endif
Index: trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90	(revision 3149)
@@ -89,8 +89,7 @@
  REAL :: ice_bottom_depth
 
+write(*,*) "Update soil propreties"
+        
 ! 1. Modification of the regolith thermal inertia.
-
-
-
  do islope = 1,nslope
    regolith_inertia(:,islope) = inertiedat_PEM(:,1)
@@ -112,5 +111,4 @@
 
 ! 2. Build new Thermal Inertia 
-
 do  islope=1,nslope
    do ig = 1,ngrid
Index: trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90	(revision 3149)
+++ trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90	(revision 3149)
@@ -0,0 +1,151 @@
+MODULE stopping_crit_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Criterions to check if the PEM needs to call the PCM
+!!! Author: RV & LL, 02/2023
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)
+
+use time_evol_mod, only: water_ice_criterion
+use comslope_mod,  only: subslope_dist, nslope
+
+implicit none
+
+!=======================================================================
+!
+! Routine to check if the h2o ice criterion to stop the PEM is reached
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,                       intent(in) :: ngrid     ! # of physical grid points
+real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
+real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of water ice
+real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
+!   OUTPUT
+integer, intent(inout) :: stopPEM ! Stopping criterion code
+!   local:
+!   ------
+integer :: i, islope ! Loop
+real    :: surf_now  ! Current surface of h2o ice
+
+!=======================================================================
+! Computation of the present surface of h2o ice
+surf_now = 0.
+do i = 1,ngrid
+    do islope = 1,nslope
+        if (h2o_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
+    enddo
+enddo
+
+! Check of the criterion
+if (surf_now < surf_ini*(1. - water_ice_criterion)) then
+    stopPEM = 1
+    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
+    write(*,*) "surf_now < surf_ini*(1. - water_ice_criterion)", surf_now < surf_ini*(1. - water_ice_criterion)
+    write(*,*) "Current surface of water ice sublimating =", surf_now
+    write(*,*) "Initial surface of water ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
+else if (surf_now > surf_ini*(1. + water_ice_criterion)) then
+    stopPEM = 1
+    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
+    write(*,*) "surf_now > surf_ini*(1. + water_ice_criterion)", surf_now > surf_ini*(1. + water_ice_criterion)
+    write(*,*) "Current surface of water ice sublimating =", surf_now
+    write(*,*) "Initial surface of water ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
+endif
+
+if (abs(surf_ini) < 1.e-5) stopPEM = 0
+
+END SUBROUTINE stopping_crit_h2o_ice
+
+!=======================================================================
+
+SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
+
+use time_evol_mod, only: co2_ice_criterion, ps_criterion
+use comslope_mod,  only: subslope_dist
+
+implicit none
+
+!=======================================================================
+!
+! Routine to check if the co2 and pressure criteria to stop the PEM are reached
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
+real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
+real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of water ice
+real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
+real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
+real,                          intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
+!   OUTPUT
+integer, intent(inout) :: stopPEM ! Stopping criterion code
+
+!   local:
+!   ------
+integer :: i, islope ! Loop
+real    :: surf_now  ! Current surface of co2 ice
+
+!=======================================================================
+! Computation of the present surface of co2 ice
+surf_now = 0.
+do i = 1,ngrid
+    do islope = 1,nslope
+        if (co2_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
+    enddo
+enddo
+
+! Check of the criterion
+if (surf_now < surf_ini*(1. - co2_ice_criterion)) then
+    stopPEM = 3
+    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
+    write(*,*) "surf_now < surf_ini*(1. - co2_ice_criterion)", surf_now < surf_ini*(1. - co2_ice_criterion)
+    write(*,*) "Current surface of co2 ice sublimating =", surf_now
+    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
+else if (surf_now > surf_ini*(1. + co2_ice_criterion)) then
+    stopPEM = 3
+    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
+    write(*,*) "surf_now > surf_ini*(1. + co2_ice_criterion)", surf_now > surf_ini*(1. + co2_ice_criterion)
+    write(*,*) "Current surface of co2 ice sublimating =", surf_now
+    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
+    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
+endif
+
+if (abs(surf_ini) < 1.e-5) stopPEM = .false.
+
+if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
+    stopPEM = 4
+    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
+    write(*,*) "global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)", global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)
+    write(*,*) "Current global pressure =", global_avg_press_new
+    write(*,*) "PCM global pressure =", global_avg_press_PCM
+    write(*,*) "Percentage of change accepted =", ps_criterion*100.
+else if (global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)) then
+    stopPEM = 4
+    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
+    write(*,*) "global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)", global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)
+    write(*,*) "Current global pressure =", global_avg_press_new
+    write(*,*) "PCM global pressure =", global_avg_press_PCM
+    write(*,*) "Percentage of change accepted =", ps_criterion*100.
+endif
+
+END SUBROUTINE stopping_crit_co2
+
+END MODULE
Index: trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3145)
+++ trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3149)
@@ -12,5 +12,5 @@
 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_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
 logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
 
