Index: trunk/LMDZ.COMMON/libf/evolution/abort_pem.F
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/abort_pem.F	(revision 2888)
+++ trunk/LMDZ.COMMON/libf/evolution/abort_pem.F	(revision 2888)
@@ -0,0 +1,59 @@
+!
+! $Id: abort_pem.F 1425 2010-09-02 13:45:23Z lguez $
+!
+c
+c
+      SUBROUTINE abort_pem(modname, message, ierr)
+     
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin_dump
+      USE ioipsl_getincom
+#endif
+
+#ifdef CPP_XIOS
+    ! ug Pour les sorties XIOS
+      USE wxios
+#endif
+
+#include "iniprint.h"
+ 
+C
+C Stops the simulation cleanly, closing files and printing various
+C comments
+C
+C  Input: modname = name of calling program
+C         message = stuff to print
+C         ierr    = severity of situation ( = 0 normal )
+
+      character(len=*), intent(in):: modname
+      integer, intent(in):: ierr
+      character(len=*), intent(in):: message
+
+      write(lunout,*) 'in abort_gcm'
+
+#ifdef CPP_XIOS
+    !Fermeture propre de XIOS
+      CALL wxios_close()
+#endif
+
+#ifdef CPP_IOIPSL
+      call histclo
+      call restclo
+#endif
+      call getin_dump
+c     call histclo(2)
+c     call histclo(3)
+c     call histclo(4)
+c     call histclo(5)
+      write(lunout,*) 'Stopping in ', modname
+      write(lunout,*) 'Reason = ',message
+      if (ierr .eq. 0) then
+        write(lunout,*) 'Everything is cool'
+        stop
+      else
+        write(lunout,*) 'Houston, we have a problem ', ierr
+        stop 1
+      endif
+      END
Index: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2888)
@@ -1,4 +1,7 @@
   module adsorption_mod  
   implicit none
+  LOGICAL adsorption_pem ! True by default, to compute adsorption/desorption. Read in  pem.def
+  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
+  real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
   contains
 
@@ -10,4 +13,27 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+  subroutine ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
+
+  implicit none
+  integer,intent(in) :: ngrid ! number of atmospheric columns
+  integer,intent(in) :: nslope ! number of slope within a mesh 
+  integer,intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
+    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+    allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+  end subroutine ini_adsorption_h_PEM
+
+  subroutine end_adsorption_h_PEM
+
+  implicit none
+    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
+    if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
+  end subroutine end_adsorption_h_PEM
+
+
+
+
+
 
   subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
Index: trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2888)
@@ -31,5 +31,9 @@
 
          enddo
-
+         if (ig.eq.400) then 
+         write(*,*) 'ig,islope,diff',ig,islope,diff_rho(:)
+         write(*,*) 'rho surf',rhowatersurf_ave(ig,islope)
+         write(*,*) 'rho soil',rhowatersoil_ave(ig,:,islope)
+         endif
          if(diff_rho(1) > 0) then                                   ! ice is at the surface
            ice_table(ig,islope) = 0.
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2888)
@@ -15,13 +15,11 @@
   real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients [SI]
   real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients [SI]
-  real,save,allocatable :: alph_PEM(:,:,:)      ! (FC) alpha_k coefficients [SI]
-  real,save,allocatable :: beta_PEM(:,:,:)      ! beta_k coefficients [SI]
+  real,save,allocatable :: alph_PEM(:,:)      ! (FC) alpha_k coefficients [SI]
+  real,save,allocatable :: beta_PEM(:,:)      ! beta_k coefficients [SI]
   real,save :: mu_PEM                           ! mu coefficient [SI]
-  real,parameter :: fluxgeo = 30e-3             ! Geothermal flux [W/m^2]
-  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
-  real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
+  real,save :: fluxgeo                               ! Geothermal flux [W/m^2]
   LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
   real,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
-
+  real,save :: water_reservoir_nom
 contains
 
@@ -40,9 +38,8 @@
     allocate(coefq_PEM(0:nsoilmx_PEM-1))
     allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
-    allocate(alph_PEM(ngrid,nsoilmx_PEM-1,nslope))
-    allocate(beta_PEM(ngrid,nsoilmx_PEM-1,nslope))
+    allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
+    allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
     allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 
-    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
-    allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+    allocate(water_reservoir(ngrid))
     allocate(water_reservoir(ngrid))
   end subroutine ini_comsoil_h_PEM
@@ -64,6 +61,5 @@
     if (allocated(beta_PEM)) deallocate(beta_PEM)
     if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
-    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
-    if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
+    if (allocated(water_reservoir)) deallocate(water_reservoir)
     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 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 2888)
@@ -14,9 +14,9 @@
 #endif
   
-  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, alpha_criterion, &
+  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, ice_criterion, ps_criterion, &
                 Max_iter_pem, evol_orbit_pem
-  USE comsoil_h_pem, only: soil_pem
-
-  
+  USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom
+  USE adsorption_mod,only: adsorption_pem
+  CHARACTER(len=20),parameter :: modname ='conf_pem' 
 !PEM parameter
   year_bp_ini=0.
