source: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 @ 3578

MODULE layering_mod

!=======================================================================
! Purpose: Manage layered deposits in the PEM with a linked list data structure
!
! Author: JBC
! Date: 08/03/2024
!=======================================================================

use glaciers_mod, only: rho_co2ice, rho_h2oice

implicit none

!---- Declarations
! Numerical parameter
real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
integer         :: nb_str_max
logical         :: layering_algo

! Physical parameters
real, parameter :: d_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_dust = 2500.             ! Density of dust [kg.m-3]
real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]

! Lag layer parameters -> see Levrard et al. 2007
real, parameter :: hmin_lag = 0.5  ! Minimal height of the lag deposit to reduce the sblimation rate
real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate

! Stratum = node of the linked list
type :: stratum
    real                   :: thickness      ! Layer thickness [m]
    real                   :: top_elevation  ! Layer top_elevation (top height from the surface) [m]
    real                   :: co2ice_volfrac ! CO2 ice volumetric fraction
    real                   :: h2oice_volfrac ! H2O ice volumetric fraction
    real                   :: dust_volfrac   ! Dust volumetric fraction
    real                   :: air_volfrac    ! Air volumetric 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

! Array of pointers to "stratum"
type ptrarray
    type(stratum), pointer :: p
end type ptrarray

!=======================================================================
contains
! Procedures to manipulate the data structure:
!     > print_stratum
!     > add_stratum
!     > insert_stratum
!     > rm_stratum
!     > get_stratum
!     > ini_layering
!     > del_layering
!     > print_layering
!     > get_nb_str_max
!     > stratif2array
!     > array2stratif
!     > print_stratif
! Procedures for the algorithm to evolve the layering:
!     > make_layering
!     > subl_co2ice_layering
!     > subl_h2oice_layering
!     > lift_dust_lags
!=======================================================================
! 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)') 'Top elevation  : ', this%top_elevation
write(*,'(a20,es13.6)') 'CO2 ice (m3/m3): ', this%co2ice_volfrac
write(*,'(a20,es13.6)') 'H2O ice (m3/m3): ', this%h2oice_volfrac
write(*,'(a20,es13.6)') 'Dust (m3/m3)   : ', this%dust_volfrac
write(*,'(a20,es13.6)') 'Air (m3/m3)    : ', this%air_volfrac

END SUBROUTINE print_stratum

!=======================================================================
! To add a stratum to the top of a layering
SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,air)

implicit none

!---- Arguments
type(layering), intent(inout)  :: this
real, intent(in), optional :: thickness, top_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%top_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(top_elevation)) str%top_elevation = top_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 volumetric 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,top_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, top_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,top_elevation,co2ice,h2oice,dust,air)
else
    ! Creation of the stratum
    allocate(str)
    nullify(str%up,str%down)
    str%thickness = 1.
    str%top_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(top_elevation)) str%top_elevation = top_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 volumetric 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 top_elevation for the upper strata
    current => str%up
    do while (associated(current))
        current%top_elevation = current%down%top_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%top_elevation <= 0.) return

! Decrement the number of strata
this%nb_str = this%nb_str - 1

