Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3553)
@@ -516,2 +516,5 @@
 == 10/12/2024 == JBC
 Using the correct subroutine to output 'ice_porefilling'.
+
+== 16/12/2024 == JBC
+Addition of the shifting of the soil temperature profile to follow the surface evolution due to ice condensation/sublimation + small cleanings/improvements.
Index: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3553)
@@ -216,3 +216,119 @@
 END SUBROUTINE ini_tsoil_pem
 
+!=======================================================================
+
+SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil)
+
+use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
+use math_mod,      only: solve_steady_heat
+
+implicit none
+
+!-----------------------------------------------------------------------
+!  Author: JBC
+!  Purpose: Shifting the soil temperature profile to follow the surface evolution due to ice condensation/sublimation
+!-----------------------------------------------------------------------
+! Inputs:
+! -------
+integer,                       intent(in) :: ngrid       ! number of (horizontal) grid-points
+integer,                       intent(in) :: nsoil       ! number of soil layers
+integer,                       intent(in) :: nslope      ! number of sub-slopes
+real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m]
+real, dimension(ngrid,nslope), intent(in) :: zlag        ! newly built lag thickness [m]
+real, dimension(ngrid,nslope), intent(in) :: tsurf       ! surface temperature [K]
+! Outputs:
+! --------
+real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
+! Local:
+! ------
+integer                             :: ig, isoil, islope, ishift, ilag, iz
+real                                :: a, z, zshift_surfloc
+real, dimension(ngrid,nsoil,nslope) :: tsoil_old
+
+write(*,*)"Shifting soil temperature to surface"
+tsoil_old = tsoil
+
+do ig = 1,ngrid
+    do islope = 1,nslope
+        zshift_surfloc = zshift_surf(ig,islope)
+
+        if (zshift_surfloc >= 0.) then ! In case of the surface is higher than initially
+            if (zshift_surfloc < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
+                ! Nothing to do; we keep the soil temperature profile
+            else if (zshift_surfloc >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
+                tsoil(ig,:,islope) = tsurf(ig,islope)
+            else
+                ishift = 0
+                do while (mlayer_PEM(ishift) <= zshift_surfloc)
+                    ishift = ishift + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
+                enddo
+                ! The "new soil" temperature is set to tsurf
+                tsoil(ig,:ishift,islope) = tsurf(ig,islope)
+
+                do isoil = ishift + 1,nsoil
+                    ! Position in the old discretization of the depth
+                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
+                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
+                    iz = 0
+                    do while (mlayer_PEM(iz) <= z)
+                        iz = iz + 1
+                    enddo
+                    ! Interpolation of the temperature profile from the old discretization
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
+                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                enddo
+            endif
+        else ! In case of the surface is lower than initially
+            if (abs(zshift_surfloc) < mlayer_PEM(0)) then ! Surface change is too small to be taken into account
+                ! Nothing to do; we keep the soil temperature profile
+            else if (abs(zshift_surfloc) >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization
+                call solve_steady_heat(nsoil,mlayer_PEM,layer_PEM,mthermdiff_PEM(ig,:),thermdiff_PEM(ig,:),tsurf(ig,islope),fluxgeo,tsoil(ig,:,islope))
+            else
+                if (zlag(ig,islope) < mlayer_PEM(0)) then ! The lag is too thin to be taken into account
+                    ilag = 0
+                else
+                    ilag = 0
+                    do while (mlayer_PEM(ilag) <= zlag(ig,islope))
+                        ilag = ilag + 1 ! mlayer indices begin at 0 so this the good index for tsoil!
+                    enddo
+                    ! Position of the lag bottom in the old discretization of the depth
+                    z = zlag(ig,islope) - zshift_surfloc
+                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
+                    iz = 0
+                    do while (mlayer_PEM(iz) <= z)
+                        iz = iz + 1
+                    enddo
+                    ! The "new lag" temperature is set to the ice temperature just below
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
+                    tsoil(ig,:ilag,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                endif
+
+                ishift = nsoil - 1
+                z = mlayer_PEM(nsoil - 1) + zshift_surfloc
+                do while (mlayer_PEM(ishift) >= z)
+                    ishift = ishift - 1
+                enddo
+                ishift = ishift + 1 ! Adding 1 is needed to match the good index for tsoil!
+                do isoil = ilag + 1,ishift
+                    ! Position in the old discretization of the depth
+                    z = mlayer_PEM(isoil - 1) - zshift_surfloc
+                    ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
+                    iz = 0
+                    do while (mlayer_PEM(iz) <= z)
+                        iz = iz + 1
+                    enddo
+                    ! Interpolation of the temperature profile from the old discretization
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
+                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                enddo
+
+                ! The "new deepest layers" temperature is set by solving the steady heat equation
+                call solve_steady_heat(nsoil - ishift + 1,mlayer_PEM(ishift - 1:),layer_PEM(ishift:),mthermdiff_PEM(ig,ishift - 1:),thermdiff_PEM(ig,ishift:),tsoil(ig,ishift,islope),fluxgeo,tsoil(ig,ishift:,islope))
+            endif
+        endif
+    enddo
+enddo
+
+END SUBROUTINE shift_tsoil2surf
+
 END MODULE compute_soiltemp_mod
Index: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3553)
@@ -7,7 +7,8 @@
 !=======================================================================
 
-SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
+SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
 
 use time_evol_mod, only: dt
+use glaciers_mod,  only: rho_co2ice
 
 implicit none
@@ -18,15 +19,16 @@
 !
 !=======================================================================
-!   arguments:
-!   ----------
-!   INPUT
+! 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) :: d_co2ice ! Evolution of perennial ice over one year
+! OUTPUT
+real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! Previous and actual density of CO2 ice
+real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of perennial ice over one year
+real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
 
-!   local:
-!   ------
+! local:
+! ------
 real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
 !=======================================================================
@@ -41,12 +43,15 @@
 end where
 
+zshift_surf = d_co2ice*dt/rho_co2ice
+
 END SUBROUTINE evol_co2_ice
 
 !=======================================================================
 
-SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
+SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
 
-use time_evol_mod,          only: dt
-use comslope_mod,           only: subslope_dist, def_slope_mean
+use time_evol_mod, only: dt
+use comslope_mod,  only: subslope_dist, def_slope_mean
+use glaciers_mod,  only: rho_h2oice
 
 #ifndef CPP_STD
@@ -63,7 +68,7 @@
 !
 !=======================================================================
-!   arguments:
-!   ----------
-!   INPUT
+! 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)
@@ -71,14 +76,15 @@
 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) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year)
