MODULE layering_mod

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


implicit none

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

! Physical parameters
real, parameter :: tend_dust = 5.78e-2          ! Tendency of dust [kg.m-2.y-1]
real, parameter :: dry_lag_porosity = 0.5       ! Porosity of dust lag
real, parameter :: dry_regolith_porosity = 0.45 ! Porosity of regolith
real, parameter :: co2ice_porosity = 0.3        ! Porosity of CO2 ice
real, parameter :: h2oice_porosity = 0.3        ! Porosity of H2O ice
real, parameter :: rho_co2ice = 1650.           ! Density of CO2 ice [kg.m-3]
real, parameter :: rho_h2oice = 920.            ! Density of H2O ice [kg.m-3]
real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]

! Stratum = node of the linked list
type :: stratum
    real               :: thickness          ! Layer thickness [m]
    real               :: elevation          ! Layer elevation (top height from the surface) [m]
    real               :: co2ice_volfrac     ! CO2 ice volumic fraction
    real               :: h2oice_volfrac     ! H2O ice volumic fraction
    real               :: dust_volfrac       ! Dust volumic fraction
    real               :: air_volfrac        ! Air volumic fraction inside pores
    type(stratum), pointer :: up => null()   ! Upper stratum (next node)
    type(stratum), pointer :: down => null() ! Lower stratum (previous node)
end type stratum

! Layering = linked list
type layering
    integer                :: nb_str = 0       ! Number of strata
    type(stratum), pointer :: bottom => null() ! First stratum (head of the list)
    type(stratum), pointer :: top => null()    ! Last stratum (tail of the list)
end type layering

!=======================================================================
contains
!=======================================================================
! To display a stratum
SUBROUTINE print_stratum(this)

implicit none

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

!---- Code
write(*,'(a20,es13.6)') 'Thickness      : ', this%thickness
write(*,'(a20,es13.6)') 'Elevation      : ', this%elevation
write(*,'(a20,es13.6)') 'CO2 ice content: ', this%co2ice_volfrac
write(*,'(a20,es13.6)') 'H2O ice content: ', this%h2oice_volfrac
write(*,'(a20,es13.6)') 'Dust content   : ', this%dust_volfrac
write(*,'(a20,es13.6)') 'Air content    : ', this%air_volfrac

END SUBROUTINE print_stratum

!=======================================================================
! To display a layering
SUBROUTINE print_lay(this)

implicit none

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

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

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

i = this%nb_str
do while (associated(current))
    write(*,'(a8,i4)') 'Stratum ', i
    call print_stratum(current)
    write(*,*)
    current => current%down
    i = i - 1
enddo

END SUBROUTINE print_lay

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

implicit none

!---- Arguments
type(layering), intent(inout)  :: this
real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air

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

!---- Code
! Creation of the stratum
allocate(str)
nullify(str%up,str%down)
str%thickness = 1.
str%elevation = 1.
str%co2ice_volfrac = 0.
str%h2oice_volfrac = 0.
str%dust_volfrac = 0.
str%air_volfrac = 1.
if (present(thickness)) str%thickness = thickness
if (present(elevation)) str%elevation = elevation
if (present(co2ice)) str%co2ice_volfrac = co2ice
if (present(h2oice)) str%h2oice_volfrac = h2oice
if (present(dust)) str%dust_volfrac = dust
if (present(air)) str%air_volfrac = air

! Verification of volume fraction
if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
     error stop 'add_stratum: properties for the new stratum are not possible (sum of volumic fraction /= 1)!'

! Increment the number of strata
this%nb_str = this%nb_str + 1

! Add the stratum to the layering
if (.not. associated(this%bottom)) then ! If it is the very first one
    this%bottom => str
    this%top => str
else
    str%down => this%top
    this%top%up => str
    this%top => str
endif

END SUBROUTINE add_stratum

!=======================================================================
! To insert a stratum after another one in a layering
SUBROUTINE insert_stratum(this,str_prev,thickness,elevation,co2ice,h2oice,dust,air)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str_prev
real, intent(in), optional :: thickness, elevation, co2ice, h2oice, dust, air

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

!---- Code
if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
    call add_stratum(this,thickness,elevation,co2ice,h2oice,dust,air)