@@ -26,6 +26,9 @@
   CALL getin('dt_pem', dt_pem)
 
-  alpha_criterion=0.2
-  CALL getin('alpha_criterion', alpha_criterion)
+  ice_criterion=0.2
+  CALL getin('ice_criterion', ice_criterion)
+
+  ps_criterion = 0.15
+  CALL getin('ps_criterion',ps_criterion)
 
   evol_orbit_pem=.false.
@@ -38,4 +41,25 @@
   CALL getin('soil_pem', soil_pem)
 
+  adsorption_pem = .true.
+  CALL getin('adsorption_pem',adsorption_pem)
+
+  fluxgeo = 0.
+  CALL getin('Fluxgeo_PEM',fluxgeo)
+  print*,'Flux Geothermal is set to',fluxgeo
+
+  if ((not(soil_pem)).and.adsorption_pem) then
+       print*,'Adsorption must be used when soil_pem = T'
+       call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
+        endif
+  
+
+  if ((not(soil_pem)).and.fluxgeo.gt.0.) then
+       print*,'Soil is not activated but Flux Geo > 0.'
+       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.'",1)
+        endif
+  
+
+   water_reservoir_nom = 1e4
+  CALL getin('water_reservoir_nom',water_reservoir_nom)
   END SUBROUTINE conf_pem
 
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90	(revision 2886)
+++ 	(revision )
@@ -1,67 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE criterion_co2_ice_stop(cell_area,initial_co2_ice,co2ice,STOPPING,ngrid,latitude,n_band_lat)
-
-  USE temps_mod_evol, ONLY: alpha_criterion
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!  Routine that checks if the criterion to stop the PEM is reached
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-  INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
-  REAL,    intent(in) ::  co2ice(ngrid)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) ::  latitude(ngrid)          ! physical point field : Latitude
-  REAL,    intent(in) ::  initial_co2_ice(n_band_lat)  ! Initial/Actual surface of water ice
-
-!   OUTPUT
-  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
-
-!   local:
-!   -----
-  INTEGER :: i,j,n_band_lat                    ! Loop
-  REAL :: present_co2(n_band_lat)  ! Initial/Actual surface of water ice
-  REAL :: pi
-
-!=======================================================================
-
-      pi=4.D0*DATAN(1.D0)
-
-!   initialisation to false
-    STOPPING=.FALSE.
-
-     do j=1,n_band_lat
-        present_co2(j)=0.
-     enddo
-
-  do i=1,ngrid
-            j=floor((latitude(i)+(pi/2))/(pi)*n_band_lat)+1
-      if(j.GT.n_band_lat) then
-          j=n_band_lat
-      endif
-      present_co2(j)=present_co2(j)+co2ice(i)*cell_area(i)
-  enddo
-  
-!   check of the criterion
-  do j=1,n_band_lat
-    if(present_co2(j).LT.initial_co2_ice(j)*(1-alpha_criterion) .OR. &
-       present_co2(j).GT.initial_co2_ice(j)*(1+alpha_criterion)) then
-         STOPPING=.TRUE. 
-    endif
-  enddo
-
-END SUBROUTINE criterion_co2_ice_stop
-
-
-
-
-
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90	(revision 2886)
+++ 	(revision )
@@ -1,80 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE criterion_ice_stop_slope(cell_area,ini_surf,qsurf,STOPPING,STOPPING_ps,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
-
-  USE temps_mod_evol, ONLY: alpha_criterion
-  use comslope_mod, ONLY: subslope_dist
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!  Routine that checks if the criterion to stop the PEM is reached
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
-  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf
-  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
-  REAL,    intent(in) :: global_ave_press_GCM
-  REAL,    intent(in) :: global_ave_press_new
-
-!   OUTPUT
-  LOGICAL, intent(out) :: STOPPING,STOPPING_ps              ! Logical : is the criterion reached?
-
-!   local:
-!   -----
-  INTEGER :: i,islope   ! Loop
-  REAL :: present_surf  ! Initial/Actual surface of water ice
-
-!=======================================================================
-
-!   initialisation to false
-    STOPPING=.FALSE.
-    STOPPING_ps=.FALSE.
-
-!   computation of the actual surface
-  present_surf=0.
-  do i=1,ngrid
-   do islope=1,nslope
-      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
-      endif
-   enddo
-  enddo
-  
-!   check of the criterion
-  if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+alpha_criterion)) then
-    STOPPING=.TRUE. 
-    print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
-    print *, "Current surface of co2 ice sublimating=", present_surf
-    print *, "Initial surface of co2 ice sublimating=", ini_surf
-    print *, "Percentage of change accepted=", alpha_criterion*100
-    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
-  endif
-
-  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
-       STOPPING=.FALSE.
-  endif
-
-  if(global_ave_press_new.LT.global_ave_press_GCM*(0.9) .OR. &
-     global_ave_press_new.GT.global_ave_press_GCM*(1.1)) then
-    STOPPING_ps=.TRUE. 
-    print *, "Reason of stopping : The global pressure reach the threshold:"
-    print *, "Current global pressure=", global_ave_press_new
-    print *, "GCM global pressure=", global_ave_press_GCM
-    print *, "Percentage of change accepted=", 0.1*100
-    print *, "global_ave_press_new<global_ave_press_GCM*(0.9)", (global_ave_press_new.LT.global_ave_press_GCM*(0.9))
-  endif
-
-END SUBROUTINE criterion_ice_stop_slope
-
-
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90	(revision 2886)
+++ 	(revision )
@@ -1,69 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE criterion_ice_stop_water_slope(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
-
-  USE temps_mod_evol, ONLY: alpha_criterion
-	  use comslope_mod, ONLY: subslope_dist,nslope
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!  Routine that checks if the criterion to stop the PEM is reached
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-  INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
-  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf
-  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
-
-!   OUTPUT
-  LOGICAL, intent(out) :: STOPPING              ! Logical : 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 water ice sublimating
-  present_surf=0.
-  do i=1,ngrid
-    do islope=1, nslope
-      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
-      endif
-    enddo
-  enddo
-  
-!   check of the criterion
-  if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+alpha_criterion)) then
-    STOPPING=.TRUE. 
-    print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
-    print *, "Current surface of water ice sublimating=", present_surf
-    print *, "Initial surface of water ice sublimating=", ini_surf
-    print *, "Percentage of change accepted=", alpha_criterion*100
-    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
-  endif
-
-  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
-    STOPPING=.FALSE.
-  endif
-END SUBROUTINE criterion_ice_stop_water_slope
-
-
-
-
-
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 2888)
+++ trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 2888)
@@ -0,0 +1,162 @@
+  module criterion_pem_stop_mod
+  implicit none
+
+  contains
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Criterions to check if the PEM needs to call the GCM !!!
+!!! Author: RV & LL, 02/2023
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
+
+  USE temps_mod_evol, ONLY: 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 grid physical grid points 
+  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
+  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
+  REAL,    intent(in) :: ini_surf
+  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
+
+!   OUTPUT
+  LOGICAL, intent(out) :: STOPPING              ! Logical : 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 water ice sublimating
+  present_surf=0.
+  do i=1,ngrid
+    do islope=1, nslope
+      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
+         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
+      endif
+    enddo
+  enddo
+  
+!   check of the criterion
+  if(present_surf.LT.ini_surf*(1-ice_criterion) .OR. &
+     present_surf.GT.ini_surf*(1+ice_criterion)) then
+    STOPPING=.TRUE. 
+    print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
+    print *, "Current surface of water ice sublimating=", present_surf
+    print *, "Initial surface of water ice sublimating=", ini_surf
+    print *, "Percentage of change accepted=", ice_criterion*100
+    print *, "present_surf<ini_surf*(1-ice_criterion)", (present_surf.LT.ini_surf*(1-ice_criterion))
+  endif
+
+  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
+    STOPPING=.FALSE.
+  endif
+END SUBROUTINE criterion_waterice_stop
+
+! ------------------------------------------------------------------------------------------------
+
+
+SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
+
+  USE temps_mod_evol, ONLY: ice_criterion,ps_criterion
+  use comslope_mod, ONLY: subslope_dist
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+!  Routine that checks if the criterion to stop the PEM is reached
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+
+!   INPUT
+  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
+  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
+  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
+  REAL,    intent(in) :: ini_surf
+  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
+  REAL,    intent(in) :: global_ave_press_GCM
+  REAL,    intent(in) :: global_ave_press_new
+
+!   OUTPUT
+  LOGICAL, intent(out) :: STOPPING_ice              ! Logical : is the criterion for ice reached?
+  LOGICAL, intent(out) :: STOPPING_ps               ! Logical : 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 actual surface
+  present_surf=0.
+  do i=1,ngrid
+   do islope=1,nslope
+      if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
+         present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
+      endif
+   enddo
+  enddo
+  
+!   check of the criterion
+  if(present_surf.LT.ini_surf*(1-ice_criterion) .OR. &
+     present_surf.GT.ini_surf*(1+ice_criterion)) then
+    STOPPING_ice=.TRUE. 
+    print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
+    print *, "Current surface of co2 ice sublimating=", present_surf
+    print *, "Initial surface of co2 ice sublimating=", ini_surf
+    print *, "Percentage of change accepted=", ice_criterion*100
+    print *, "present_surf<ini_surf*(1-ice_criterion)", (present_surf.LT.ini_surf*(1-ice_criterion))
+  endif
+
+  if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
+       STOPPING_ice=.FALSE.
+  endif
+
+  if(global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion) .OR. &
+     global_ave_press_new.GT.global_ave_press_GCM*(1+ps_criterion)) then
+    STOPPING_ps=.TRUE. 
+    print *, "Reason of stopping : The global pressure reach the threshold:"
+    print *, "Current global pressure=", global_ave_press_new
+    print *, "GCM global pressure=", global_ave_press_GCM
+    print *, "Percentage of change accepted=", ps_criterion*100
+    print *, "global_ave_press_new<global_ave_press_GCM*(ps_criterion)", (global_ave_press_new.LT.global_ave_press_GCM*(1-ps_criterion))
+  endif
+
+END SUBROUTINE criterion_co2_stop
+
+
+END MODULE
+
+
+
+
+
+
Index: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2888)
@@ -7,5 +7,5 @@
   USE temps_mod_evol, ONLY: dt_pem
   use comslope_mod, ONLY: subslope_dist