! Remove the stratum from the layering
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 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 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 display a layering
SUBROUTINE print_layering(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_layering

!=======================================================================
! To get the maximum number of "stratum" in the stratification (layerings)
FUNCTION get_nb_str_max(stratif,ngrid,nslope) RESULT(nb_str_max)

implicit none

!---- Arguments
type(layering), dimension(ngrid,nslope), intent(in) :: stratif
integer,                                 intent(in) :: ngrid, nslope
integer :: nb_str_max

!---- Local variables
integer :: ig, islope

!---- Code
nb_str_max = 0
do islope = 1,nslope
    do ig = 1,ngrid
        nb_str_max = max(stratif(ig,islope)%nb_str,nb_str_max)
    enddo
enddo

END FUNCTION get_nb_str_max

!=======================================================================
! To convert the stratification data structure (layerings) into an array able to be outputted
SUBROUTINE stratif2array(stratif,ngrid,nslope,stratif_array)

implicit none

!---- Arguments
integer,                                 intent(in) :: ngrid, nslope
type(layering), dimension(ngrid,nslope), intent(in) :: stratif
real, dimension(:,:,:,:), allocatable, intent(inout) :: stratif_array

!---- Local variables
type(stratum), pointer :: current
integer                :: ig, islope, k

!---- Code
stratif_array = 0.
do islope = 1,nslope
    do ig = 1,ngrid
        current => stratif(ig,islope)%bottom
        k = 1
        do while (associated(current))
            stratif_array(ig,islope,k,1) = current%thickness
            stratif_array(ig,islope,k,2) = current%top_elevation
            stratif_array(ig,islope,k,3) = current%co2ice_volfrac
            stratif_array(ig,islope,k,4) = current%h2oice_volfrac
            stratif_array(ig,islope,k,5) = current%dust_volfrac
            stratif_array(ig,islope,k,6) = current%air_volfrac
            current => current%up
            k = k + 1
        enddo
    enddo
enddo

END SUBROUTINE stratif2array

!=======================================================================
! To convert the stratification array into the stratification data structure (layerings)
SUBROUTINE array2stratif(stratif_array,ngrid,nslope,stratif)

implicit none

!---- Arguments
integer,                               intent(in) :: ngrid, nslope
real, dimension(:,:,:,:), allocatable, intent(in) :: stratif_array
type(layering), dimension(ngrid,nslope), intent(inout) :: stratif

!---- Local variables
integer :: ig, islope, k

!---- Code
do islope = 1,nslope
    do ig = 1,ngrid
        stratif(ig,islope)%bottom%thickness      = stratif_array(ig,islope,1,1)
        stratif(ig,islope)%bottom%top_elevation  = stratif_array(ig,islope,1,2)
        stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3)
        stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4)
        stratif(ig,islope)%bottom%dust_volfrac   = stratif_array(ig,islope,1,5)
        stratif(ig,islope)%bottom%air_volfrac    = stratif_array(ig,islope,1,6)
        do k = 2,size(stratif_array,3)
            if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit
            call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6))
        enddo
    enddo
enddo

END SUBROUTINE array2stratif

!=======================================================================
! To display a stratification (layerings)
SUBROUTINE print_stratif(this,ngrid,nslope)

implicit none

!---- Arguments
integer,                                 intent(in) :: ngrid, nslope
type(layering), dimension(ngrid,nslope), intent(in) :: this

!---- Local variables
integer :: ig, islope

!---- Code
do ig = 1,ngrid
    write(*,'(a10,i4)') 'Grid cell ', ig
    do islope = 1,nslope
        write(*,'(a13,i1)') 'Slope number ', islope
        call print_layering(this(ig,islope))
        write(*,*) '~~~~~~~~~~'
    enddo
    write(*,*) '--------------------'
enddo

END SUBROUTINE print_stratif

!=======================================================================
! Layering algorithm
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(out) :: zshift_surf, zlag

!---- Local variables
real                   :: dh_co2ice, dh_h2oice, dh_dust
real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
real                   :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot
real                   :: h_toplag
type(stratum), pointer :: currentloc

!---- Code
dh_co2ice = d_co2ice/rho_co2ice
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
    if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
    write(*,'(a)') ' Stratification -> Dust lifting'
    h2lift_tot = abs(dh_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(*,'(a)') ' Stratification -> Dust lifting: no dust lag!'
            stopPEM = 8
            exit
        endif
    enddo
    if (h2lift_tot > 0.) then
        write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!'
        stopPEM = 8
    endif
    zshift_surf = dh_dust + h2lift_tot

else

!------ Dust sedimentation only
    if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
        write(*,'(a)') ' Stratification -> Dust sedimentation'
        ! New stratum at the layering top by sedimentation of dust
        thickness = dh_dust/(1. - dry_regolith_porosity)
        if (thickness > 0.) then ! Only if the stratum is building up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity)
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%top_elevation = this%top%top_elevation + thickness
             endif
        endif
        zshift_surf = thickness

!------ Condensation of CO2 ice + H2O ice
    else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then
        write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation'
        ! New stratum at the layering top by condensation of CO2 and H2O ice
        thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust
        if (thickness > 0.) then ! Only if the stratum is building up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness))
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%top_elevation = this%top%top_elevation + thickness
             endif
        endif
        zshift_surf = thickness