-integer,                       intent(inout) :: stopPEM      ! Stopping criterion code
+! 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) :: d_h2oice    ! Evolution of perennial ice over one year (kg/m^2/year)
+integer,                       intent(inout) :: stopPEM     ! Stopping criterion code
+real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
 
-!   local:
-!   ------
+! local:
+! ------
 integer                       :: i, islope                                           ! Loop variables
 real                          :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o
-real, dimension(ngrid,nslope) :: new_tend                                      ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
+real, dimension(ngrid,nslope) :: new_tend                                            ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
 !=======================================================================
 if (ngrid /= 1) then ! Not in 1D
@@ -164,4 +170,6 @@
 endif
 
+zshift_surf = d_h2oice*dt/rho_h2oice
+
 END SUBROUTINE evol_h2o_ice
 
Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3553)
@@ -10,7 +10,10 @@
 
 ! Thresholds for ice management
-real, save :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
-real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2]
-real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2]
+real :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
+real :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2]
+real :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2]
+
+real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
+real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
 
 !=======================================================================
Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3553)
@@ -8,4 +8,5 @@
 !=======================================================================
 
+use glaciers_mod, only: rho_co2ice, rho_h2oice
 
 implicit none
@@ -23,6 +24,4 @@
 real, parameter :: co2ice_porosity = 0.3        ! Porosity of CO2 ice
 real, parameter :: h2oice_porosity = 0.3        ! Porosity of H2O ice
-real, parameter :: rho_co2ice = 1650.           ! Density of CO2 ice [kg.m-3]
-real, parameter :: rho_h2oice = 920.            ! Density of H2O ice [kg.m-3]
 real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
 real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]
@@ -452,14 +451,18 @@
 !=======================================================================
 ! Layering algorithm
-SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
-
-implicit none
-
-!---- Arguments
+SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,zshift_surf,new_lag1,new_lag2,zlag,stopPEM,current1,current2)
+
+implicit none
+
+!---- Arguments
+! Inputs
+real, intent(in) :: d_co2ice, d_h2oice, d_dust
+
+! Outputs
 type(layering),         intent(inout) :: this
 type(stratum), pointer, intent(inout) :: current1, current2
 logical,                intent(inout) :: new_str, new_lag1, new_lag2
 integer,                intent(inout) :: stopPEM
-real,     intent(in) :: d_co2ice, d_h2oice, d_dust
+real, intent(out) :: zshift_surf, zlag
 
 !---- Local variables
@@ -474,4 +477,6 @@
 dh_h2oice = d_h2oice/rho_h2oice
 dh_dust = d_dust/rho_dust
+zshift_surf = 0.
+zlag = 0.
 
 if (dh_dust < 0.) then ! Dust lifting only
@@ -504,4 +509,5 @@
         stopPEM = 8
     endif
+    zshift_surf = dh_dust + h2lift_tot
 
 else
@@ -521,4 +527,5 @@
              endif
         endif
+        zshift_surf = thickness
 
 !------ Condensation of CO2 ice + H2O ice
@@ -536,4 +543,5 @@
              endif
         endif
+        zshift_surf = thickness
 
 !------ Sublimation of CO2 ice + Condensation of H2O ice
@@ -563,5 +571,5 @@
                 h2subl = h2subl_tot
                 h2subl_tot = 0.
-                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
@@ -575,5 +583,5 @@
                 h2subl = h_co2ice_old
                 h2subl_tot = h2subl_tot - h2subl
-                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
                 current1 => current1%down
                 new_lag1 = .true.
@@ -595,4 +603,5 @@
              endif
         endif
+        zshift_surf = zshift_surf + thickness
 
 !------ Sublimation of CO2 ice + H2O ice
@@ -622,5 +631,5 @@
                 h2subl = h2subl_tot
                 h2subl_tot = 0.
-                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
             else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
@@ -634,5 +643,5 @@
                 h2subl = h_co2ice_old
                 h2subl_tot = h2subl_tot - h2subl
-                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag)
                 current1 => current1%down
                 new_lag1 = .true.
@@ -666,5 +675,5 @@
                 h2subl = h2subl_tot
                 h2subl_tot = 0.
-                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
+                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag)
             else if (h_h2oice_old < eps) then ! There is nothing to sublimate
                 ! Is the considered stratum a dust lag?
@@ -678,5 +687,5 @@
                 h2subl = h_h2oice_old
                 h2subl_tot = h2subl_tot - h2subl
-                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
+                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag)
                 current2 => current2%down
                 new_lag2 = .true.
@@ -698,5 +707,5 @@
 !=======================================================================
 ! To sublimate CO2 ice in the layering
-SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
+SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
 
 implicit none
@@ -706,4 +715,5 @@
 type(stratum), pointer, intent(inout) :: str
 logical,                intent(inout) :: new_lag
+real,                   intent(inout) :: zshift_surf, zlag
 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
 
@@ -719,4 +729,5 @@
 
 str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
+zshift_surf = zshift_surf - h2subl/(1. - co2ice_porosity) - hsubl_dust
 if (str%thickness < tol) then
     ! Remove the sublimating stratum if there is no ice anymore
@@ -748,4 +759,5 @@
     endif
 endif
+zlag = zlag + h_lag
 
 END SUBROUTINE subl_co2ice_layering
@@ -753,5 +765,5 @@
 !=======================================================================
 ! To sublimate H2O ice in the layering
-SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)
+SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
 
 implicit none
@@ -761,4 +773,5 @@
 type(stratum), pointer, intent(inout) :: str
 logical,                intent(inout) :: new_lag
+real,                   intent(inout) :: zshift_surf, zlag
 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
 
@@ -774,4 +787,5 @@
 
 str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
+zshift_surf = zshift_surf - h2subl/(1. - h2oice_porosity) - hsubl_dust
 if (str%thickness < tol) then
     ! Remove the sublimating stratum if there is no ice anymore
@@ -803,4 +817,5 @@
     endif
 endif
+zlag = zlag + h_lag
 
 END SUBROUTINE subl_h2oice_layering
Index: trunk/LMDZ.COMMON/libf/evolution/math_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/math_mod.F90	(revision 3553)
@@ -28,10 +28,15 @@
 implicit none
 