-
+  use criterion_pem_stop_mod, only: criterion_waterice_stop
       IMPLICIT NONE
 
@@ -80,5 +80,5 @@
     print *, "Tendencies on ice increasing=", pos_tend
     print *, "This can be due to the absence of water ice in the PCM run!!"
-      call criterion_ice_stop_water_slope(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
+      call criterion_waterice_stop(cell_area,1.,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
       do i=1,ngrid
          do islope=1,nslope
Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 2888)
+++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 2888)
@@ -0,0 +1,76 @@
+module ice_table_mod
+
+  implicit none
+
+  contains
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the ice table in two ways: dynamic and at equilibrium
+!!! Author:  LL, 02/2023
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
+   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,ice_table)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the ice table depth knowing the yearly average water
+!!! density at the surface and at depth.
+!!! Computations are made following the methods in Schorgofer et al., 2005
+!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
+!!! 
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+#ifndef CPP_STD
+    USE comsoil_h_PEM, only: layer_PEM                             ! Depth of the vertical grid 
+    implicit none 
+
+    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM 
+    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
+    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
+    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
+    real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
+    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root 
+    integer ig, islope,isoil                                       ! loop variables
+    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
+
+
+     do ig = 1,ngrid
+      if(watercaptag(ig)) then
+         ice_table(ig,:) = 0.
+      else
+        do islope = 1,nslope
+           ice_table(ig,islope) = -10.
+
+         do isoil = 1,nsoil_PEM
+           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
+
+         enddo
+         if(diff_rho(1) > 0) then                                   ! ice is at the surface
+           ice_table(ig,islope) = 0.
+         else
+           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
+             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
+               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
+               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
+               ice_table(ig,islope) = -z2/z1
+               exit
+             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
+           enddo
+          endif ! diff_rho(1) > 0
+        enddo
+       endif ! watercaptag
+      enddo
+!=======================================================================
+      RETURN
+#endif
+      END
+
+
+
+end module
Index: trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90	(revision 2888)
@@ -1,14 +1,19 @@
 SUBROUTINE info_run_PEM(year_iter, criterion_stop)
 
