source: trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 @ 4074

MODULE layered_deposits
!-----------------------------------------------------------------------
! NAME
!     layered_deposits
!
! DESCRIPTION
!     Manage layered deposits in the PEM with a linked list data structure.
!     Provides data types and subroutines to manipulate layers of ice and dust.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4, tol
use geometry, only: ngrid, nslope
use surf_ice, only: rho_co2ice, rho_h2oice

! DECLARATION
! -----------
implicit none

! VARIABLES
! ---------
! Numerical parameter
logical(k4),  protected :: do_layering

! Physical parameters
logical(k4), protected :: impose_dust_ratio            ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio')
real(dp),    protected :: dust2ice_ratio               ! Dust-to-ice ratio when ice condenses (10% by default)
real(dp),    protected :: d_dust                       ! Tendency of dust [kg.m-2.y-1]
real(dp),    parameter :: dry_lag_porosity  = 0.2_dp   ! Porosity of dust lag
real(dp),    parameter :: wet_lag_porosity  = 0.15_dp  ! Porosity of water ice lag
real(dp),    parameter :: regolith_porosity = 0.4_dp   ! Porosity of regolith
real(dp),    parameter :: breccia_porosity  = 0.1_dp   ! Porosity of breccia
real(dp),    parameter :: co2ice_porosity   = 0.05_dp  ! Porosity of CO2 ice
real(dp),    parameter :: h2oice_porosity   = 0.1_dp   ! Porosity of H2O ice
real(dp),    parameter :: rho_dust          = 2500._dp ! Density of dust [kg.m-3]
real(dp),    parameter :: rho_rock          = 3200._dp ! Density of rock [kg.m-3]
real(dp),    parameter :: h_patchy_dust     = 5.e-4_dp ! Height under which a dust lag is considered patchy [m]
real(dp),    parameter :: h_patchy_ice      = 5.e-4_dp ! Height under which a H2O ice lag is considered patchy [m]

! Parameters for CO2 flux correction (inspired by Levrard et al. 2007)
real(dp), parameter :: hmin_lag_co2 = 0.5_dp ! Minimal height of the lag deposit to reduce the sublimation rate [m]
real(dp), parameter :: fred_sublco2 = 0.1_dp ! Reduction factor of sublimation rate

! DATA STRUCTURES
! ---------------
! Stratum = node of the linked list
type :: stratum
    real(dp)               :: top_elevation   ! Layer top_elevation (top height from the surface) [m]
    real(dp)               :: h_co2ice        ! Height of CO2 ice [m]
    real(dp)               :: h_h2oice        ! Height of H2O ice [m]
    real(dp)               :: h_dust          ! Height of dust [m]
    real(dp)               :: h_pore          ! Height of pores [m]
    real(dp)               :: 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(di)            :: 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(dp)               :: 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
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE set_layered_deposits_config(do_layering_in,impose_dust_ratio_in,d_dust_in,dust2ice_ratio_in)
!-----------------------------------------------------------------------
! NAME
!     set_layered_deposits_config
!
! DESCRIPTION
!     Setter for 'layered_deposits' configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: real2str, bool2str
use display, only: print_msg

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
logical(k4), intent(in) :: do_layering_in, impose_dust_ratio_in
real(dp),    intent(in) :: d_dust_in, dust2ice_ratio_in