-integer, intent(in) :: nz      ! number of layer
-real,    intent(in) :: z(nz)   ! depth layer
-real,    intent(in) :: y(nz)   ! function which needs to be derived
-real,    intent(in) :: y0,ybot ! boundary conditions
-real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth
-! local
+! Inputs
+!-------
+integer,             intent(in) :: nz       ! number of layer
+real, dimension(nz), intent(in) :: z        ! depth layer
+real, dimension(nz), intent(in) :: y        ! function which needs to be derived
+real,                intent(in) :: y0, ybot ! boundary conditions
+! Outputs
+!--------
+real, dimension(nz), intent(out) :: dzY ! derivative of y w.r.t depth
+! Local variables
+!----------------
 integer :: j
 real    :: hm, hp, c1, c2, c3
@@ -69,7 +74,14 @@
 implicit none
 
-integer, intent(in) :: nz
-real,    intent(in) :: z(nz), y(nz), y0, yNp1
-real, intent(out) :: yp2(nz)
+! Inputs
+!-------
+integer,             intent(in) :: nz
+real, dimension(nz), intent(in) :: z, y
+real,                intent(in) :: y0, yNp1
+! Outputs
+!--------
+real, dimension(nz), intent(out) :: yp2
+! Local variables
+!----------------
 integer :: j
 real    :: hm, hp, c1, c2, c3
@@ -102,9 +114,15 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  implicit none
-
-integer, intent(in) :: nz, j
-real,    intent(in) :: z(nz), y(nz)
+implicit none
+
+! Inputs
+!-------
+integer,             intent(in) :: nz, j
+real, dimension(nz), intent(in) :: z, y
+! Outputs
+!--------
 real, intent(out) :: dy_zj
+! Local viariables
+!-----------------
 real :: h1, h2, c1, c2, c3
 
@@ -135,9 +153,15 @@
 implicit none
 
-integer, intent(in) :: nz, i1, i2
-real,    intent(in) :: y(nz), z(nz)
+! Inputs
+!-------
+integer,             intent(in) :: nz, i1, i2
+real, dimension(nz), intent(in) :: y, z
+! Outputs
+!--------
 real, intent(out) :: integral
-integer i
-real dz(nz)
+! Local viariables
+!-----------------
+integer             :: i
+real, dimension(nz) :: dz
 
 dz(1) = (z(2) - 0.)/2
@@ -163,7 +187,11 @@
 implicit none
 
+! Inputs
+!-------
 real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
 real, intent(in) :: z1, z2 ! depth [m]
-real, intent(out) :: zr    ! depth at which we have zero
+! Outputs
+!--------
+real, intent(out) :: zr ! depth at which we have zero
 
 zr = (y1*z2 - y2*z1)/(y1 - y2)
@@ -171,3 +199,136 @@
 END SUBROUTINE findroot
 
