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

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.4      ! Porosity of dust lag
real, parameter :: dry_regolith_porosity = 0.4 ! Porosity of regolith
real, parameter :: breccia_porosity = 0.1      ! Porosity of breccia
real, parameter :: co2ice_porosity = 0.        ! Porosity of CO2 ice
real, parameter :: h2oice_porosity = 0.        ! 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 sublimation rate
real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate

! Stratum = node of the linked list
type :: stratum
    real                   :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
    real                   :: h_co2ice        ! Height of CO2 ice [m]
    real                   :: h_h2oice        ! Height of H2O ice [m]
    real                   :: h_dust          ! Height of dust [m]
    real                   :: h_pore          ! Height of pores [m]
    real                   :: poreice_volfrac ! Pore-filled ice volumetric fraction
    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_deposits
! Procedures to get information about a stratum:
!     > thickness
!     > is_co2ice_str
!     > is_h2oice_str
!     > is_dust_str
!     > is_dust_lag
!     > compute_h_toplag
! Procedures for the algorithm to evolve the layering:
!     > make_layering
!     > lift_dust_lag
!     > subl_ice_str
!=======================================================================
! To display a stratum
SUBROUTINE print_stratum(this)

implicit none

!---- Arguments
type(stratum), intent(in) :: this

!---- Code
write(*,'(a,es13.6)') 'Top elevation       [m]: ', this%top_elevation
write(*,'(a,es13.6)') 'Height of CO2 ice   [m]: ', this%h_co2ice
write(*,'(a,es13.6)') 'Height of H2O ice   [m]: ', this%h_h2oice
write(*,'(a,es13.6)') 'Height of dust      [m]: ', this%h_dust
write(*,'(a,es13.6)') 'Height of pores     [m]: ', this%h_pore
write(*,'(a,es13.6)') 'Pore-filled ice [m3/m3]: ', this%poreice_volfrac

END SUBROUTINE print_stratum

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

implicit none

!---- Arguments
type(layering), intent(inout)  :: this
real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice

!---- Local variables
type(stratum), pointer :: str

!---- Code
! Creation of the stratum
allocate(str)
nullify(str%up,str%down)
str%top_elevation = 1.
str%h_co2ice = 0.
str%h_h2oice = 0.
str%h_dust = 1.
str%h_pore = 0.
str%poreice_volfrac = 0.
if (present(top_elevation)) str%top_elevation = top_elevation
if (present(co2ice)) str%h_co2ice = co2ice
if (present(h2oice)) str%h_h2oice = h2oice
if (present(dust)) str%h_dust = dust
if (present(pore)) str%h_pore = pore
if (present(poreice)) str%poreice_volfrac = poreice

! 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,top_elevation,co2ice,h2oice,dust,pore,poreice)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str_prev
real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice

!---- 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,top_elevation,co2ice,h2oice,dust,pore,poreice)
else
    ! Creation of the stratum
    allocate(str)
    nullify(str%up,str%down)
    str%top_elevation = 1.
    str%h_co2ice = 0.
    str%h_h2oice = 0.
    str%h_dust = 1.
    str%h_pore = 0.
    str%poreice_volfrac = 0.
    if (present(top_elevation)) str%top_elevation = top_elevation
    if (present(co2ice)) str%h_co2ice = co2ice
    if (present(h2oice)) str%h_h2oice = h2oice
    if (present(dust)) str%h_dust = dust
    if (present(pore)) str%h_pore = pore
    if (present(poreice)) str%poreice_volfrac = poreice

    ! 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 + thickness(current)
        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 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)

use comsoil_h_PEM, only: soil_pem, nsoilmx_PEM, layer_PEM, index_breccia, index_bedrock

implicit none

!---- Arguments
type(layering), intent(inout) :: this

!---- Local variables
integer :: i
real    :: h_soil, h_pore, h_tot