! CODE
! ----
do_layering = do_layering_in
impose_dust_ratio = impose_dust_ratio_in
d_dust = d_dust_in
dust2ice_ratio = dust2ice_ratio_in
call print_msg('do_layering       = '//bool2str(do_layering))
call print_msg('impose_dust_ratio = '//bool2str(impose_dust_ratio))
call print_msg('d_dust            = '//real2str(d_dust))
call print_msg('dust2ice_ratio    = '//real2str(dust2ice_ratio))

END SUBROUTINE set_layered_deposits_config
!=======================================================================

!=======================================================================
SUBROUTINE print_stratum(this)
!-----------------------------------------------------------------------
! NAME
!     print_stratum
!
! DESCRIPTION
!     Print a stratum.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: real2str
use display, only: print_msg

! DECLARATION
! -----------
implicit none

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

! CODE
! ----
call print_msg('Top elevation       [m]: '//real2str(this%top_elevation))
call print_msg('Height of CO2 ice   [m]: '//real2str(this%h_co2ice))
call print_msg('Height of H2O ice   [m]: '//real2str(this%h_h2oice))
call print_msg('Height of dust      [m]: '//real2str(this%h_dust))
call print_msg('Height of pores     [m]: '//real2str(this%h_pore))
call print_msg('Pore-filled ice [m3/m3]: '//real2str(this%poreice_volfrac))

END SUBROUTINE print_stratum
!=======================================================================

!=======================================================================
SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
!-----------------------------------------------------------------------
! NAME
!     add_stratum
!
! DESCRIPTION
!     Add a stratum to the top of a layering.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), intent(inout)        :: this
real(dp),       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._dp
str%h_co2ice        = 0._dp
str%h_h2oice        = 0._dp
str%h_dust          = 1._dp
str%h_pore          = 0._dp
str%poreice_volfrac = 0._dp
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
end if

END SUBROUTINE add_stratum
!=======================================================================

!=======================================================================
SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice)
!-----------------------------------------------------------------------
! NAME
!     insert_stratum
!
! DESCRIPTION
!     Insert a stratum after a given previous stratum in the layering.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering),         intent(inout)        ::  this
type(stratum), pointer, intent(inout)        ::  str_prev
real(dp),               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._dp
    str%h_co2ice        = 0._dp
    str%h_h2oice        = 0._dp
    str%h_dust          = 1._dp
    str%h_pore          = 0._dp
    str%poreice_volfrac = 0._dp
    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
    end do
end if

END SUBROUTINE insert_stratum
!=======================================================================

!=======================================================================
SUBROUTINE del_stratum(this,str)
!-----------------------------------------------------------------------
! NAME
!     del_stratum
!
! DESCRIPTION
!     Delete a stratum from the layering and fix links.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
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 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
end if
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
end if

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

!=======================================================================
FUNCTION get_stratum(lay,i) RESULT(str)
!-----------------------------------------------------------------------
! NAME
!     get_stratum
!
! DESCRIPTION
!     Return the i-th stratum in the layering.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), intent(in) :: lay
integer(di),    intent(in) :: i
type(stratum), pointer     :: str

! LOCAL VARIABLES
! ---------------
integer(di) :: k

! CODE
! ----
if (i < 1 .or. i > lay%nb_str) call stop_clean(__FILE__,__LINE__,'requested number is beyond the number of strata of the layering!',1)
k = 1
str => lay%bottom
do while (k /= i)
    str => str%up
    k = k + 1
end do

END FUNCTION get_stratum
!=======================================================================

!=======================================================================
SUBROUTINE ini_layering(this)
!-----------------------------------------------------------------------
! NAME
!     del_layering
!
! DESCRIPTION
!     Initialize the layering base structure.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,     only: do_soil, layer, index_breccia, index_bedrock
use geometry, only: nsoil

! DECLARATION
! -----------
implicit none

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

! LOCAL VARIABLES
! ---------------
integer(di) :: i
real(dp)    :: h_soil, h_pore, h_tot

! CODE
! ----
! Creation of strata at the bottom of the layering to describe the sub-surface
if (do_soil) then
    do i = nsoil,index_bedrock,-1
        h_soil = layer(i) - layer(i - 1) ! No porosity
        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,0._dp,0._dp)
    end do
    do i = index_bedrock - 1,index_breccia,-1
        h_tot = layer(i) - layer(i - 1)
        h_pore = h_tot*breccia_porosity
        h_soil = h_tot - h_pore
        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,0._dp)
    end do
    do i = index_breccia - 1,2,-1
        h_tot = layer(i) - layer(i - 1)
        h_pore = h_tot*regolith_porosity
        h_soil = h_tot - h_pore
        call add_stratum(this,-layer(i - 1),0._dp,0._dp,h_soil,h_pore,0._dp)
    end do
    h_pore = layer(1)*regolith_porosity
    h_soil = layer(1) - h_pore
    call add_stratum(this,0._dp,0._dp,0._dp,h_soil,h_pore,0._dp)
end if

END SUBROUTINE ini_layering
!=======================================================================

!=======================================================================
SUBROUTINE del_layering(this)
!-----------------------------------------------------------------------
! NAME
!     del_layering
!
! DESCRIPTION
!     Delete all strata in a layering and reset counters.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
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
end do
this%nb_str = 0

END SUBROUTINE del_layering
!=======================================================================

!=======================================================================
SUBROUTINE print_layering(this)
!-----------------------------------------------------------------------
! NAME
!     print_layering
!
! DESCRIPTION
!     Display all strata from top to bottom for a layering.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: int2str
use display, only: print_msg

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), intent(in) :: this

! LOCAL VARIABLES
! ---------------
type(stratum), pointer :: current
integer(di)            :: i

! CODE
! ----
current => this%top

i = this%nb_str
do while (associated(current))
    call print_msg('Stratum '//int2str(i))
    call print_stratum(current)
    call print_msg('')
    current => current%down
    i = i - 1
end do

END SUBROUTINE print_layering
!=======================================================================

!=======================================================================
FUNCTION get_nb_str_max(layerings_map) RESULT(nb_str_max)
!-----------------------------------------------------------------------
! NAME
!     get_nb_str_max
!
! DESCRIPTION
!     Return the maximum number of strata across all grid/slope layerings.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), dimension(:,:), intent(in) :: layerings_map
integer(di)                                :: nb_str_max

! LOCAL VARIABLES
! ---------------
integer(di) :: 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)
    end do
end do

END FUNCTION get_nb_str_max
!=======================================================================

!=======================================================================
SUBROUTINE map2array(layerings_map,layerings_array)
!-----------------------------------------------------------------------
! NAME
!     map2array
!
! DESCRIPTION
!     Convert the layerings map into an allocatable array for output.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), dimension(:,:),     intent(in)  :: layerings_map
real(dp),       dimension(:,:,:,:), intent(out) :: layerings_array

! LOCAL VARIABLES
! ---------------
type(stratum), pointer :: current
integer(di)            :: ig, islope, k

! CODE
! ----
layerings_array = 0.
do islope = 1,nslope
    do ig = 1,ngrid
        current => layerings_map(ig,islope)%bottom
        k = 1
        do while (associated(current))
            layerings_array(ig,islope,k,1) = current%top_elevation
            layerings_array(ig,islope,k,2) = current%h_co2ice
            layerings_array(ig,islope,k,3) = current%h_h2oice
            layerings_array(ig,islope,k,4) = current%h_dust
            layerings_array(ig,islope,k,5) = current%h_pore
            layerings_array(ig,islope,k,6) = current%poreice_volfrac
            current => current%up
            k = k + 1
        end do
    end do
end do

END SUBROUTINE map2array
!=======================================================================

!=======================================================================
SUBROUTINE array2map(layerings_array,layerings_map)
!-----------------------------------------------------------------------
! NAME
!     array2map
!
! DESCRIPTION
!     Convert the stratification array back into the layerings map.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp),       dimension(:,:,:,:), intent(in)    :: layerings_array
type(layering), dimension(:,:),     intent(inout) :: layerings_map

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, islope, k