else
    ! Creation of the stratum
    allocate(str)
    nullify(str%up,str%down)
    str%thickness = 1.
    str%elevation = 1.
    str%co2ice_volfrac = 0.
    str%h2oice_volfrac = 0.
    str%dust_volfrac = 0.
    str%air_volfrac = 1.
    if (present(thickness)) str%thickness = thickness
    if (present(elevation)) str%elevation = elevation
    if (present(co2ice)) str%co2ice_volfrac = co2ice
    if (present(h2oice)) str%h2oice_volfrac = h2oice
    if (present(dust)) str%dust_volfrac = dust
    if (present(air)) str%air_volfrac = air

    ! Verification of volume fraction
    if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%air_volfrac)) > tol) &
         error stop 'insert_stratum: properties for the new stratum are not possible (sum of volumic fraction /= 1)!'

    ! Increment the number of strata
    this%nb_str = this%nb_str + 1

    ! Insert the stratum into the layering at the right place
    str%down => str_prev
    str%up => str_prev%up
    str_prev%up => str
    str%up%down => str

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

END SUBROUTINE insert_stratum

!=======================================================================
! To remove a stratum in a layering
SUBROUTINE rm_stratum(this,str)

implicit none

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

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

!---- Code
! Protect the 3 sub-surface strata from removing
if (str%elevation <= 0.) return

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

! Remove the stratum from the stratification
str%down%up => str%up
if (associated(str%up)) then ! If it is not the last one at the top of the layering
    str%up%down => str%down
else
    this%top => str%down
endif
new_str => str%down
nullify(str%up,str%down)
deallocate(str)
nullify(str)
str => new_str

END SUBROUTINE rm_stratum

!=======================================================================
! To initialize the layering
SUBROUTINE ini_layering(this)

implicit none

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

!---- Code
! Creation of three strata at the bottom of the layering
! 1) Dry porous deep regolith
call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)
! 2) Ice-cemented regolith (ice table)
call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.)
! 3) Dry porous regolith
call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity)

END SUBROUTINE ini_layering

!=======================================================================
! To delete the entire layering
SUBROUTINE del_layering(this)

implicit none

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

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

!---- Code
current => this%bottom
do while (associated(current))
    next => current%up
    deallocate(current)
    nullify(current)
    current => next
enddo
this%nb_str = 0

END SUBROUTINE del_layering

!=======================================================================
! To get a specific stratum in a layering
FUNCTION get_stratum(lay,i) RESULT(str)

implicit none

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

!---- Local variables
integer :: k

!---- Code
if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
k = 1
str => lay%bottom
do while (k /= i)
    str => str%up
    k = k + 1
enddo

END FUNCTION get_stratum

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

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str
logical,                intent(inout) :: new_lag
real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old

!---- Local variables
real               :: r_subl, hsubl_dust, h_lag
type(stratum), pointer :: current

!---- Code
! How much dust does CO2 ice sublimation involve to create the lag?
r_subl = h2subl/(1. - co2ice_porosity) &
        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
hsubl_dust = r_subl*str%dust_volfrac*str%thickness
! Update of properties in the sublimating stratum
str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
str%elevation = str%elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust
str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness
str%h2oice_volfrac = h_h2oice_old/str%thickness
str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
! Correct the value of elevation for the upper strata
current => str%up
do while (associated(current))
    current%elevation = current%down%elevation + current%thickness
    current => current%up
enddo
! Remove the sublimating stratum if there is no ice anymore
if (str%thickness < tol) call rm_stratum(this,str)
! New stratum = dust lag
h_lag = hsubl_dust/(1. - dry_lag_porosity)
if (h_lag > 0.) then ! Only if the dust lag is building up
    if (new_lag) then
        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
        new_lag = .false.
    else
        str%up%thickness = str%up%thickness + h_lag
        str%up%elevation = str%up%elevation + h_lag
    endif
endif

END SUBROUTINE subl_co2ice_layering

!=======================================================================
! To sublimate H2O ice in the layering
SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: str
logical,                intent(inout) :: new_lag
real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old

!---- Local variables
real               :: r_subl, hsubl_dust, h_lag
type(stratum), pointer :: current

!---- Code
! How much dust does CO2 ice sublimation involve to create the lag?
r_subl = h2subl/(1. - h2oice_porosity) &
        /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity)))
