Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3285)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3286)
@@ -250,2 +250,10 @@
 == 13/03/2024 == LL
 Correction of small mistake when computing the total mass of CO2/H2O stored in the soil; the wrong layer index was used
+
+== 28/03/2024 == JBC
+Introduction of the module aiming to manage the layered deposits through a linked list data structure. It contains:
+    - the object 'stratum' which holds several properties for one layer as well as pointers to the lower/upper 'stratum';
+    - the object 'layering' which brings the linked 'stratum' objects together as a comprehensive structure;
+    - properties for CO2/H2O ice, regolith and dust;
+    - several subroutines to manipulate the 'layering' and 'stratum' (initilize, finalize, display, modify, etc);
+    - subroutines of the algorithm to make the layering evolve according to the input tendencies (new layer, change of content in a layer, creation of dust lag layer, etc).
Index: /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3286)
+++ /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3286)
@@ -0,0 +1,609 @@
+MODULE layering_mod
+
+!=======================================================================
+! Purpose: Manage layered deposits in the PEM with a linked list data structure
+!
+! Author: JBC
+! Date: 08/03/2024
+!=======================================================================
+
+
+implicit none
+
+!---- Declarations
+! Numerical parameter
+real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
+
+! Physical parameters
+real, parameter :: tend_dust = 5.78e-2          ! Tendency of dust [kg.m-2.y-1]
+real, parameter :: dry_lag_porosity = 0.5       ! Porosity of dust lag
+real, parameter :: dry_regolith_porosity = 0.45 ! Porosity of regolith
+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]
+
+! Stratum = node of the linked list
+type :: stratum
+    real               :: thickness          ! Layer thickness [m]
+    real               :: elevation          ! Layer elevation (top height from the surface) [m]
+    real               :: co2ice_volfrac     ! CO2 ice volumic fraction
+    real               :: h2oice_volfrac     ! H2O ice volumic fraction
+    real               :: dust_volfrac       ! Dust volumic fraction
+    real               :: air_volfrac        ! Air volumic fraction inside pores
+    type(stratum), pointer :: up => null()   ! Upper stratum (next node)
+    type(stratum), pointer :: down => null() ! Lower stratum (previous node)
+end type stratum
+
+! Layering = linked list
+type layering
+    integer                :: nb_str = 0       ! Number of strata
+    type(stratum), pointer :: bottom => null() ! First stratum (head of the list)
+    type(stratum), pointer :: top => null()    ! Last stratum (tail of the list)
+end type layering
+
+!=======================================================================
+contains
+!=======================================================================
+! To display a stratum
+SUBROUTINE print_stratum(this)
+
+implicit none
+
+!---- Arguments
+type(stratum), intent(in) :: this
+
+!---- Code
+write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
+write(*,'(a20,es13.6)') 'Elevation      : ', this%elevation
+write(*,'(a20,es13.6)') 'CO2 ice content: ', this%co2ice_volfrac
+write(*,'(a20,es13.6)') 'H2O ice content: ', this%h2oice_volfrac
+write(*,'(a20,es13.6)') 'Dust content   : ', this%dust_volfrac
+write(*,'(a20,es13.6)') 'Air content    : ', this%air_volfrac
+
+END SUBROUTINE print_stratum
+
+!=======================================================================
+! To display a layering
+SUBROUTINE print_lay(this)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(in) :: this
+
+!---- Local variables
+type(stratum), pointer :: current
+integer                :: i
+
+!---- Code
+current => this%top
+
+i = this%nb_str
+do while (associated(current))
+    write(*,'(a8,i4)') 'Stratum ', i
+    call print_stratum(current)
+    write(*,*)
+    current => current%down
+    i = i - 1
+enddo
+
+END SUBROUTINE print_lay
+
+!=======================================================================
+! To add a stratum to the top of a layering
+SUBROUTINE add_stratum(this,thickness,elevation,co2ice,h2oice,dust,air)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(inout)  :: this
+real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air
+
+!---- Local variables
+type(stratum), pointer :: str
+
+!---- Code
+! Creation of the stratum
+allocate(str)
+nullify(str%up,str%down)
+str%thickness = 1.
+str%elevation = 1.
+str%co2ice_volfrac = 0.
+str%h2oice_volfrac = 0.
+str%dust_volfrac = 0.
+str%air_volfrac = 1.
+if (present(thickness)) str%thickness = thickness
+if (present(elevation)) str%elevation = elevation
+if (present(co2ice)) str%co2ice_volfrac = co2ice
+if (present(h2oice)) str%h2oice_volfrac = h2oice
+if (present(dust)) str%dust_volfrac = dust
+if (present(air)) str%air_volfrac = air
+
+! Verification of volume fraction
+if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
+     error stop 'add_stratum: properties for the new stratum are not possible (sum of volumic fraction /= 1)!'
+
+! Increment the number of strata
+this%nb_str = this%nb_str + 1
+
+! Add the stratum to the layering
+if (.not. associated(this%bottom)) then ! If it is the very first one
+    this%bottom => str
+    this%top => str
+else
+    str%down => this%top
+    this%top%up => str
+    this%top => str
+endif
+
+END SUBROUTINE add_stratum
+
+!=======================================================================
+! To insert a stratum after another one in a layering
+SUBROUTINE insert_stratum(this,str_prev,thickness,elevation,co2ice,h2oice,dust,air)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str_prev
+real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air
+
+!---- Local variables
+type(stratum), pointer :: str, current
+
+!---- Code
+if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
+    call add_stratum(this,thickness,elevation,co2ice,h2oice,dust,air)
+else
+    ! Creation of the stratum
+    allocate(str)
+    nullify(str%up,str%down)
+    str%thickness = 1.
+    str%elevation = 1.
+    str%co2ice_volfrac = 0.
+    str%h2oice_volfrac = 0.
+    str%dust_volfrac = 0.
+    str%air_volfrac = 1.
+    if (present(thickness)) str%thickness = thickness
+    if (present(elevation)) str%elevation = elevation
+    if (present(co2ice)) str%co2ice_volfrac = co2ice
+    if (present(h2oice)) str%h2oice_volfrac = h2oice
+    if (present(dust)) str%dust_volfrac = dust
+    if (present(air)) str%air_volfrac = air
+
+    ! Verification of volume fraction
+    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
+         error stop 'insert_stratum: properties for the new stratum are not possible (sum of volumic fraction /= 1)!'
+
+    ! Increment the number of strata
+    this%nb_str = this%nb_str + 1
+
+    ! Insert the stratum into the layering at the right place
+    str%down => str_prev
+    str%up => str_prev%up
+    str_prev%up => str
+    str%up%down => str
+
+    ! Correct the value of elevation for the upper strata
+    current => str%up
+    do while (associated(current))
+        current%elevation = current%down%elevation + current%thickness
+        current => current%up
+    enddo
+endif
+
+END SUBROUTINE insert_stratum
+
+!=======================================================================
+! To remove a stratum in a layering
+SUBROUTINE rm_stratum(this,str)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+
+!---- Local variables
+type(stratum), pointer :: new_str
+
+!---- Code
+! Protect the 3 sub-surface strata from removing
+if (str%elevation <= 0.) return
+
+! Decrement the number of strata
+this%nb_str = this%nb_str - 1
+
+! Remove the stratum from the stratification
+str%down%up => str%up
+if (associated(str%up)) then ! If it is not the last one at the top of the layering
+    str%up%down => str%down
+else
+    this%top => str%down
+endif
+new_str => str%down
+nullify(str%up,str%down)
+deallocate(str)
+nullify(str)
+str => new_str
+
+END SUBROUTINE rm_stratum
+
+!=======================================================================
+! To initialize the layering
+SUBROUTINE ini_layering(this)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(inout) :: this
+
+!---- Code
+! Creation of three strata at the bottom of the layering
+! 1) Dry porous deep regolith
+call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
+! 2) Ice-cemented regolith (ice table)
+call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
+! 3) Dry porous regolith
+call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
+
+END SUBROUTINE ini_layering
+
+!=======================================================================
+! To delete the entire layering
+SUBROUTINE del_layering(this)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(inout) :: this
+
+!---- Local variables
+type(stratum), pointer :: current, next
+
+!---- Code
+current => this%bottom
+do while (associated(current))
+    next => current%up
+    deallocate(current)
+    nullify(current)
+    current => next
+enddo
+this%nb_str = 0
+
+END SUBROUTINE del_layering
+
+!=======================================================================
+! To get a specific stratum in a layering
+FUNCTION get_stratum(lay,i) RESULT(str)
+
+implicit none
+
+!---- Arguments
+type(layering), intent(in) :: lay
+integer,        intent(in) :: i
+type(stratum), pointer     :: str
+
+!---- Local variables
+integer :: k
+
+!---- Code
+if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
+k = 1
+str => lay%bottom
+do while (k /= i)
+    str => str%up
+    k = k + 1
+enddo
+
+END FUNCTION get_stratum
+
+!=======================================================================
+! 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)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+logical,                intent(inout) :: new_lag
+real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
+
+!---- Local variables
+real               :: r_subl, hsubl_dust, h_lag
+type(stratum), pointer :: current
+
+!---- Code
+! How much dust does CO2 ice sublimation involve to create the lag?
+r_subl = h2subl/(1. - co2ice_porosity) &
+        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
+hsubl_dust = r_subl*str%dust_volfrac*str%thickness
+! Update of properties in the sublimating stratum
+str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
+str%elevation = str%elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
+str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
+str%h2oice_volfrac = h_h2oice_old/str%thickness
+str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
+str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+! Correct the value of elevation for the upper strata
+current => str%up
+do while (associated(current))
+    current%elevation = current%down%elevation + current%thickness
+    current => current%up
+enddo
+! Remove the sublimating stratum if there is no ice anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+! New stratum = dust lag
+h_lag = hsubl_dust/(1. - dry_lag_porosity)
+if (h_lag > 0.) then ! Only if the dust lag is building up
+    if (new_lag) then
+        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
+        new_lag = .false.
+    else
+        str%up%thickness = str%up%thickness + h_lag
+        str%up%elevation = str%up%elevation + h_lag
+    endif
+endif
+
+END SUBROUTINE subl_co2ice_layering
+
+!=======================================================================
+! 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)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+logical,                intent(inout) :: new_lag
+real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
+
+!---- Local variables
+real               :: r_subl, hsubl_dust, h_lag
+type(stratum), pointer :: current
+
+!---- Code
+! How much dust does CO2 ice sublimation involve to create the lag?
+r_subl = h2subl/(1. - h2oice_porosity) &
+        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
+hsubl_dust = r_subl*str%dust_volfrac*str%thickness
+! Update of properties in the sublimating stratum
+str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
+str%elevation = str%elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
+str%co2ice_volfrac = h_co2ice_old/str%thickness
+str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
+str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
+str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
+! Correct the value of elevation for the upper strata
+current => str%up
+do while (associated(current))
+    current%elevation = current%down%elevation + current%thickness
+    current => current%up
+enddo
+! Remove the sublimating stratum if there is no ice anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+! New stratum = dust lag
+h_lag = hsubl_dust/(1. - dry_lag_porosity)
+if (h_lag > 0.) then ! Only if the dust lag is building up
+    if (new_lag) then
+        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
+        new_lag = .false.
+    else
+        str%up%thickness = str%up%thickness + h_lag
+        str%up%elevation = str%up%elevation + h_lag
+    endif
+endif
+
+END SUBROUTINE subl_h2oice_layering
+
+!=======================================================================
+! To lift dust in the layering
+SUBROUTINE lift_dust_lags(this,str,h2lift)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: str
+real, intent(in) :: h2lift
+
+!---- Code
+! Update of properties in the eroding dust lag
+str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
+str%elevation = str%elevation - h2lift/(1. - str%air_volfrac)
+! Remove the eroding dust lag if there is no dust anymore
+if (str%thickness < tol) call rm_stratum(this,str)
+
+END SUBROUTINE lift_dust_lags
+
+!=======================================================================
+! Layering algorithm
+SUBROUTINE make_layering(this,k,htend_co2ice,htend_h2oice,htend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)
+
+implicit none
+
+!---- Arguments
+type(layering),         intent(inout) :: this
+type(stratum), pointer, intent(inout) :: current1, current2
+logical,                intent(inout) :: new_str, new_lag1, new_lag2, stopPEM
+real, intent(in) :: htend_co2ice, htend_h2oice, htend_dust
+integer,  intent(in) :: k
+
+!---- Local variables
+real :: h_co2ice_old, h_h2oice_old, h_dust_old, thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
+
+!---- Code
+if (htend_dust < 0.) then ! Dust lifting
+    if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
+    write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting'
+    h2lift_tot = abs(htend_dust)
+    do while (h2lift_tot > 0. .and. associated(current1))
+        ! Is the considered stratum a dust lag?
+        if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we can lift dust
+            h_dust_old = current1%dust_volfrac*current1%thickness
+            ! How much dust can be lifted in the considered dust lag?
+            if (h2lift_tot - h_dust_old <= eps) then ! There is enough to lift
+                h2lift = h2lift_tot
+                h2lift_tot = 0.
+                call lift_dust_lags(this,current1,h2lift)
+            else ! Only a fraction can be lifted and so we move to the underlying stratum
+                h2lift = h_dust_old
+                h2lift_tot = h2lift_tot - h2lift
+                call lift_dust_lags(this,current1,h2lift)
+                current1 => current1%down
+            endif
+        else ! No, we stop
+            write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting: no dust lag!'
+            stopPEM = .true.
+            exit
+        endif
+    enddo
+    if (h2lift_tot > 0.) then
+        write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough dust in the lag layers to complete the lifting!'
+        stopPEM = .true.
+    endif
+
+else
+
+!------ Dust sedimentation
+    if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
+        write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust sedimentation'
+        ! New stratum at the layering top by sedimentation of dust
+        thickness = htend_dust/(1. - dry_regolith_porosity)
+        if (thickness > 0.) then ! Only if the stratum is builiding up
+             if (new_str) then
+                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
+                 new_str = .false.
+             else
+                 this%top%thickness = this%top%thickness + thickness
+                 this%top%elevation = this%top%elevation + thickness
+             endif
+        endif
+
+!------ Condensation of CO2 ice + H2O ice
+    else if ((htend_co2ice >= 0. .and. htend_h2oice > 0.) .or. (htend_co2ice > 0. .and. htend_h2oice >= 0.)) then
+        write(*,'(a5,i3,a)') 'Year ', k, ' -> CO2 and H2O ice condensation'
+        ! New stratum at the layering top by condensation of CO2 and H2O ice
+        thickness = htend_co2ice/(1. - co2ice_porosity) + htend_h2oice/(1. - h2oice_porosity) + htend_dust
+        if (thickness > 0.) then ! Only if the stratum is builiding up
+             if (new_str) then
+                 call add_stratum(this,thickness,this%top%elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
+                 new_str = .false.
+             else
+                 this%top%thickness = this%top%thickness + thickness
+                 this%top%elevation = this%top%elevation + thickness
+             endif
+        endif
+
+!------ Sublimation of CO2 ice + Condensation of H2O ice
+    else if (htend_co2ice < 0. .and. htend_h2oice >= 0. ) then
+        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 ice + Condensation of H2O ice'
+        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
+        h2subl_tot = abs(htend_co2ice)
+        do while (h2subl_tot > 0. .and. associated(current1))
+            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
+            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
+            h_dust_old = current1%dust_volfrac*current1%thickness
+            ! How much CO2 ice can sublimate in the considered stratum?
+            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
+                h2subl = h2subl_tot
+                h2subl_tot = 0.
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
+                current1 => current1%down
+                new_lag1 = .true.
+            else ! Only a fraction can sublimate and so we move to the underlying stratum
+                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)
+                current1 => current1%down
+                new_lag1 = .true.
+            endif
+        enddo
+        if (h2subl_tot > 0.) then
+            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
+            stopPEM = .true.
+        endif
+        ! New stratum at the layering top by condensation of H2O ice
+        thickness = htend_h2oice/(1. - h2oice_porosity) + htend_dust
+        if (thickness > 0.) then ! Only if the stratum is builiding up
+             if (new_str) then
+                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
+                 new_str = .false.
+             else
+                 this%top%thickness = this%top%thickness + thickness
+                 this%top%elevation = this%top%elevation + thickness
+             endif
+        endif
+
+!------ Sublimation of CO2 ice + H2O ice
+    else if (htend_co2ice < 0. .and. htend_h2oice < 0.) then
+        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 and H2O ice'
+        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
+        h2subl_tot = abs(htend_co2ice)
+        do while (h2subl_tot > 0. .and. associated(current1))
+            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
+            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
+            h_dust_old = current1%dust_volfrac*current1%thickness
+            ! How much CO2 ice can sublimate in the considered stratum?
+            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
+                h2subl = h2subl_tot
+                h2subl_tot = 0.
+                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
+            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
+                current1 => current1%down
+                new_lag1 = .true.
+            else ! Only a fraction can sublimate and so we move to the underlying stratum
+                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)
+                current1 => current1%down
+                new_lag1 = .true.
+            endif
+        enddo
+        if (h2subl_tot > 0.) then
+            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
+            stopPEM = .true.
+        endif
+        ! H2O ice sublimation in the considered stratum + New stratum for dust lag
+        h2subl_tot = abs(htend_h2oice)
+        do while (h2subl_tot > 0. .and. associated(current2))
+            h_co2ice_old = current2%co2ice_volfrac*current2%thickness
+            h_h2oice_old = current2%h2oice_volfrac*current2%thickness
+            h_dust_old = current2%dust_volfrac*current2%thickness
+            ! How much H2O ice can sublimate in the considered stratum?
+            if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate
+                h2subl = h2subl_tot
+                h2subl_tot = 0.
+                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
+            else if (h_h2oice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
+                current2 => current2%down
+                new_lag2 = .true.
+            else ! Only a fraction can sublimate and so we move to the underlying stratum
+                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)
+                current2 => current2%down
+                new_lag2 = .true.
+            endif
+        enddo
+        if (h2subl_tot > 0.) then
+            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough H2O ice in the layering to complete the sublimation!'
+            stopPEM = .true.
+        endif
+
+!------ Condensation of CO2 ice + Sublimation of H2O ice
+    else if (htend_co2ice >= 0. .and. htend_h2oice < 0.) then
+        error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
+    endif
+endif
+
+END SUBROUTINE make_layering
+
+END MODULE layering_mod