! CODE
! ----
do islope = 1,nslope
    do ig = 1,ngrid
        do k = 1,size(layerings_array,3)
            call add_stratum(layerings_map(ig,islope),layerings_array(ig,islope,k,1),layerings_array(ig,islope,k,2),layerings_array(ig,islope,k,3),layerings_array(ig,islope,k,4),layerings_array(ig,islope,k,5),layerings_array(ig,islope,k,6))
        end do
    end do
end do

END SUBROUTINE array2map
!=======================================================================

!=======================================================================
SUBROUTINE layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
!-----------------------------------------------------------------------
! NAME
!     layering2surfice
!
! DESCRIPTION
!     Convert the layerings map into surface ice arrays.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), dimension(:,:), intent(in)  :: layerings_map
real(dp),       dimension(:,:), intent(out) :: h2o_ice, co2_ice, h2oice_depth

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, islope

! CODE
! ----
do ig = 1,ngrid
    do islope = 1,nslope
        if (is_co2ice_str(layerings_map(ig,islope)%top)) then
            co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice
        else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
            h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice
        else
            call subsurface_ice_layering(layerings_map(ig,islope),h2oice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope))
        end if
    end do
end do

END SUBROUTINE layering2surfice
!=======================================================================

!=======================================================================
SUBROUTINE surfice2layering(h2o_ice,co2_ice,layerings_map)
!-----------------------------------------------------------------------
! NAME
!     surfice2layering
!
! DESCRIPTION
!     Convert surface ice arrays into the layerings map.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp),       dimension(:,:), intent(in)    :: h2o_ice, co2_ice
type(layering), dimension(:,:), intent(inout) :: layerings_map