hsubl_dust = r_subl*str%dust_volfrac*str%thickness
! Update of properties in the sublimating stratum
str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
str%elevation = str%elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust
str%co2ice_volfrac = h_co2ice_old/str%thickness
str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness
str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness
str%air_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac)
! Correct the value of elevation for the upper strata
current => str%up
do while (associated(current))
    current%elevation = current%down%elevation + current%thickness
    current => current%up
enddo
! Remove the sublimating stratum if there is no ice anymore
if (str%thickness < tol) call rm_stratum(this,str)
! New stratum = dust lag
h_lag = hsubl_dust/(1. - dry_lag_porosity)
if (h_lag > 0.) then ! Only if the dust lag is building up
    if (new_lag) then
        call insert_stratum(this,str,h_lag,str%elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity)
        new_lag = .false.
    else
        str%up%thickness = str%up%thickness + h_lag
        str%up%elevation = str%up%elevation + h_lag
    endif
endif

END SUBROUTINE subl_h2oice_layering

!=======================================================================
! To lift dust in the layering
SUBROUTINE lift_dust_lags(this,str,h2lift)

implicit none

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

!---- Code
! Update of properties in the eroding dust lag
str%thickness = str%thickness - h2lift/(1. - str%air_volfrac)
str%elevation = str%elevation - h2lift/(1. - str%air_volfrac)
! Remove the eroding dust lag if there is no dust anymore
if (str%thickness < tol) call rm_stratum(this,str)

END SUBROUTINE lift_dust_lags

!=======================================================================
! Layering algorithm
SUBROUTINE make_layering(this,k,htend_co2ice,htend_h2oice,htend_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2)

implicit none

!---- Arguments
type(layering),         intent(inout) :: this
type(stratum), pointer, intent(inout) :: current1, current2
logical,                intent(inout) :: new_str, new_lag1, new_lag2, stopPEM
real, intent(in) :: htend_co2ice, htend_h2oice, htend_dust
integer,  intent(in) :: k

!---- Local variables
real :: h_co2ice_old, h_h2oice_old, h_dust_old, thickness, h2subl, h2subl_tot, h2lift, h2lift_tot

!---- Code
if (htend_dust < 0.) then ! Dust lifting
    if (abs(htend_co2ice) > eps .or. abs(htend_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!'
    write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting'
    h2lift_tot = abs(htend_dust)
    do while (h2lift_tot > 0. .and. associated(current1))
        ! Is the considered stratum a dust lag?
        if (abs(1. - (current1%dust_volfrac + current1%air_volfrac)) < tol) then ! Yes, we can lift dust
            h_dust_old = current1%dust_volfrac*current1%thickness
            ! How much dust can be lifted in the considered dust lag?
            if (h2lift_tot - h_dust_old <= eps) then ! There is enough to lift
                h2lift = h2lift_tot
                h2lift_tot = 0.
                call lift_dust_lags(this,current1,h2lift)
            else ! Only a fraction can be lifted and so we move to the underlying stratum
                h2lift = h_dust_old
                h2lift_tot = h2lift_tot - h2lift
                call lift_dust_lags(this,current1,h2lift)
                current1 => current1%down
            endif
        else ! No, we stop
            write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust lifting: no dust lag!'
            stopPEM = .true.
            exit
        endif
    enddo
    if (h2lift_tot > 0.) then
        write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough dust in the lag layers to complete the lifting!'
        stopPEM = .true.
    endif

else

!------ Dust sedimentation
    if (abs(htend_co2ice) < eps .and. abs(htend_h2oice) < eps) then
        write(*,'(a5,i3,a)') 'Year ', k, ' -> Dust sedimentation'
        ! New stratum at the layering top by sedimentation of dust
        thickness = htend_dust/(1. - dry_regolith_porosity)
        if (thickness > 0.) then ! Only if the stratum is builiding up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,0.,htend_dust/thickness,dry_regolith_porosity)
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%elevation = this%top%elevation + thickness
             endif
        endif

!------ Condensation of CO2 ice + H2O ice
    else if ((htend_co2ice >= 0. .and. htend_h2oice > 0.) .or. (htend_co2ice > 0. .and. htend_h2oice >= 0.)) then
        write(*,'(a5,i3,a)') 'Year ', k, ' -> CO2 and H2O ice condensation'
        ! New stratum at the layering top by condensation of CO2 and H2O ice
        thickness = htend_co2ice/(1. - co2ice_porosity) + htend_h2oice/(1. - h2oice_porosity) + htend_dust
        if (thickness > 0.) then ! Only if the stratum is builiding up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%elevation + thickness,htend_co2ice/thickness,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_co2ice/thickness + htend_h2oice/thickness + htend_dust/thickness))
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%elevation = this%top%elevation + thickness
             endif
        endif