!---- Code
! Creation of strata at the bottom of the layering to describe the sub-surface
if (soil_pem) then
    do i = nsoilmx_PEM,index_bedrock,-1
        h_soil = layer_PEM(i) - layer_PEM(i - 1) ! No porosity
        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,0.,0.)
    enddo
    do i = index_bedrock - 1,index_breccia,-1
        h_tot = layer_PEM(i) - layer_PEM(i - 1)
        h_pore = h_tot*breccia_porosity
        h_soil = h_tot - h_pore
        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
    enddo
    do i = index_breccia - 1,2,-1
        h_tot = layer_PEM(i) - layer_PEM(i - 1)
        h_pore = h_tot*dry_regolith_porosity
        h_soil = h_tot - h_pore
        call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.)
    enddo
    h_pore = layer_PEM(1)*dry_regolith_porosity
    h_soil = layer_PEM(1) - h_pore
    call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.)
endif

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(*,'(a,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%top_elevation
            stratif_array(ig,islope,k,2) = current%h_co2ice
            stratif_array(ig,islope,k,3) = current%h_h2oice
            stratif_array(ig,islope,k,4) = current%h_dust
            stratif_array(ig,islope,k,5) = current%h_pore
            stratif_array(ig,islope,k,6) = current%poreice_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
        do k = 1,size(stratif_array,3)
            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 deposits
SUBROUTINE print_deposits(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_deposits

!=======================================================================
!=======================================================================
! To compute the thickness of a stratum
FUNCTION thickness(this) RESULT(h_str)

implicit none

!---- Arguments
type(stratum), intent(in) :: this
real :: h_str

!---- Code
h_str = this%h_co2ice + this%h_h2oice + this%h_dust + this%h_pore

END FUNCTION thickness

!=======================================================================
! To know whether the stratum can be considered as a CO2 ice layer
FUNCTION is_co2ice_str(str) RESULT(co2ice_str)

implicit none

!---- Arguments
type(stratum), intent(in) :: str
logical :: co2ice_str

!---- Code
co2ice_str = .false.
if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.

END FUNCTION is_co2ice_str

!=======================================================================
! To know whether the stratum can be considered as a H2O ice layer
FUNCTION is_h2oice_str(str) RESULT(h2oice_str)

implicit none

!---- Arguments
type(stratum), intent(in) :: str
logical :: h2oice_str

!---- Code
h2oice_str = .false.
if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.

END FUNCTION is_h2oice_str

!=======================================================================
! To know whether the stratum can be considered as a dust layer
FUNCTION is_dust_str(str) RESULT(dust_str)

implicit none

!---- Arguments
type(stratum), intent(in) :: str
logical :: dust_str

!---- Code
dust_str = .false.
if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.

END FUNCTION is_dust_str

!=======================================================================
! To know whether the stratum can be considered as a dust lag
FUNCTION is_dust_lag(str) RESULT(dust_lag)

implicit none

!---- Arguments
type(stratum), intent(in) :: str
logical :: dust_lag

!---- Code
dust_lag = .false.
if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.

END FUNCTION is_dust_lag

!=======================================================================
! To get the top lag height
FUNCTION thickness_toplag(this) RESULT(h_toplag)

implicit none

!---- Arguments
type(layering), intent(in) :: this
real :: h_toplag

!---- Local variables
type(stratum), pointer :: current

!---- Code
h_toplag = 0.
current => this%top
! Is the considered stratum a dust lag?
do while (current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol .and. associated(current))
    h_toplag = h_toplag + thickness(current)
    current => current%down
enddo

END FUNCTION thickness_toplag

!=======================================================================
!=======================================================================
! Layering algorithm
SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)

implicit none

!---- Arguments
! Inputs

! Outputs
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: current
real,                   intent(inout) :: d_co2ice, d_h2oice
logical,                intent(inout) :: new_str, new_lag
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                   :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
real                   :: h2lift, h2lift_tot, h_pore, h_tot, 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 = this%top%top_elevation
zlag = 0.

! DUST SEDIMENTATION with possibly ice condensation
if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
    write(*,*) '> Layering: dust sedimentation'
    h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity)
    h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
    ! New stratum at the top of the layering
    if (new_str) then
        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.)
        new_str = .false.
    else
        this%top%top_elevation = this%top%top_elevation + h_tot
        this%top%h_co2ice = this%top%h_co2ice + dh_co2ice
        this%top%h_h2oice = this%top%h_h2oice + dh_h2oice
        this%top%h_dust = this%top%h_dust + dh_dust
        this%top%h_pore = this%top%h_pore + h_pore
    endif