! LOCAL VARIABLES
! ---------------
integer(di) :: i, islope

! CODE
! ----
do islope = 1,nslope
    do i = 1,ngrid
        layerings_map(i,islope)%top%h_h2oice = h2o_ice(i,islope)/rho_h2oice
        layerings_map(i,islope)%top%h_co2ice = co2_ice(i,islope)/rho_co2ice
    end do
end do

END SUBROUTINE surfice2layering
!=======================================================================

!=======================================================================
SUBROUTINE print_layerings_map(this)
!-----------------------------------------------------------------------
! NAME
!     print_layerings_map
!
! DESCRIPTION
!     Display the entire layerings map for all grids and slopes.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: int2str
use display, only: print_msg

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering), dimension(:,:), intent(in) :: this

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, islope

! CODE
! ----
do ig = 1,ngrid
    call print_msg('Grid cell '//int2str(ig))
    do islope = 1,nslope
        call print_msg('Slope number '//int2str(islope))
        call print_layering(this(ig,islope))
        call print_msg('~~~~~~~~~~')
    end do
    call print_msg('--------------------')
end do

END SUBROUTINE print_layerings_map
!=======================================================================

!=======================================================================
FUNCTION thickness(this) RESULT(h_str)
!-----------------------------------------------------------------------
! NAME
!     thickness
!
! DESCRIPTION
!     Compute the total thickness of a stratum.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

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

! LOCAL VARIABLES
! ---------------
real(dp) :: h_str

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

END FUNCTION thickness
!=======================================================================

!=======================================================================
FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
!-----------------------------------------------------------------------
! NAME
!     is_co2ice_str
!
! DESCRIPTION
!     Check if a stratum can be considered as a CO2 ice layer.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(stratum), intent(in) :: str

! LOCAL VARIABLES
! ---------------
logical(k4) :: 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
!=======================================================================

!=======================================================================
FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
!-----------------------------------------------------------------------
! NAME
!     is_h2oice_str
!
! DESCRIPTION
!     Check if a stratum can be considered as a H2O ice layer.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(stratum), intent(in) :: str

! LOCAL VARIABLES
! ---------------
logical(k4) :: 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
!=======================================================================

!=======================================================================
FUNCTION is_dust_str(str) RESULT(dust_str)
!-----------------------------------------------------------------------
! NAME
!     is_dust_str
!
! DESCRIPTION
!     Check if a stratum can be considered as a dust layer.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(stratum), intent(in) :: str

! LOCAL VARIABLES
! ---------------
logical(k4) :: 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
!=======================================================================

!=======================================================================
FUNCTION is_dust_lag(str) RESULT(dust_lag)
!-----------------------------------------------------------------------
! NAME
!     is_dust_lag
!
! DESCRIPTION
!     Check if a stratum is a dust lag (no ice content).
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(stratum), intent(in) :: str

! LOCAL VARIABLES
! ---------------
logical(k4) :: 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
!=======================================================================

!=======================================================================
SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice)
!-----------------------------------------------------------------------
! NAME
!     subsurface_ice_layering
!
! DESCRIPTION
!     Get data about possible subsurface ice in a layering.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

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

! LOCAL VARIABLES
! ---------------
type(stratum), pointer :: current

! CODE
! ----
h_toplag = 0._dp
h2o_ice = 0._dp
co2_ice = 0._dp
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
    end if
end do

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

! Is the stratum below the dust lag made of ice?
if (current%h_h2oice > 0._dp .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._dp
        h2o_ice = current%h_h2oice*rho_h2oice
    end if
else if (current%h_co2ice > 0._dp .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
    h_toplag = 0._dp ! Because it matters only for H2O ice
    if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice*rho_co2ice ! But the dust lag is too thin to considered ice as subsurface ice
end if

END SUBROUTINE subsurface_ice_layering
!=======================================================================

!=======================================================================
SUBROUTINE lift_dust(this,str,h2lift)
!-----------------------------------------------------------------------
! NAME
!     lift_dust
!
! DESCRIPTION
!     Lift dust in a stratum by wind erosion or other processes.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str
real(dp),               intent(in)    :: h2lift

! CODE
! ----
! Update of properties in the eroding dust lag
str%top_elevation = str%top_elevation - h2lift*(1._dp + 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
!=======================================================================

!=======================================================================
SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
!-----------------------------------------------------------------------
! NAME
!     sublimate_ice
!
! DESCRIPTION
!     Handle ice sublimation in a stratum and create lag layer.
!
! AUTHORS & DATE
!     JB Clement, 2024
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

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

! LOCAL VARIABLES
! ---------------
real(dp)               :: h2subl, h_ice, h_pore, h_pore_new, h_tot
real(dp)               :: 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._dp
if (h2subl_co2ice > 0._dp .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._dp + 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
end do

! 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._dp - dry_lag_porosity)
else
    h_pore_new = hlag_h2oice*wet_lag_porosity/(1._dp - wet_lag_porosity)
end if
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._dp,hlag_h2oice,hlag_dust,h_pore_new,0._dp)
    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
end if

zlag = zlag + h_tot

END SUBROUTINE sublimate_ice
!=======================================================================

!=======================================================================
SUBROUTINE evolve_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
!-----------------------------------------------------------------------
! NAME
!     evolve_layering
!
! DESCRIPTION
!     Main algorithm for layering evolution. Manages ice condensation/sublimation
!     and dust sedimentation/lifting. Handles multiple scenarios: building new
!     strata from ice/dust, retreating strata from sublimation/lifting, and mixed
!     scenarios with CO2 ice sublimation combined with H2O ice condensation.
!
! AUTHORS & DATE
!     JB Clement, 09/2024
!
! NOTES
!     Creates dust lag layers when ice sublimation occurs.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean
use utility,  only: real2str
use display,  only: print_err

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
type(layering),         intent(inout) ::  this
type(stratum), pointer, intent(inout) ::  current
real(dp),               intent(inout) ::  d_co2ice, d_h2oice
logical(k4),            intent(inout) ::  new_str, new_lag
real(dp),               intent(out)   ::  zshift_surf, zlag

! LOCAL VARIABLES
! ---------------
real(dp)               :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
real(dp)               :: h_co2ice_old, h_h2oice_old
real(dp)               :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
real(dp)               :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag, R_subl
logical(k4)            :: unexpected
type(stratum), pointer :: currentloc

! CODE
! ----
dh_co2ice = d_co2ice/rho_co2ice
dh_h2oice = d_h2oice/rho_h2oice
dh_dust = 0._dp
if (dh_h2oice > 0._dp) then ! To put a dust sedimentation tendency only when ice is condensing
    if (impose_dust_ratio) then
        if (dh_co2ice > 0._dp) then
            dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice)
        else
            dh_dust = dust2ice_ratio*dh_h2oice
        end if
    else
        dh_dust = d_dust/rho_dust
    end if
end if
zshift_surf = this%top%top_elevation
zlag = 0._dp
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._dp .and. dh_h2oice >= 0._dp .and. dh_dust >= 0._dp) 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._dp - 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._dp - h2oice_porosity)
    else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant
        h_pore = dh_dust*regolith_porosity/(1._dp - regolith_porosity)
    else
        unexpected = .true.
    end if
    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._dp)
        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
    end if