!------ Sublimation of CO2 ice + Condensation of H2O ice
    else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
        write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(dh_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 is CO2 ice sublimation inhibited by the top dust lag?
            h_toplag = 0.
            if (associated(current1%up)) then ! If there is an upper stratum
                currentloc => current1%up
                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
                    h_toplag = h_toplag + currentloc%thickness
                    if (.not. associated(currentloc%up)) exit
                    currentloc => currentloc%up
                enddo
            endif
            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)

            ! 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,zshift_surf,new_lag1,zlag)
            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                ! Is the considered stratum a dust lag?
                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                    current1 => current1%down
                    new_lag1 = .true.
                else ! No, so we stop here
                    exit
                endif
            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,zshift_surf,new_lag1,zlag)
                current1 => current1%down
                new_lag1 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
            stopPEM = 8
        endif
        ! New stratum at the layering top by condensation of H2O ice
        thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust
        if (thickness > 0.) then ! Only if the stratum is building up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness))
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%top_elevation = this%top%top_elevation + thickness
             endif
        endif
        zshift_surf = zshift_surf + thickness

!------ Sublimation of CO2 ice + H2O ice
    else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then
        write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(dh_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 is CO2 ice sublimation inhibited by the top dust lag?
            h_toplag = 0.
            if (associated(current1%up)) then ! If there is an upper stratum
                currentloc => current1%up
                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
                    h_toplag = h_toplag + currentloc%thickness
                    if (.not. associated(currentloc%up)) exit
                    currentloc => currentloc%up
                enddo
            endif
            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)

            ! 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,zshift_surf,new_lag1,zlag)
            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
                ! Is the considered stratum a dust lag?
                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                    current1 => current1%down
                    new_lag1 = .true.
                else ! No, so we stop here
                    exit
                endif
            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,zshift_surf,new_lag1,zlag)
                current1 => current1%down
                new_lag1 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!'
            stopPEM = 8
        endif
        ! H2O ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(dh_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 is H2O ice sublimation inhibited by the top dust lag?
            h_toplag = 0.
            if (associated(current2%up)) then ! If there is an upper stratum
                currentloc => current2%up
                do while (abs(1. - (currentloc%dust_volfrac + currentloc%air_volfrac)) < tol)
                    h_toplag = h_toplag + currentloc%thickness
                    if (.not. associated(currentloc%up)) exit
                    currentloc => currentloc%up
                enddo
            endif
            h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag)

            ! 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,zshift_surf,new_lag2,zlag)
            else if (h_h2oice_old < eps) then ! There is nothing to sublimate
                ! Is the considered stratum a dust lag?
                if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we move to the underlying stratum
                    current1 => current1%down
                    new_lag1 = .true.
                else ! No, so we stop here
                    exit
                endif
            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,zshift_surf,new_lag2,zlag)
                current2 => current2%down
                new_lag2 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!'
            stopPEM = 8
        endif

!------ Condensation of CO2 ice + Sublimation of H2O ice
    else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then
        error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
    endif
endif

END SUBROUTINE make_layering

!=======================================================================
! To sublimate CO2 ice in the layering
SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
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

!---- 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

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
    call rm_stratum(this,str)
else
    ! Update of properties in the sublimating stratum
    str%top_elevation = str%top_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 top_elevation for the upper strata
    current => str%up
    do while (associated(current))
        current%top_elevation = current%down%top_elevation + current%thickness
        current => current%up
    enddo
endif

! 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%top_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%top_elevation = str%up%top_elevation + h_lag
    endif
endif
zlag = zlag + h_lag

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,zshift_surf,new_lag,zlag)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
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

!---- 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

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
    call rm_stratum(this,str)
else
    ! Update of properties in the sublimating stratum
    str%top_elevation = str%top_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 top_elevation for the upper strata
    current => str%up
    do while (associated(current))
        current%top_elevation = current%down%top_elevation + current%thickness
        current => current%up
    enddo
endif

! 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%top_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%top_elevation = str%up%top_elevation + h_lag
    endif
endif
zlag = zlag + h_lag

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%top_elevation = str%top_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

END MODULE layering_mod