!------ Sublimation of CO2 ice + Condensation of H2O ice
    else if (htend_co2ice < 0. .and. htend_h2oice >= 0. ) then
        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 ice + Condensation of H2O ice'
        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(htend_co2ice)
        do while (h2subl_tot > 0. .and. associated(current1))
            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
            h_dust_old = current1%dust_volfrac*current1%thickness
            ! How much CO2 ice can sublimate in the considered stratum?
            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
                h2subl = h2subl_tot
                h2subl_tot = 0.
                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
                current1 => current1%down
                new_lag1 = .true.
            else ! Only a fraction can sublimate and so we move to the underlying stratum
                h2subl = h_co2ice_old
                h2subl_tot = h2subl_tot - h2subl
                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
                current1 => current1%down
                new_lag1 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
            stopPEM = .true.
        endif
        ! New stratum at the layering top by condensation of H2O ice
        thickness = htend_h2oice/(1. - h2oice_porosity) + htend_dust
        if (thickness > 0.) then ! Only if the stratum is builiding up
             if (new_str) then
                 call add_stratum(this,thickness,this%top%elevation + thickness,0.,htend_h2oice/thickness,htend_dust/thickness,1. - (htend_h2oice/thickness + htend_dust/thickness))
                 new_str = .false.
             else
                 this%top%thickness = this%top%thickness + thickness
                 this%top%elevation = this%top%elevation + thickness
             endif
        endif

!------ Sublimation of CO2 ice + H2O ice
    else if (htend_co2ice < 0. .and. htend_h2oice < 0.) then
        write(*,'(a5,i3,a)') 'Year ', k, ' -> Sublimation of CO2 and H2O ice'
        ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(htend_co2ice)
        do while (h2subl_tot > 0. .and. associated(current1))
            h_co2ice_old = current1%co2ice_volfrac*current1%thickness
            h_h2oice_old = current1%h2oice_volfrac*current1%thickness
            h_dust_old = current1%dust_volfrac*current1%thickness
            ! How much CO2 ice can sublimate in the considered stratum?
            if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate
                h2subl = h2subl_tot
                h2subl_tot = 0.
                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
            else if (h_co2ice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
                current1 => current1%down
                new_lag1 = .true.
            else ! Only a fraction can sublimate and so we move to the underlying stratum
                h2subl = h_co2ice_old
                h2subl_tot = h2subl_tot - h2subl
                call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag1)
                current1 => current1%down
                new_lag1 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough CO2 ice in the layering to complete the sublimation!'
            stopPEM = .true.
        endif
        ! H2O ice sublimation in the considered stratum + New stratum for dust lag
        h2subl_tot = abs(htend_h2oice)
        do while (h2subl_tot > 0. .and. associated(current2))
            h_co2ice_old = current2%co2ice_volfrac*current2%thickness
            h_h2oice_old = current2%h2oice_volfrac*current2%thickness
            h_dust_old = current2%dust_volfrac*current2%thickness
            ! How much H2O ice can sublimate in the considered stratum?
            if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate
                h2subl = h2subl_tot
                h2subl_tot = 0.
                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
            else if (h_h2oice_old < eps) then ! There is nothing to sublimate so we move to the underlying stratum
                current2 => current2%down
                new_lag2 = .true.
            else ! Only a fraction can sublimate and so we move to the underlying stratum
                h2subl = h_h2oice_old
                h2subl_tot = h2subl_tot - h2subl
                call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,new_lag2)
                current2 => current2%down
                new_lag2 = .true.
            endif
        enddo
        if (h2subl_tot > 0.) then
            write(*,'(a5,i3,a)') 'Year ', k, ' -> Not enough H2O ice in the layering to complete the sublimation!'
            stopPEM = .true.
        endif

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

END SUBROUTINE make_layering

END MODULE layering_mod