!~~~~~~~~~~
 ! RETREATING A STRATUM (ice sublimation and/or dust lifting)
else if (dh_co2ice <= 0._dp .and. dh_h2oice <= 0._dp .and. dh_dust <= 0._dp) then
        h2subl_co2ice_tot = abs(dh_co2ice)
        h2subl_h2oice_tot = abs(dh_h2oice)
        h2lift_tot = abs(dh_dust)
        do while ((h2subl_co2ice_tot > 0._dp .or. h2subl_h2oice_tot > 0._dp .or. h2lift_tot > 0._dp) .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._dp
                currentloc => current%up
                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
                    h_toplag = h_toplag + thickness(currentloc)
                    currentloc => currentloc%up
                end do
                if (currentloc%h_h2oice >= h_patchy_ice) then
                    R_subl = 0._dp
                else if (0._dp < 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)
                end if
            end if
           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl

            ! Is there CO2 ice to sublimate?
            h2subl_co2ice = 0._dp
            if (h_co2ice_old > 0._dp .and. h2subl_co2ice_tot > 0._dp) then
                h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
                h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
            end if
            ! Is there H2O ice to sublimate?
            h2subl_h2oice = 0._dp
            if (h_h2oice_old > 0._dp .and. h2subl_h2oice_tot > 0._dp) then
                h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old)
                h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice
            end if
            ! Ice sublimation
            if (h2subl_co2ice > 0._dp .or. h2subl_h2oice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)

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

            ! 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._dp .or. h2lift_tot > 0._dp) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol) then
                current => current%down
                new_lag = .true.
            else
                exit
            end if
        end do
        if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
        if (h2subl_h2oice_tot > 0._dp) dh_h2oice = 0._dp ! No H2O ice available anymore so the tendency is set to 0
        if (h2lift_tot > 0._dp) dh_dust = 0._dp ! No dust available anymore so the tendency is set to 0