-      IMPLICIT NONE
+   IMPLICIT NONE
 
   INTEGER, intent(in) :: year_iter, criterion_stop                  ! # of year and reason to stop
 
- open(1,file='info_run_PEM.txt', status='new')
-
- write(1,*)  year_iter, criterion_stop
-
- close(1)
+  logical :: exist
 
 
+  inquire(file='info_run_PEM.txt', exist=exist)
+  if (exist) then
+    open(12, file='info_run_PEM.txt', status="old", position="append", action="write")
+  else
+    open(12, file='info_run_PEM.txt', status="new", action="write")
+  end if
+  write(12, *) year_iter, criterion_stop
+  close(12)
+
 END SUBROUTINE info_run_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2888)
@@ -84,13 +84,15 @@
       USE geometry_mod, only: latitude_deg 
       use conf_pem_mod, only: conf_pem
-      use pemredem, only:  pemdem1
+      use pemredem, only:  pemdem0,pemdem1
       use co2glaciers_mod,only: co2glaciers_evol
+      use criterion_pem_stop_mod,only: criterion_waterice_stop,criterion_co2_stop
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
       use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
-                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia          
+                              TI_PEM,inertiedat_PEM,         & ! soil thermal inertia          
                               tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
-                              fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys,  & ! geothermal flux, mass of co2 & h2o in the regolith
-                              water_reservoir
-      use adsorption_mod, only : regolith_adsorption
+                              fluxgeo,water_reservoir                             ! geothermal flux
+      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
 
 !!! For orbit parameters
@@ -98,4 +100,5 @@
       use orbit_param_criterion_mod, only : orbit_param_criterion
       use recomp_orb_param_mod, only: recomp_orb_param
+      use ice_table_mod, only: computeice_table_equilibrium
 
 
@@ -189,5 +192,5 @@
       LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
       LOGICAL :: STOPPING_pressure
-      INTEGER :: criterion_stop            !Which criterion is reached : 1=surf h20, 2=surf_co2, 3=ps, 4=orb.param
+      INTEGER :: criterion_stop  ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
 
 
@@ -216,4 +219,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
       REAL ,allocatable :: watercap_slope(:,:)                           ! Physics x Nslope: watercap per slope 
+      REAL ,allocatable :: watercap_slope_saved                          ! Value saved at the previous time step
       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
@@ -229,5 +233,5 @@
       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
       REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
-      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial co2 ice 
+      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
       REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
       REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
@@ -257,6 +261,4 @@
      REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
      REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
-     REAL,ALLOCATABLE  :: alph_locslope(:,:)  ! Physic x Soil: Intermediate to compute Tsoil [1]
-     REAL,ALLOCATABLE  :: beta_locslope(:,:)  ! Physic x Soil : Intermediate tocompute Tsoil [K]
      REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
      REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