!-----------------------------------------------------------------------
! DUST LIFTING
!(with possibly ice sublimation???)
!else if (dh_dust < dh_co2ice .and. dh_dust < dh_h2oice .and. dh_co2ice <= 0. .and. dh_h2oice <= 0.) then
else if (dh_dust < 0. .and. abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then
    write(*,*) '> Layering: dust lifting'
    h2lift_tot = abs(dh_dust)
    do while (h2lift_tot > 0. .and. associated(current))
        ! Is the considered stratum a dust lag?
        !if ((current%h_co2ice > eps .and. abs(dh_co2ice) < tol) .or. (current%h_h2oice > eps .and. abs(dh_h2oice) < tol)) then ! No, there is ice cementing the dust
        if (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) then ! No, there is ice cementing the dust
            dh_dust = 0. ! The tendency is set to 0
            exit
        else ! Yes, we can lift dust
            h2lift = min(h2lift_tot,current%h_dust) ! How much dust can be lifted?
            h2lift_tot = h2lift_tot - h2lift
            call lift_dust_lag(this,current,h2lift)
            if (h2lift_tot > 0.) current => current%down ! Still dust to be lifted so we move to the underlying stratum
        endif
    enddo
    if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
!-----------------------------------------------------------------------
! ICE CONDENSATION
else if ((dh_co2ice > dh_dust .or. dh_h2oice > dh_dust) .and. dh_dust >= 0.) then
    write(*,*) '> Layering: ice condensation'
    ! Which porosity is considered?
    if (dh_co2ice > dh_h2oice) then ! CO2 ice is dominant
        h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity)
    else ! H2O ice is dominant
        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
    endif
    h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore
    ! New stratum at the top of the layering
    if (new_str) then
        call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.)
        new_str = .false.
    else
        this%top%top_elevation = this%top%top_elevation + h_tot
        this%top%h_co2ice = this%top%h_co2ice + dh_co2ice
        this%top%h_h2oice = this%top%h_h2oice + dh_h2oice
        this%top%h_dust = this%top%h_dust + dh_dust
        this%top%h_pore = this%top%h_pore + h_pore
    endif
!-----------------------------------------------------------------------
! ICE SUBLIMATION
else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. abs(dh_dust) < tol) then
    write(*,*) '> Layering: ice sublimation'
    h2subl_co2ice_tot = abs(dh_co2ice)
    h2subl_h2oice_tot = abs(dh_h2oice)
    do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. associated(current))
        h_co2ice_old = current%h_co2ice
        h_h2oice_old = current%h_h2oice
        h_dust_old = current%h_dust

!~         ! How much is CO2 ice sublimation inhibited by the top dust lag?
!~         h_toplag = 0.
!~         if (associated(current%up)) then ! If there is an upper stratum
!~             currentloc => current%up
!~             ! Is the considered stratum a dust lag?
!~             do while (.not. (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) .and. associated(currentloc%up))
!~                 h_toplag = h_toplag + thickness(currentloc)
!~                 currentloc => currentloc%up
!~             enddo
!~         endif
!~         h2subl_co2ice_tot = h2subl_co2ice_tot*fred_subl**(h_toplag/hmin_lag)

        ! Is there CO2 ice to sublimate?
        h2subl_co2ice = 0.
        if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then
            h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
            h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
        endif
        ! Is there H2O ice to sublimate?
        h2subl_h2oice = 0.
        if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then
            h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old)
            h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice
        endif
        ! Sublimation
        if (h2subl_co2ice > 0. .or. h2subl_h2oice >0.) call subl_ice_str(this,current,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)
        ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum
        if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then 
            current => current%down
            new_lag = .true.
        else
            exit
        endif
    enddo
    if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0
    if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