+!=======================================================================
+
+SUBROUTINE solve_tridiag(a,b,c,d,n,x,error)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Solve a tridiagonal system Ax = d using the Thomas' algorithm
+!!!          a: sub-diagonal
+!!!          b: main diagonal
+!!!          c: super-diagonal
+!!!          d: right-hand side
+!!!          x: solution
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+
+! Inputs
+!-------
+integer,                intent(in) :: n
+real, dimension(n),     intent(in) :: b, d
+real, dimension(n - 1), intent(in) :: a, c
+! Outputs
+!--------
+real, dimension(n), intent(out) :: x
+integer,            intent(out) :: error
+! Local viariables
+!-----------------
+integer            :: i
+real               :: m
+real, dimension(n) :: cp, dp
+
+! Check stability: diagonally dominant condition
+error = 0
+if (abs(b(1)) < abs(c(1))) then
+    error = 1
+    return
+endif
+do i = 2,n - 1
+    if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then
+        error = 1
+        return
+    endif
+enddo
+if (abs(b(n)) < abs(a(n - 1))) then
+    error = 1
+    return
+endif
+
+! Initialization
+cp(1) = c(1)/b(1)
+dp(1) = d(1)/b(1)
+
+! Forward sweep
+do i = 2,n - 1
+    m = b(i) - a(i - 1)*cp(i - 1)
+    cp(i) = c(i)/m
+    dp(i) = (d(i) - a(i - 1)*dp(i - 1))/m
+enddo
+m = b(n) - a(n - 1)*cp(n - 1)
+dp(n) = (d(n) - a(n - 1)*dp(n - 1))/m
+
+! Backward substitution
+x(n) = dp(n)
+do i = n - 1,1,-1
+    x(i) = dp(i) - cp(i)*x(i + 1)
+enddo
+
+END SUBROUTINE solve_tridiag
+
+!=======================================================================
+
+SUBROUTINE solve_steady_heat(n,z,mz,kappa,mkappa,T_left,q_right,T)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Solve 1D steady-state heat equation with space-dependent diffusivity
+!!!          Left boudary condition : prescribed temperature T_left
+!!!          Right boudary condition: prescribed thermal flux q_right
+!!!
+!!!          z     : grid points
+!!!          mz    : mid-grid points
+!!!          kappa : thermal diffusivity at grid points
+!!!          mkappa: thermal diffusivity at mid-grid points
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use abort_pem_mod, only: abort_pem
+
+implicit none
+
+! Inputs
+!-------
+integer,                intent(in) :: n
+real, dimension(n),     intent(in) :: z, kappa
+real, dimension(n - 1), intent(in) :: mz, mkappa
+real,                   intent(in) :: T_left, q_right
+! Outputs
+!--------
+real, dimension(n), intent(out) :: T
+! Local viariables
+!-----------------
+integer                :: i, error
+real, dimension(n)     :: b, d
+real, dimension(n - 1) :: a, c
+
+! Initialization
+a = 0.; b = 0.; c = 0.; d = 0.
+
+! Left boundary condition (Dirichlet: prescribed temperature)
+b(1) = 1.
+d(1) = T_left
+
+! Internal points
+do i = 2,n - 1
+    a(i - 1) = -mkappa(i - 1)/((mz(i) - mz(i - 1))*(z(i) - z(i - 1)))
+    c(i) = -mkappa(i)/((mz(i) - mz(i - 1))*(z(i + 1) - z(i)))
+    b(i) = -(a(i - 1) + c(i))
+enddo
+
+! Right boundary condition (Neumann: prescribed temperature)
+a(n - 1) = kappa(n - 1)/(z(n) - z(n - 1))
+b(n) = -kappa(n)/(z(n) - z(n - 1))
+d(n) = q_right
+
+! Solve the tridiagonal linear system with the Thomas' algorithm
+call solve_tridiag(a,b,c,d,n,T,error)
+if (error /= 0) call abort_pem("solve_steady_heat","Unstable solving!",1)
+
+END SUBROUTINE solve_steady_heat
+
 end module math_mod
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3553)
@@ -44,5 +44,5 @@
 use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
 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, &
+use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, mlayer_PEM, &
                                       TI_PEM,               & ! soil thermal inertia
                                       tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
@@ -50,5 +50,5 @@
 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
+                                      co2_adsorbed_phys, h2o_adsorbed_phys        ! Mass of co2 and h2O adsorbed
 use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
 use orbit_param_criterion_mod,  only: orbit_param_criterion
@@ -65,6 +65,6 @@
 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 compute_soiltemp_mod,       only: compute_tsoil_pem
+use recomp_tend_co2_slope_mod,  only: recomp_tend_co2
+use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
 use writediagpem_mod,           only: writediagpem, writediagsoilpem
 use co2condens_mod,             only: CO2cond_ps
@@ -217,8 +217,8 @@
 real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
 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]
-real                                  :: totmassco2_adsorbded             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
-real                                  :: totmassh2o_adsorbded             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
+real, dimension(:),       allocatable :: delta_co2_adsorbed              ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
+real, dimension(:),       allocatable :: delta_h2o_adsorbed              ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
+real                                  :: totmassco2_adsorbed             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
+real                                  :: totmassh2o_adsorbed             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
 logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
 logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
@@ -228,4 +228,6 @@
 real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
 real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
+real, dimension(:,:),     allocatable :: zshift_surf                      ! Elevation shift for the surface [m]
+real, dimension(:,:),     allocatable :: zlag                             ! Newly built lag thickness [m]
 
 ! Some variables for the PEM run
@@ -541,12 +543,11 @@
 allocate(vmr_co2_PCM(ngrid,timelen))
 allocate(ps_timeseries(ngrid,timelen))