@@ -278,5 +280,7 @@
     INTEGER ::  year_iter                ! number of iteration
     INTEGER :: year_iter_max             ! maximum number of iterations before stopping
-    REAL :: timestep                     ! timestep [s] 
+    REAL :: timestep                     ! timestep [s]
+    REAL :: watercap_sum                 ! total mass of water cap [kg/m^2] 
+
     REAL WC_sum
 
@@ -424,9 +428,10 @@
        enddo
      enddo
-
+   
      call surfini(ngrid,qsurf)
+     call surfini(ngrid,qsurf)
 
 #else
-         call phys_state_var_init(nq)
+        call phys_state_var_init(nq)
          IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
          call initracer(ngrid,nq)
@@ -615,5 +620,6 @@
       call end_comsoil_h_PEM
       call ini_comsoil_h_PEM(ngrid,nslope)
-
+      call end_adsorption_h_PEM
+      call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
       allocate(ice_depth(ngrid,nslope))
       ice_depth(:,:) = 0.
@@ -800,13 +806,15 @@
       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
       watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
-      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded)
-
-      do ig=1,ngrid
-         do islope=1,nslope
-           qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)
-         enddo
+      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
+
+    do ig = 1,ngrid
+      do islope = 1,nslope
+          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
       enddo
-
-    if(soil_pem) then
+    enddo
+
+
+
+    if(adsorption_pem) then
      totmassco2_adsorbded = 0.
      totmassh2o_adsorbded = 0.
@@ -878,5 +886,5 @@
        print *, 'Global average pressure new time step',global_ave_press_new
      
-     if(soil_pem) then
+     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
@@ -986,5 +994,5 @@
 
 !------------------------
-
+     
 ! II_d.1 Update Tsurf
 
@@ -1063,7 +1071,4 @@
   allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
   allocate(Tsurf_locslope(ngrid))
-  allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
-  allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
-
   print *,"Updating soil temperature"
 
@@ -1071,32 +1076,26 @@
   do islope = 1,nslope
      TI_locslope(:,:) = TI_PEM(:,:,islope)
-     alph_locslope(:,:) = alph_PEM(:,:,islope)
-     beta_locslope(:,:) = beta_PEM(:,:,islope)
      do t = 1,timelen
        Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
        Tsurf_locslope(:) =  tsurf_phys_GCM_timeseries(:,islope,t)
-
-       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
-       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
-
-
+       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
        tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
        do ig = 1,ngrid
          do isoil = 1,nsoilmx_PEM
-          watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
+          watersoil_density_phys_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
           if(isnan(Tsoil_locslope(ig,isoil))) then
-             write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
-             stop 
+            call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1) 
           endif
          enddo
        enddo
-     enddo
-     alph_PEM(:,:,islope) = alph_locslope(:,:) 
-     beta_PEM(:,:,islope) = beta_locslope(:,:) 
+
+
+     enddo
   enddo
-
- 
      tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
      watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
+
+
   print *, "Update of soil temperature done"
 
@@ -1104,23 +1103,21 @@
   deallocate(Tsoil_locslope)
   deallocate(Tsurf_locslope)
-  deallocate(alph_locslope)
-  deallocate(beta_locslope)
-
   write(*,*) "Compute ice table"
 
 ! II_d.3 Update the ice table
-     call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
+     call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
        
       print *, "Update soil propreties"
 ! II_d.4 Update the soil thermal properties
-      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
+      call update_soil(ngrid,nslope,nsoilmx,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
         ice_depth,TI_PEM)
 
 ! II_d.5 Update the mass of the regolith adsorbded
-
+     if(adsorption_pem) then
     call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
                              tsoil_PEM,TI_PEM,ps_phys_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
                              h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
-
+ 
+     endif
   endif !soil_pem
 
@@ -1151,13 +1148,17 @@
 
 !------------------------
-      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
-
-      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
+      call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
+
+      call criterion_co2_stop(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid, &
+                                   initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
 
       year_iter=year_iter+dt_pem
+
+      
 
       print *, "Checking all the stopping criterion."
       if (STOPPING_water) then
-        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 
+        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
+        criterion_stop=1
         criterion_stop=1
       endif
@@ -1183,12 +1184,12 @@
       endif
 
-
-
-      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
+     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure)  then
         exit
       else
         print *, "We continue!"
-        print *, "Number of iteration of the PEM : year_iter=", year_iter  
+        print *, "Number of iteration of the PEM : year_iter=", year_iter
       endif
+
+
 
       global_ave_press_old=global_ave_press_new
@@ -1221,4 +1222,77 @@
 #endif
       ENDDO ! of DO ig=1,ngrid
