MODULE layering_mod

!=======================================================================
! Purpose: Manage layered layerings_map 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
logical         :: impose_dust_ratio = .false. ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
real            :: dust2ice_ratio    = 0.01    ! Dust-to-ice ratio when ice condenses (1% by default)
real            :: d_dust            = 1.e-3   ! Tendency of dust [kg.m-2.y-1]
real, parameter :: dry_lag_porosity  = 0.4     ! Porosity of dust lag
real, parameter :: wet_lag_porosity  = 0.1     ! Porosity of water ice lag
real, parameter :: 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]
real, parameter :: h_patchy_dust     = 5.e-4   ! Height under which a dust lag is considered patchy [m]
real, parameter :: h_patchy_ice      = 5.e-4   ! Height under which a H2O ice lag is considered patchy [m]

! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
real, parameter :: hmin_lag_co2 = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate [m]
real, parameter :: fred_sublco2 = 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)
    !real                   :: h_toplag         ! Height of dust layer at the top [m]
    !type(stratum), pointer :: subice => null() ! Shallower stratum containing ice
    !type(stratum), pointer :: surf => null()   ! Surface stratum
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
!     > del_stratum
!     > get_stratum
!     > ini_layering
!     > del_layering
!     > print_layering
!     > get_nb_str_max
!     > stratif2array
!     > array2stratif
!     > print_layerings_map
! Procedures to get information about a stratum:
!     > thickness
!     > is_co2ice_str
!     > is_h2oice_str
!     > is_dust_str
!     > is_dust_lag
!     > subsurface_ice_layering
! Procedures for the algorithm to evolve the layering:
!     > make_layering
!     > lift_dust
!     > 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 del_stratum(this,str)

implicit none

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

!---- Local variables
type(stratum), pointer :: new_str ! To make the argument point towards something

!---- Code
! Protect the sub-surface strata from deletion
if (str%top_elevation <= 0.) return

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

! Making the new links
if (associated(str%down)) then ! If there is a lower stratum
    str%down%up => str%up
else ! If it was the first stratum (bottom of layering)
    this%bottom => str%up
endif
if (associated(str%up)) then ! If there is an upper stratum
    str%up%down => str%down
else ! If it was the last stratum (top of layering)
    this%top => str%down
endif

! Remove the stratum from the layering
new_str => str%down
nullify(str%up,str%down)
deallocate(str)
nullify(str)
str => new_str

END SUBROUTINE del_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*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)*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(layerings_map,ngrid,nslope) RESULT(nb_str_max)

implicit none

!---- Arguments
type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
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(layerings_map(ig,islope)%nb_str,nb_str_max)
    enddo
enddo

END FUNCTION get_nb_str_max

!=======================================================================
! To convert the layerings map into an array able to be outputted
SUBROUTINE stratif2array(layerings_map,ngrid,nslope,stratif_array)

implicit none