-allocate(min_co2_ice(ngrid,nslope,2))
-allocate(min_h2o_ice(ngrid,nslope,2))
+allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
 allocate(tsurf_avg_yr1(ngrid,nslope))
 allocate(tsurf_avg(ngrid,nslope))
 allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
 allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
-allocate(q_co2_PEM_phys(ngrid,timelen))
-allocate(q_h2o_PEM_phys(ngrid,timelen))
+allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
+allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
 allocate(watersurf_density_avg(ngrid,nslope))
 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
@@ -554,10 +555,8 @@
 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
-allocate(delta_co2_adsorbded(ngrid))
+allocate(delta_co2_adsorbed(ngrid))
 allocate(co2ice_disappeared(ngrid,nslope))
-allocate(icetable_thickness_old(ngrid,nslope))
-allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
-allocate(delta_h2o_icetablesublim(ngrid))
-allocate(delta_h2o_adsorbded(ngrid))
+allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
+allocate(delta_h2o_icetablesublim(ngrid),delta_h2o_adsorbed(ngrid))
 allocate(vmr_co2_PEM_phys(ngrid,timelen))
 
@@ -649,6 +648,6 @@
               icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
               ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
-              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,         &
-              delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
+              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,         &
+              delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
 
 ! We save the places where h2o ice is sublimating
@@ -680,18 +679,18 @@
 
 if (adsorption_pem) then
-    totmassco2_adsorbded = 0.
-    totmassh2o_adsorbded = 0.
+    totmassco2_adsorbed = 0.
+    totmassh2o_adsorbed = 0.
     do ig = 1,ngrid
         do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
                 if (l==1) then
-                   totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
+                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(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))* &
+                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                 else
-                   totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
+                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
                                        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) - layer_PEM(l-1))* &
+                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                 endif
@@ -699,6 +698,6 @@
         enddo
     enddo
-    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
-    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
+    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed
+    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
 endif ! adsorption
 deallocate(tsurf_avg_yr1)
@@ -750,5 +749,5 @@
     if (adsorption_pem) then
         do i = 1,ngrid
-            global_avg_press_new = global_avg_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_adsorbed(i)/Total_surface
         enddo
     endif
@@ -819,12 +818,14 @@
 !    II_b Evolution of ice
 !------------------------
-    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
-    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
+    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
+    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
     if (layering_algo) then
         do islope = 1,nslope
             do ig = 1,ngrid
-                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
+                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
             enddo
         enddo
+    else
+        zlag = 0.
     endif
 
@@ -887,4 +888,5 @@
                 tsurf_avg(ig,islope) = ave/timelen
             endif
+            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBLIMATION??? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         enddo
     enddo
@@ -901,6 +903,8 @@
 
     if (soil_pem) then
-
-! II_d.2 Update soil temperature
+! II_d.2 Shifting soil temperature to surface
+        call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
+
+! II_d.3 Update soil temperature
         write(*,*)"Updating soil temperature"
         allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
@@ -926,5 +930,5 @@
         deallocate(Tsoil_locslope)
 
-! II_d.3 Update the ice table
+! II_d.4 Update the ice table
         if (icetable_equilibrium) then
             write(*,*) "Compute ice table (equilibrium method)"
@@ -947,27 +951,27 @@
         endif
 
-! II_d.4 Update the soil thermal properties
+! II_d.5 Update the soil thermal properties
         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
 
-! II_d.5 Update the mass of the regolith adsorbed
+! II_d.6 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                           &
                                      h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
-                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
-
-            totmassco2_adsorbded = 0.
-            totmassh2o_adsorbded = 0.
+                                     h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
+
+            totmassco2_adsorbed = 0.
+            totmassh2o_adsorbed = 0.
             do ig = 1,ngrid
                 do islope = 1,nslope
                     do l = 1,nsoilmx_PEM
                         if (l == 1) then
-                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
+                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(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))* &
+                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         else
-                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
+                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        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) - layer_PEM(l - 1))* &
+                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         endif
@@ -975,6 +979,6 @@
                 enddo
             enddo
-            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
-            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
+            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbed
+            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbed
         endif
     endif !soil_pem
@@ -1005,6 +1009,6 @@
             if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
             if (adsorption_pem) then
-                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
-                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
+                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
+                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
             endif
         endif