+
+
+! H2O ice
+
+
+
+     DO ig=1,ngrid
+       if(watercaptag(ig)) then
+        
+         watercap_sum=0.
+         DO islope=1,nslope
+           watercap_slope_saved = watercap_slope(ig,islope)
+           if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then
+             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
+           else
+             watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
+             qsurf_slope(ig,igcm_h2o_ice,islope)=0.
+           endif
+           watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+!           watercap_slope(ig,islope)=0.
+         enddo
+           water_reservoir(ig)=water_reservoir(ig)+watercap_sum
+       endif
+     enddo
+
+
+      DO ig = 1,ngrid
+        qsurf(ig,igcm_h2o_ice) = 0.
+        DO islope = 1,nslope
+          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
+                                   * subslope_dist(ig,islope) /          &
+                                  cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
+     DO ig=1,ngrid
+!       DO islope=1,nslope
+         if(qsurf(ig,igcm_h2o_ice).GT.500) then
+            watercaptag(ig)=.true.
+       DO islope=1,nslope
+            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
+            water_reservoir(ig)=water_reservoir(ig)+250
+       ENDDO
+         endif
+!       enddo
+     enddo
+
+     DO ig=1,ngrid
+         if(water_reservoir(ig).LT. 10) then
+            watercaptag(ig)=.false.
+            qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
+       DO islope=1,nslope
+            qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
+       ENDDO
+         endif
+
+     enddo
+
+      DO ig = 1,ngrid
+        watercap(ig) = 0.
+        DO islope = 1,nslope
+          watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
+                                   * subslope_dist(ig,islope) /          &
+                                  cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
+
+
+
+
+      
+
 
 ! III_a.2  Tsoil update (for startfi)
@@ -1460,9 +1534,11 @@
 
 !------------------------
-
-         call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
+        call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
+                         float(day_ini),0.,nslope,def_slope,subslope_dist)
+
+
+        call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
                       tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
-
-     call info_run_PEM(year_iter, criterion_stop)
+        call info_run_PEM(year_iter, criterion_stop)
 
       print *, "restartfi_PEM.nc has been written"
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2888)
@@ -1,15 +1,16 @@
    SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
                       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)
+                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
+                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    
    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM, water_reservoir
+   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM,water_reservoir_nom
    use comsoil_h, only:  volcapa,inertiedat 
-   use adsorption_mod, only : regolith_adsorption
+   use adsorption_mod, only : regolith_adsorption,adsorption_pem
    USE temps_mod_evol, ONLY: year_PEM
+   USE ice_table_mod, only: computeice_table_equilibrium
 #ifndef CPP_STD   
    use surfdat_h, only: watercaptag
 #endif
-
 
   implicit none
@@ -45,4 +46,5 @@
   real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of h2o adsorbed [kg/m^2]
   real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+  real,intent(inout) :: water_reservoir(ngrid)                        ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
 ! local
    real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
@@ -70,9 +72,10 @@
 !!! Purpose: read start_pem. Need a specific iostart_PEM
 !!!
-!!! Order: 1. Previous year of the PEM run
-!!!        2. Thermal Inertia 
-!!!        3. Soil Temperature
-!!!        4. Ice table
-!!!        5. Mass of CO2 adsorbed
+!!! Order: 0. Previous year of the PEM run
+!!!        1. Thermal Inertia 
+!!!        2. Soil Temperature
+!!!        3. Ice table
+!!!        4. Mass of CO2 & H2O adsorbed
+!!!        5. Water reservoir
 !!!
 !!! /!\ This order must be respected !
@@ -86,7 +89,7 @@
 inquire(file=filename,exist =  startpem_file)
 
-write(*,*)'Start PEM=',startpem_file
-
-
+write(*,*)'Is start PEM?',startpem_file
+
+startpem_file = .false.
 !1. Run
 
@@ -213,6 +216,14 @@
             tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
       enddo
-    enddo
-
+          do iloop = nsoil_GCM+1,nsoil_PEM
+             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
+          enddo
+      enddo
+!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+!     do iloop = 1,2000000
+!      write(*,*) 'iloop,islope',islope,iloop
+!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+!     enddo
+    
    else
 ! predictor corrector: restart from year 1 of the GCM and build the evolution of
@@ -220,8 +231,8 @@
 
     tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
-    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+    call soil_pem_routine(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_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
 
 
@@ -243,5 +254,5 @@
 
     
-ENDDO
+ENDDO ! islope
 
 print *,'PEMETAT0: SOIL TEMP  DONE'
@@ -253,6 +264,6 @@
 
 
-     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
-     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
+     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
+     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
 
 
@@ -262,4 +273,5 @@
 
 !4. CO2 & H2O Adsorption
+ if(adsorption_pem) then
   DO islope=1,nslope
    write(num,fmt='(i2.2)') islope
@@ -290,6 +302,27 @@
     endif
 print *,'PEMETAT0: CO2 & H2O adsorption done  '
-
+    endif
 endif ! soil_pem
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!. 5 water reservoir
+
+#ifndef CPP_STD   
+   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'
+      do ig=1,ngrid
+        if(watercaptag(ig)) then
+           water_reservoir=water_reservoir_nom
+        else
+           water_reservoir=0.
+        endif
+      enddo
+   endif
+#endif
+
 
 call close_startphy
@@ -396,6 +429,16 @@
             tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
       enddo
-     enddo
-    
+     
+           do iloop = nsoil_GCM+1,nsoil_PEM
+             tsoil_PEM(ig,iloop,islope) =  tsoil_PEM(ig,nsoil_GCM,islope)
+          enddo
+       enddo
+!     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+!     do iloop = 1,2000000
+!      write(*,*) 'islope,iloop',islope,iloop
+!      call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+!     enddo
+  
+   
     do it = 1,timelen
         do isoil = nsoil_GCM+1,nsoil_PEM
@@ -412,6 +455,6 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table    
-     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
-     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
+     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
+     call update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
 
 
@@ -425,5 +468,5 @@
 !d) Regolith adsorbed
 