!~~~~~~~~~~
! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
else if (dh_co2ice <= 0._dp .and. dh_h2oice >= 0._dp) then
    if (dh_dust >= 0._dp) then
        dh_dust_co2 = 0._dp
        dh_dust_h2o = dh_dust
    else
        dh_dust_co2 = dh_dust
        dh_dust_h2o = 0._dp
    end if
    if (abs(dh_co2ice) > dh_h2oice) then ! CO2 ice sublimation is dominant
        ! CO2 ICE SUBLIMATION: retreating a stratum
        h2subl_co2ice_tot = abs(dh_co2ice)
        h2lift_tot = abs(dh_dust_co2)
        do while ((h2subl_co2ice_tot > 0._dp .or. h2lift_tot > 0._dp) .and. associated(current))
            h_co2ice_old = current%h_co2ice

            ! 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._dp
                currentloc => current%up
                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
                    h_toplag = h_toplag + thickness(currentloc)
                    currentloc => currentloc%up
                end do
                if (currentloc%h_h2oice >= h_patchy_ice) then
                    R_subl = 0._dp
                else if (0._dp < 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)
                end if
            end if
           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl

            ! Is there CO2 ice to sublimate?
            h2subl_co2ice = 0._dp
            if (h_co2ice_old > 0._dp .and. h2subl_co2ice_tot > 0._dp) then
                h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
                h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
            end if
            ! Ice sublimation
            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag)

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

            ! 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._dp .or. h2lift_tot > 0._dp) .and. current%h_co2ice < tol) then
                current => current%down
                new_lag = .true.
            else
                exit
            end if
        end do
        if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
        if (h2lift_tot > 0._dp) dh_dust = 0._dp ! No dust available anymore so the tendency is set to 0

        ! H2O ICE CONDENSATION: building a stratum
        h_pore = dh_h2oice*h2oice_porosity/(1._dp - h2oice_porosity)
        h_tot = dh_h2oice + dh_dust_h2o + h_pore
        ! New stratum at the top of the layering
        if (new_str) then
            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,0._dp)
            new_str = .false.
        else
            this%top%top_elevation = this%top%top_elevation + h_tot
            this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
            this%top%h_dust        = this%top%h_dust        + dh_dust_h2o
            this%top%h_pore        = this%top%h_pore        + h_pore
        end if
    else ! H2O ice condensation is dominant
        ! H2O ICE CONDENSATION: building a stratum
        h_pore = dh_h2oice*h2oice_porosity/(1._dp - h2oice_porosity)
        h_tot = dh_h2oice + dh_dust_h2o + h_pore
        ! New stratum at the top of the layering
        if (new_str) then
            call add_stratum(this,this%top%top_elevation + h_tot,0._dp,dh_h2oice,dh_dust_h2o,h_pore,0._dp)
            new_str = .false.
        else
            this%top%top_elevation = this%top%top_elevation + h_tot
            this%top%h_h2oice      = this%top%h_h2oice      + dh_h2oice
            this%top%h_dust        = this%top%h_dust        + dh_dust_h2o
            this%top%h_pore        = this%top%h_pore        + h_pore
        end if

        ! CO2 ICE SUBLIMATION: retreating a stratum
        h2subl_co2ice_tot = abs(dh_co2ice)
        h2lift_tot = abs(dh_dust_co2)
        do while ((h2subl_co2ice_tot > 0._dp .or. h2lift_tot > 0._dp) .and. associated(current))
            h_co2ice_old = current%h_co2ice

            ! How much is CO2 ice sublimation inhibited by the top dust lag?
            R_subl = 1._dp
            if (associated(current%up)) then ! If there is an upper stratum
                h_toplag = 0._dp
                currentloc => current%up
                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
                    h_toplag = h_toplag + thickness(currentloc)
                    currentloc => currentloc%up
                end do
                if (currentloc%h_h2oice >= h_patchy_ice) then
                    R_subl = 0._dp
                else if (0._dp < 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)
                end if
            end if
           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl

            ! Is there CO2 ice to sublimate?
            h2subl_co2ice = 0._dp
            if (h_co2ice_old > 0._dp .and. h2subl_co2ice_tot > 0._dp) then
                h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old)
                h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice
            end if
            ! Ice sublimation
            if (h2subl_co2ice > 0._dp) call sublimate_ice(this,current,h2subl_co2ice,0._dp,new_lag,zlag)

            ! Is there dust to lift?
            h2lift = 0._dp
            if (is_dust_lag(current) .and. h2lift_tot > 0._dp) then
                h2lift = min(h2lift_tot,current%h_dust)
                h2lift_tot = h2lift_tot - h2lift
                ! Dust lifting
                if (h2lift > 0._dp) 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)
                end if
            end if

            ! 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._dp .or. h2lift_tot > 0._dp) .and. current%h_co2ice < tol) then
                current => current%down
                new_lag = .true.
            else
                exit
            end if
        end do
        if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
        if (h2lift_tot > 0._dp) dh_dust = 0._dp ! No dust available anymore so the tendency is set to 0
    end if
!~~~~~~~~~~
! NOT EXPECTED SITUATION
else
    unexpected = .true.
end if

if (unexpected) then
    call print_err('    dh_co2ice [m/y] = '//real2str(dh_co2ice))
    call print_err('    dh_h2oice [m/y] = '//real2str(dh_h2oice))
    call print_err('    dh_dust   [m/y] = '//real2str(dh_dust))
    call stop_clean(__FILE__,__LINE__,'situation for the layering construction not managed!',1)
end if

zshift_surf = this%top%top_elevation - zshift_surf

END SUBROUTINE evolve_layering
!=======================================================================

END MODULE layered_deposits