@@ -1015,5 +1019,5 @@
 !    II_f Update the tendencies
 !------------------------
-    call recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new)
+    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new)
 
 !------------------------
@@ -1021,5 +1025,4 @@
 !    II_g Checking the stopping criterion
 !------------------------
-
     write(*,*) "Checking the stopping criteria..."
     call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
@@ -1210,5 +1213,5 @@
 call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
-             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
+             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif)
 write(*,*) "restartpem.nc has been written"
 
@@ -1231,6 +1234,6 @@
 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
-deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
-deallocate(icetable_thickness_old,ice_porefilling_old)
+deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
+deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag)
 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
 !----------------------------- END OUTPUT ------------------------------
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3553)
@@ -65,4 +65,5 @@
 write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
 write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp
+write(*,*) 'recomp_orb_param, New year =', Year
 
 found_year = .false.
@@ -116,5 +117,4 @@
 #endif
 
-write(*,*) 'recomp_orb_param, New year =', Year
 write(*,*) 'recomp_orb_param, New obl. =', obliquit
 write(*,*) 'recomp_orb_param, New ecc. =', e_elips
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3543)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3553)
@@ -7,6 +7,6 @@
 !=======================================================================
 
-SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_phys_ini,co2ice_slope,emissivity_slope, &
-                                 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
+SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
+                           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
@@ -24,14 +24,14 @@
 !   INPUT
 integer,                        intent(in) :: timelen, ngrid, nslope
-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_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) :: d_co2ice_phys_ini ! physical point field: Evolution of perennial ice over one year
-real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope (kg/m^2)
-real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
+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_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) :: d_co2ice_ini         ! physical point field: Evolution of perennial ice over one year
+real, dimension(ngrid,nslope),  intent(in) :: co2ice               ! CO2 ice per mesh and sub-grid slope (kg/m^2)
+real, dimension(ngrid,nslope),  intent(in) :: emissivity           ! Emissivity per mesh and sub-grid slope(1)
 !   OUTPUT
-real, dimension(ngrid,nslope), intent(inout) ::  d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
+real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
 
 !   local:
@@ -45,18 +45,55 @@
 do i = 1,ngrid
     do islope = 1,nslope
-        coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
+        coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
         ave = 0.
-        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
-            do t=1,timelen
+        if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
+            do t = 1,timelen
                 ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
-                      - (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
+                      - (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
-            if (ave < 1e-4) ave = 0.
-            d_co2ice_phys(i,islope) = d_co2ice_phys_ini(i,islope) - coef*ave/timelen
+            if (ave < 1.e-4) ave = 0.
+            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*ave/timelen
         endif
     enddo
 enddo
 
-END SUBROUTINE recomp_tend_co2_slope
+END SUBROUTINE recomp_tend_co2
+!=======================================================================
+
+SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp)
+
+implicit none
+
+!=======================================================================
+!
+!  Routine that compute the evolution of the tendencie for h2o ice
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,            intent(in) :: timelen, ngrid, nslope
+real, dimension(:), intent(in) :: PCM_temp, PEM_temp
+
+!   OUTPUT
+real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year
+
+!   local:
+!   ------
+real :: coef, ave, Rz_old, Rz_new, R_dec, soil_psv_old, soil_psv_new, hum_dec, h2oice_depth_new, h2oice_depth_old
+
+write(*,*) "Update of the H2O tendency due to lag layer"
+
+! Flux correction due to lag layer
+!~ Rz_old = h2oice_depth_old*0.0325/4.e-4              ! resistance from PCM
+!~ Rz_new = h2oice_depth_new*0.0325/4.e-4              ! new resistance based on new depth
+!~ R_dec = (1./Rz_old)/(1./Rz_new)                     ! decrease because of resistance
+!~ soil_psv_old = psv(max(PCM_temp(h2oice_depth_old))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the GCM run temperature at the old ice location
+!~ soil_psv_new = psv(max(PEM_temp(h2oice_depth_new))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the PEM run temperature at the new ice location
+!~ hum_dec = soil_psv_old/soil_psv_new                 ! decrease because of lower water vapor pressure at the new depth
+!~ d_h2oice = d_h2oice*R_dec*hum_dec                   ! decrease of flux
+
+END SUBROUTINE recomp_tend_h2o
 
 END MODULE recomp_tend_co2_slope_mod