-
+ if(adsorption_pem) then
     m_co2_regolith_phys(:,:,:) = 0.
     m_h2o_regolith_phys(:,:,:) = 0.
@@ -434,18 +477,31 @@
     deltam_co2_regolith_phys(:) = 0.
     deltam_h2o_regolith_phys(:) = 0.  
+ endif
+
+print *,'PEMETAT0: CO2 adsorption done  '
+
+endif !soil_pem
+
+endif ! of if (startphy_file)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
 
-print *,'PEMETAT0: CO2 adsorption done  '
-
-endif !soil_pem
-
-endif ! of if (startphy_file)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-     if(soil_pem) then
-
-   
-
+
+!. e) water reservoir
+
+#ifndef CPP_STD   
+
+      do ig=1,ngrid
+        if(watercaptag(ig)) then
+           water_reservoir=water_reservoir_nom
+        else
+           water_reservoir=0.
+        endif
+      enddo
+
+#endif
+
+if(soil_pem) then
 !! Sanity check
     DO ig = 1,ngrid
@@ -453,12 +509,10 @@
         DO iloop = 1,nsoil_PEM
            if(isnan(tsoil_PEM(ig,iloop,islope))) then 
-         write(*,*) "failed nan construction", ig, iloop, islope
-           stop
+         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1) 
            endif
         ENDDO
       ENDDO
      ENDDO
-
-     endif!soil_pem
+endif!soil_pem
 
 
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2888)
@@ -64,5 +64,4 @@
   call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
 
-  call put_var("FluxGeo","Geothermal Flux (W/m^2)",fluxgeo)
 
   ! Close file
@@ -74,4 +73,6 @@
                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
                     m_co2_regolith,m_h2o_regolith,water_reservoir)
+
+
   ! write time-dependent variable to restart file
   use iostart_PEM, only : open_restartphy, close_restartphy, & 
@@ -96,5 +97,5 @@
   real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
   real,intent(in) :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
-  real,intent(IN) :: water_reservoir(ngrid)
+  real,intent(in) :: water_reservoir(ngrid)
   integer :: islope
   CHARACTER*2 :: num  
@@ -116,5 +117,6 @@
   ! set time counter in file
   call put_var("Time","Year of simulation",year_tot)
-
+  call put_field('water_reservoir','water_reservoir', &
+           water_reservoir,year_tot)
   call put_field("water_reservoir","water_reservoir", &
                  water_reservoir,year_tot)
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2888)
@@ -69,6 +69,5 @@
   REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                               ! co2 ice, mesh averaged, of the concatenated file [kg/m^2]
   REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! h2o ice per slope of the concatenated file [kg/m^2]
-  REAL, ALLOCATABLE :: watercap_slope(:,:,:,:)
-
+  REAL, ALLOCATABLE ::  watercap_slope(:,:,:,:)
 !-----------------------------------------------------------------------
   modname="read_data_gcm"
@@ -83,4 +82,5 @@
       allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
       allocate(h2o_ice_s(iim+1,jjm+1,timelen))
+      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
 
   print *, "Opening ", fichnom, "..."
@@ -135,13 +135,13 @@
 
      print *, "Downloading data for h2o_ice_s_slope done"
+
      print *, "Downloading data for watercap_slope ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
-ENDDO
-
+DO islope=1,nslope
+       write(num,fmt='(i2.2)') islope
+       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))    
+ENDDO            
      print *, "Downloading data for watercap_slope done"
-     print *, "Downloading data for tsurf_slope ..."
+    
+ print *, "Downloading data for tsurf_slope ..."
 
 DO islope=1,nslope
Index: trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2888)
@@ -1,10 +1,9 @@
       subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
-               therm_i,                          &
-               timestep,tsurf,tsoil,alph_PEM,beta_PEM)
+               therm_i,timestep,tsurf,tsoil)
 
 
       use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
                           mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
-                          coefd_PEM, mu_PEM,fluxgeo
+                          coefd_PEM, mu_PEM,alph_PEM,beta_PEM,fluxgeo
       use comsoil_h,only: volcapa
       implicit none
@@ -34,7 +33,4 @@
 ! outputs:
       real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
-      real,intent(inout) :: alph_PEM(ngrid,nsoil-1) ! intermediate for computations [1]
-      real,intent(inout) :: beta_PEM(ngrid,nsoil-1) ! intermediate for computations [K]
-      
 ! local variables:
       integer ig,ik    
Index: trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2888)
@@ -5,5 +5,6 @@
   INTEGER   year_bp_ini     !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
   INTEGER   dt_pem          !     dt_pem  : in years, the time step used by the PEM    