!-----------------------------------------------------------------------
!~ ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
!~ else if (dh_co2ice < 0. .and. dh_h2oice > 0. .and. (abs(dh_dust) < abs(dh_co2ice) .or. abs(dh_dust) < abs(dh_h2oice))) then

! TO DO????

!-----------------------------------------------------------------------
! NOTHING HAPPENS
else if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
    ! We do nothing
!-----------------------------------------------------------------------
! UNUSUAL (WEIRD) SITUATION
else
    write(*,'(a)') 'Error: situation for the layering construction not managed!'
    write(*,'(a,es12.3)') '    > dh_co2ice [m/y] = ', dh_co2ice
    write(*,'(a,es12.3)') '    > dh_h2oice [m/y] = ', dh_h2oice
    write(*,'(a,es12.3)') '    > dh_dust   [m/y] = ', dh_dust, tol
!~     error stop
endif

zshift_surf = this%top%top_elevation - zshift_surf

END SUBROUTINE make_layering

!=======================================================================
! To lift dust in a stratum
SUBROUTINE lift_dust_lag(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%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
str%h_pore = str%h_pore - h2lift*str%h_pore/str%h_dust
str%h_dust = str%h_dust - h2lift

! Remove the eroding dust lag if there is no dust anymore
if (str%h_dust < tol) call rm_stratum(this,str)

END SUBROUTINE lift_dust_lag

!=======================================================================
! To sublimate ice in a stratum
SUBROUTINE subl_ice_str(this,str,h2subl_co2ice,h2subl_h2oice,zlag,new_lag)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str
logical,                intent(inout) :: new_lag
real,                   intent(inout) :: zlag
real, intent(in) :: h2subl_co2ice, h2subl_h2oice

!---- Local variables
real                   :: hsubl_dust, h_pore, h_ice, h_tot, h2subl
type(stratum), pointer :: current

!---- Code
h_ice = str%h_co2ice + str%h_h2oice
h2subl = h2subl_co2ice + h2subl_h2oice

! How much dust does ice sublimation involve to create the lag?
hsubl_dust = str%h_dust*h2subl/h_ice

! Update of properties in the sublimating stratum
str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust
str%h_co2ice = str%h_co2ice - h2subl_co2ice
str%h_h2oice = str%h_h2oice - h2subl_h2oice
str%h_dust = str%h_dust - hsubl_dust
str%h_pore = str%h_pore - h2subl*h_pore/h_ice

! Correct the value of top_elevation for the upper strata
current => str%up
do while (associated(current))
    current%top_elevation = current%down%top_elevation + thickness(current)
    current => current%up
enddo

! Remove the sublimating stratum if there is no ice anymore
if (str%h_co2ice < tol .and. str%h_h2oice < tol) call rm_stratum(this,str)

! New stratum above the current stratum as a dust lag
if (hsubl_dust > 0.) then ! Only if the dust lag is building up
    h_pore = hsubl_dust*dry_lag_porosity/(1. - dry_lag_porosity)
    h_tot = hsubl_dust + h_pore
    if (new_lag) then
        call insert_stratum(this,str,str%top_elevation + h_tot,0.,0.,hsubl_dust,h_pore)
        new_lag = .false.
    else
        str%up%top_elevation = str%up%top_elevation + h_tot
    endif
endif
zlag = zlag + h_tot

END SUBROUTINE subl_ice_str

END MODULE layering_mod