!---- Arguments
integer,                                 intent(in) :: ngrid, nslope
type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
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 => layerings_map(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 layerings map
SUBROUTINE array2stratif(stratif_array,ngrid,nslope,layerings_map)

implicit none

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

!---- 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(layerings_map(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 the map of layerings
SUBROUTINE print_layerings_map(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_layerings_map

!=======================================================================
!=======================================================================
! 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 data about possible subsurface ice in a layering
SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice)

implicit none

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

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

!---- Code
h_toplag = 0.
h2o_ice = 0.
co2_ice = 0.
current => this%top
! Is the considered stratum a dust lag?
if (.not. is_dust_lag(current)) return
do while (is_dust_lag(current))
    h_toplag = h_toplag + thickness(current)
    if (associated(current%down)) then
        current => current%down
    else
        exit
    endif
enddo

if (.not. associated(current%down)) then ! Bottom of the layering -> no ice below
    h_toplag = 0.
    return
endif

! Is the stratum below the dust lag made of ice?
if (current%h_h2oice > 0. .and. current%h_h2oice > current%h_co2ice) then ! Yes, there is more H2O ice than CO2 ice
    if (h_toplag < h_patchy_dust) then ! But the dust lag is too thin to considered ice as subsurface ice
        h_toplag = 0.
        h2o_ice = current%h_h2oice
    endif
else if (current%h_co2ice > 0. .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
    h_toplag = 0. ! Because it matters only for H2O ice
    if (h_toplag < h_patchy_ice) then ! But the dust lag is too thin to considered ice as subsurface ice
        co2_ice = current%h_co2ice
    endif
endif

END SUBROUTINE subsurface_ice_layering

!=======================================================================
!=======================================================================
! To lift dust in a stratum
SUBROUTINE lift_dust(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 del_stratum(this,str)

END SUBROUTINE lift_dust

!=======================================================================
! To sublimate ice in a stratum
SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,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) :: zlag
real, intent(in) :: h2subl_co2ice, h2subl_h2oice

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

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

! How much dust and H2O ice does ice sublimation involve to create the lag?
hlag_dust = str%h_dust*h2subl/h_ice
hlag_h2oice = 0.
if (h2subl_co2ice > 0. .and. h2subl_h2oice < tol) hlag_h2oice = str%h_h2oice*h2subl/h_ice

! Update of properties in the sublimating stratum
str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hlag_dust - hlag_h2oice
str%h_co2ice      = str%h_co2ice      - h2subl_co2ice
str%h_h2oice      = str%h_h2oice      - h2subl_h2oice - hlag_h2oice
str%h_dust        = str%h_dust        - hlag_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 del_stratum(this,str)

! Which porosity is considered?
if (hlag_dust >= hlag_h2oice) then
    h_pore_new = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity)
else
    h_pore_new = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity)
endif
h_tot = hlag_dust + hlag_h2oice + h_pore_new

! New stratum above the current stratum as a lag layer
if (new_lag) then
    call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore_new,0.)
    new_lag = .false.
else
    str%up%top_elevation = str%up%top_elevation + h_tot
    str%up%h_h2oice      = str%up%h_h2oice      + hlag_h2oice
    str%up%h_dust        = str%up%h_dust        + hlag_dust
    str%up%h_pore        = str%up%h_pore        + h_pore_new
endif

zlag = zlag + h_tot

END SUBROUTINE sublimate_ice

!=======================================================================
! 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, R_subl
logical                :: unexpected
type(stratum), pointer :: currentloc

!---- Code
dh_co2ice = d_co2ice/rho_co2ice
dh_h2oice = d_h2oice/rho_h2oice
! To disable dust sedimenation when there is ice sublimation (because dust tendency is fixed)
if (dh_co2ice >= 0. .and. dh_h2oice >= 0.) then
    if (impose_dust_ratio) then
        dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
    else 
        dh_dust = d_dust/rho_dust
    endif
else
    dh_dust = 0.
endif
zshift_surf = this%top%top_elevation
zlag = 0.
unexpected = .false.

!-----------------------------------------------------------------------
! NOTHING HAPPENS
if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
    ! So we do nothing
!-----------------------------------------------------------------------
! BUILDING A STRATUM (ice condensation and/or dust desimentation)
else if (dh_co2ice >= 0. .and. dh_h2oice >= 0. .and. dh_dust >= 0.) then
    ! Which porosity is considered?
    if (dh_co2ice > dh_h2oice .and. dh_co2ice > dh_dust) then ! CO2 ice is dominant
        h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity)
    else if (dh_h2oice > dh_co2ice .and. dh_h2oice > dh_dust) then ! H2O ice is dominant
        h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity)
    else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant
        h_pore = dh_dust*regolith_porosity/(1. - regolith_porosity)
    else
        unexpected = .true.
    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
!-----------------------------------------------------------------------
 ! RETREATING A STRATUM (ice sublimation and/or dust lifting)
else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then
        h2subl_co2ice_tot = abs(dh_co2ice)
        h2subl_h2oice_tot = abs(dh_h2oice)
        h2lift_tot = abs(dh_dust)
        do while ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. associated(current))
            h_co2ice_old = current%h_co2ice
            h_h2oice_old = current%h_h2oice

            ! How much is CO2 ice sublimation inhibited by the top dust lag?
            R_subl = 1.
            if (associated(current%up)) then ! If there is an upper stratum
                h_toplag = 0.
                currentloc => current%up
                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
                    h_toplag = h_toplag + thickness(currentloc)
                    currentloc => currentloc%up
                enddo
                if (currentloc%h_h2oice >= h_patchy_ice) then
                    R_subl = 0.
                else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then
                    h_toplag = h_toplag + thickness(currentloc)
                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
                else
                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
                endif
            endif
           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl

            ! 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
            ! Ice sublimation
            if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)

            ! Is there dust to lift?
            h2lift = 0.
            if (is_dust_lag(current) .and. h2lift_tot > 0.) then
                h2lift = min(h2lift_tot,current%h_dust)
                h2lift_tot = h2lift_tot - h2lift
                ! Dust lifting
                if (h2lift > 0.) call lift_dust(this,current,h2lift)
            else if (associated(current%up)) then
                if (is_dust_lag(current%up) .and. h2lift_tot > 0.) then
                    h2lift = min(h2lift_tot,current%up%h_dust)
                    h2lift_tot = h2lift_tot - h2lift
                    ! Dust lifting
                    if (h2lift > 0.) call lift_dust(this,current%up,h2lift)
                endif
            endif

            ! Still ice to be sublimated or dust to be lifted 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. .or. h2lift_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < 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
        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
!-----------------------------------------------------------------------
! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then
!~     if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then ! CO2 sublimation is dominant

!~     else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then ! H2O condensation is dominant

!~     else
        unexpected = .true.
!~     endif
!-----------------------------------------------------------------------
! NOT EXPECTED SITUATION
else
    unexpected = .true.
endif

if (unexpected) then
    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
    error stop
endif

zshift_surf = this%top%top_elevation - zshift_surf

END SUBROUTINE make_layering

END MODULE layering_mod