-  REAL      alpha_criterion !     alpha_criterion : percentage of change before stopping the PEM 
+  REAL      ice_criterion   !     ice_criterion : percentage of change of ice before stopping the PEM 
+  REAL      ps_criterion    !     ice_criterion : percentage of change of ice before stopping the PEM 
   INTEGER   year_PEM        !     year written in startfiPEM.nc
   INTEGER   Max_iter_pem    !     Maximal number of iteration when converging to a steady state, read in evol.def
Index: trunk/LMDZ.COMMON/libf/evolution/update_soil.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2886)
+++ trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2888)
@@ -1,4 +1,3 @@
-   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,&
-                          ice_depth,TI_PEM)
+   SUBROUTINE update_soil(ngrid,nslope,nsoil_GCM,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,TI_PEM)
 #ifndef CPP_STD
  USE comsoil_h, only:  inertiedat, volcapa
@@ -8,11 +7,9 @@
  implicit none
 ! Input: 
- INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
- REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
+ INTEGER,INTENT(IN) ::  ngrid, nslope,nsoil_GCM, nsoil_PEM
  REAL,INTENT(IN) :: p_avg_new
- REAL,INTENT(IN) :: co2ice(ngrid,nslope)
+ REAL,INTENT(IN) :: tendencies_waterice(ngrid,nslope)
  REAL,INTENT(IN) :: waterice(ngrid,nslope)
  REAL,INTENT(in) :: ice_depth(ngrid,nslope)
-
 ! Outputs:
 
@@ -21,17 +18,10 @@
 ! Constants:
 
- REAL ::  alpha = 0.2
- REAL ::  beta = 1.08e7
- REAL ::  To = 273.15
- REAL ::  R =  8.314
- REAL ::  L =  51058.
  REAL ::  inertie_thresold = 800. ! look for ice
  REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
  REAL ::  ice_inertia = 1200  ! Inertia of ice
- REAL ::  P610 = 610.0
- REAL ::  m_h2o = 18.01528E-3
- REAL ::  m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
- REAL ::  m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
- REAL ::  A,B,mmean             
+ REAL ::  P610 = 610.0       ! current average pressure on Mars [Pa]
+ REAL :: TI_breccia = 750.   ! THermal inertia of Breccia following Wood 2009 [SI]
+ REAL :: TI_bedrock = 2300.  ! Thermal inertia of Bedrock following Wood 2009 [SI]
 
 ! Local variables:
@@ -39,29 +29,6 @@
  INTEGER :: ig,islope,iloop,iref,k
  REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith 
- REAL :: d(ngrid,nsoil_PEM,nslope)
- REAL :: Total_surface
- INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
- INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
-
-! 0. Initialisation
-
-!  do ig = 1,ngrid
-!    do islope = 1,nslope
-!     if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then
-!        ispermanent_h2oglaciers(ig,islope) = 1
-!     else
-!        ispermanent_h2oglaciers(ig,islope) = 0
-!     endif 
-!
-!     if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then
-!        ispermanent_co2glaciers(ig,islope) = 1
-!     else
-!        ispermanent_co2glaciers(ig,islope) = 0
-!     endif
-!    enddo
-!  enddo
-
-
-
+ REAL :: delta
+ REAL :: TI_breccia_new
 !  1.Ice TI feedback
 
@@ -72,17 +39,56 @@
 ! 2. Modification of the regolith thermal inertia.
 
-!  do ig=1,ngrid
-!    do islope=1,nslope
-!     do iloop =1,n_1km
-!       d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
-!   if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
-!          TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))
-!
-!   endif 
-!     enddo
-!   enddo
-!enddo
 
-!  3. Build new TI for the PEM
+
+
+ do islope = 1,nslope
+   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
+   do ig = 1,ngrid
+
+      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
+              regolith_inertia(ig,islope) = inertie_averaged
+       endif
+        write(*,*) 'ig,islope',ig,islope,inertie_thresold,TI_PEM(ig,1,islope)
+       if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
+!          regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
+       endif 
+       TI_breccia_new = TI_breccia !*(p_avg_new/P610)**0.3
+   enddo
+ enddo
+
+
+! 3. Build new Thermal Inertia 
+
+do  islope=1,nslope
+   do ig = 1,ngrid
+     do iloop = 1,nsoil_GCM
+         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
+     enddo
+     if(regolith_inertia(ig,islope).lt.TI_breccia_new) then 
+!!! transition
+             delta = 50.
+             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
+            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
+                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia_new**2))))
+            do iloop=nsoil_GCM+2,n_1km
+              TI_PEM(ig,iloop,islope) = TI_breccia_new
+            enddo      
+      else ! we keep the high ti values
+            do iloop=nsoil_GCM+1,n_1km
+              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+            enddo 
+       endif ! TI PEM and breccia comparison
+!! transition
+       delta = 1000.
+       TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
+            ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+       do iloop=n_1km+2,nsoil_PEM
+            TI_PEM(ig,iloop,islope) = TI_bedrock
+       enddo   
+      enddo ! ig
+ENDDO ! islope
+
+!  4. Build new TI in case of ice table
 ! a) For the regolith
   do ig=1,ngrid
