source: trunk/LMDZ.COMMON/libf/evolution/ice_table.F90

MODULE ice_table
!-----------------------------------------------------------------------
! NAME
!     ice_table
!
! DESCRIPTION
!     Ice table variables and methods to compute it (dynamic and static).
!
! AUTHORS & DATE
!     L. Lange, 02/2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! MODULE VARIABLES
! ----------------
logical                             :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
logical                             :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
real, allocatable, dimension(:,:)   :: icetable_depth       ! Depth of the ice table [m]
real, allocatable, dimension(:,:)   :: icetable_thickness   ! Thickness of the ice table [m]
real, allocatable, dimension(:,:,:) :: ice_porefilling      ! Amount of porefilling in each layer in each grid [m^3/m^3]

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
!-----------------------------------------------------------------------
! NAME
!     ini_ice_table
!
! DESCRIPTION
!     Allocate module arrays for ice table fields.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
integer, intent(in) :: ngrid  ! Number of atmospheric columns
integer, intent(in) :: nslope ! Number of slope within a mesh
integer, intent(in) :: nsoil  ! Number of soil layers

! CODE
! ----
allocate(icetable_depth(ngrid,nslope))
allocate(icetable_thickness(ngrid,nslope))
allocate(ice_porefilling(ngrid,nsoil,nslope))

END SUBROUTINE ini_ice_table
!=======================================================================

!=======================================================================
SUBROUTINE end_ice_table
!-----------------------------------------------------------------------
! NAME
!     end_ice_table
!
! DESCRIPTION
!     Deallocate ice table arrays.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! CODE
! ----
if (allocated(icetable_depth)) deallocate(icetable_depth)
if (allocated(icetable_thickness)) deallocate(icetable_thickness)
if (allocated(ice_porefilling)) deallocate(ice_porefilling)

END SUBROUTINE end_ice_table
!=======================================================================

!=======================================================================
SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_depth,ice_table_thickness)
!-----------------------------------------------------------------------
! NAME
!     computeice_table_equilibrium
!
! DESCRIPTION
!     Compute the ice table depth knowing the yearly average water 
!     density at the surface and at depth. Computations are made following
!     the methods in Schorghofer et al., 2005.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 12/12/2025
!
! NOTES
!     This subroutine only gives the ice table at equilibrium and does
!     not consider exchange with the atmosphere.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use math_toolkit,         only: findroot
use soil,                 only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
use soil_thermal_inertia, only: get_ice_TI

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

! ARGUMENTS
! ---------
integer,                                intent(in)  :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
logical, dimension(ngrid),              intent(in)  :: watercaptag          ! Boolean to check the presence of a perennial glacier
real,    dimension(ngrid,nslope),       intent(in)  :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
real,    dimension(ngrid,nsoil,nslope), intent(in)  :: rhowatersoil_ave     ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
real,    dimension(ngrid,nslope),       intent(in)  :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
real,    dimension(ngrid,nslope),       intent(out) :: ice_table_depth      ! Ice table depth [m]
real,    dimension(ngrid,nslope),       intent(out) :: ice_table_thickness  ! Ice table thickness [m]

! LOCAL VARIABLES
! ---------------
integer                       :: ig, islope, isoil, isoilend
real, dimension(nsoil)        :: diff_rho                    ! Difference of water vapor density between the surface and at depth [kg/m^3]
real                          :: ice_table_end               ! Depth of the end of the ice table  [m]
real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
real                          :: stretch                     ! Stretch factor to improve the convergence of the ice table
real                          :: wice_inertia                ! Water Ice thermal Inertia [USI]

! CODE
! ----
previous_icetable_depth = ice_table_depth
do ig = 1,ngrid
    if (watercaptag(ig)) then
        ice_table_depth(ig,:) = 0.
        ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    else
        do islope = 1,nslope
            ice_table_depth(ig,islope) = -1.
            ice_table_thickness(ig,islope) = 0.
            do isoil = 1,nsoil
                diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
            enddo
            if (diff_rho(1) > 0.) then ! ice is at the surface
                ice_table_depth(ig,islope) = 0.
                do isoilend = 2,nsoil - 1
                    if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end)
                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                        exit
                    endif
                enddo
            else
                do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
                    if (diff_rho(isoil) < 0. .and. diff_rho(isoil + 1) > 0.) then
                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_depth(ig,islope))
                        ! Now let's find the end of the ice table
                        ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
                        do isoilend = isoil + 1,nsoil - 1
                            if (diff_rho(isoilend) > 0. .and. diff_rho(isoilend + 1) < 0.) then
                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end)
                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                                exit
                            endif
                        enddo
                    endif ! diff_rho(z) <0 & diff_rho(z+1) > 0
                enddo
            endif ! diff_rho(1) > 0
        enddo
    endif ! watercaptag
enddo

! Small trick to speed up the convergence, Oded's idea.
do islope = 1,nslope
    do ig = 1,ngrid
        if (ice_table_depth(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
            call get_ice_TI(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_depth(ig,islope) - &
                                             previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
            ice_table_depth(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch
        endif
    enddo
enddo

END SUBROUTINE computeice_table_equilibrium
!=======================================================================

!=======================================================================
SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil,icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_m_h2o)
!-----------------------------------------------------------------------
! NAME
!     compute_massh2o_exchange_ssi
!
! DESCRIPTION
!     Compute the mass of H2O that has sublimated from the ice table /
!     condensed.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,                  only: mlayer_PEM
use comslope_mod,          only: subslope_dist, def_slope_mean
use constants_marspem_mod, only: porosity
use phys_constants,        only: pi

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


! ARGUMENTS
! ---------
integer,                                intent(in)  :: ngrid, nslope, nsoil
real,    dimension(ngrid,nslope),       intent(in)  :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
real,    dimension(ngrid,nsoil,nslope), intent(in)  :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
real,    dimension(ngrid,nslope),       intent(in)  :: tsurf                  ! Surface temperature [K]
real,    dimension(ngrid,nsoil,nslope), intent(in)  :: tsoil                  ! Soil temperature [K]
real,    dimension(ngrid),              intent(out) :: delta_m_h2o            ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]

! LOCAL VARIABLES
! ---------------
integer :: ig, islope, isoil
real    :: rho  ! Density of water ice [kg/m^3]
real    :: Tice ! Ice temperature [k]

! CODE
! ----
! Let's compute the amount of ice that has sublimated in each subslope
if (icetable_equilibrium) then
    delta_m_h2o = 0.
    do ig = 1,ngrid
        do islope = 1,nslope
            call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),icetable_depth(ig,islope),Tice)
            delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho_ice(Tice,'h2o')*(icetable_thickness(ig,islope) - icetable_thickness_old(ig,islope)) & ! convention > 0. <=> it condenses
                              *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
        enddo
    enddo
else if (icetable_dynamic) then
    delta_m_h2o = 0.
    do ig = 1,ngrid
        do islope = 1,nslope
            do isoil = 1,nsoil
                call compute_Tice_pem(nsoil,tsoil(ig,:,islope),tsurf(ig,islope),mlayer_PEM(isoil - 1),Tice)
                delta_m_h2o(ig) = delta_m_h2o(ig) + rho_ice(Tice,'h2o')*(ice_porefilling(ig,isoil,islope) - ice_porefilling_old(ig,isoil,islope)) & ! convention > 0. <=> it condenses
                                  *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
            enddo
        enddo
    enddo
endif

END SUBROUTINE compute_massh2o_exchange_ssi
!=======================================================================

!=======================================================================
SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
!-----------------------------------------------------------------------
! NAME
!     find_layering_icetable
!
! DESCRIPTION
!     Compute layering between dry soil, pore filling ice and ice
!     sheet based on Schorghofer, Icarus (2010). Adapted from NS MSIM.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,         only: nsoilmx_PEM, mlayer_PEM
use math_toolkit, only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple

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

! ARGUMENTS
! ---------
real, dimension(nsoilmx_PEM), intent(in)    :: porefill         ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
real, dimension(nsoilmx_PEM), intent(in)    :: psat_soil        ! Soil water pressure at saturation, yearly averaged [Pa]
real,                         intent(in)    :: psat_surf        ! Surface water pressure at saturation, yearly averaged [Pa]
real,                         intent(in)    :: pwat_surf        ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
real,                         intent(in)    :: psat_bottom      ! Boundary conditions for soil vapor pressure [Pa]
real,                         intent(in)    :: B                ! Constant (Eq. 8 from  Schorgofer, Icarus (2010).)
integer,                      intent(in)    :: index_IS         ! Index of the soil layer where the ice sheet begins [1]
real,                         intent(inout) :: depth_filling    ! Depth where pore filling begins [m]
integer,                      intent(out)   :: index_filling    ! Index where the pore filling begins [1]
integer,                      intent(out)   :: index_geothermal ! Index where the ice table stops [1]
real,                         intent(out)   :: depth_geothermal ! Depth where the ice table stops [m]
real, dimension(nsoilmx_PEM), intent(out)   :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase

! LOCAL VARIABLES
! ---------------
real, dimension(nsoilmx_PEM) :: eta                          ! Constriction
real, dimension(nsoilmx_PEM) :: dz_psat                      ! First derivative of the vapor pressure at saturation
real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
integer                      :: ilay, index_tmp              ! Index for loop
real                         :: old_depth_filling            ! Depth_filling saved
real                         :: Jdry                         ! Flux trought the dry layer
real                         :: Jsat                         ! Flux trought the ice layer
real                         :: Jdry_prevlay, Jsat_prevlay   ! Same but for the previous ice layer
integer                      :: index_firstice               ! First index where ice appears (i.e., f > 0)
real                         :: massfillabove, massfillafter ! H2O mass above and after index_geothermal
real, parameter              :: pvap2rho = 18.e-3/8.314

! CODE
! ----
! 0. Compute constriction over the layer
! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
if (index_IS < 0) then
    index_tmp = nsoilmx_PEM
    do ilay = 1,nsoilmx_PEM
        eta(ilay) = constriction(porefill(ilay))
    enddo
else
    index_tmp = index_IS
    do ilay = 1,index_IS - 1
        eta(ilay) = constriction(porefill(ilay))
    enddo
    do ilay = index_IS,nsoilmx_PEM
        eta(ilay) = 0.
    enddo
endif

! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)
old_depth_filling = depth_filling

call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)

do ilay = 1,index_tmp
    Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
    Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
    if (Jdry - Jsat <= 0) then
        index_filling = ilay
        exit
    endif
enddo

if (index_filling == 1) then
    depth_filling = mlayer_PEM(1)
else if (index_filling > 1) then
    Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1)
    Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1)
    call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling)
endif

! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
! 2.0 preliminary: depth to shallowest ice (discontinuity at interface)
index_firstice = -1
do ilay = 1,nsoilmx_PEM
    if (porefill(ilay) <= 0.) then
        index_firstice = ilay  ! first point with ice
        exit
    endif
enddo

! 2.1: now we can compute
call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta)

if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))

call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat)

! 3. Ice table boundary due to geothermal heating
if (index_IS > 0) index_geothermal = -1
if (index_geothermal < 0) depth_geothermal = -1.
if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
    index_geothermal = -1
    do ilay = 2,nsoilmx_PEM
        if (dz_psat(ilay) > 0.) then  ! first point with reversed flux
            index_geothermal = ilay
            call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
            exit
        endif
    enddo
else
    index_geothermal = -1
endif
if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010)
    call  colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,index_geothermal - 1,nsoilmx_PEM,massfillabove)
    index_tmp = -1
    do ilay = index_geothermal,nsoilmx_PEM
        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full
        call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
        if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG
            if (ilay > index_geothermal) then
!                write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
                call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, &
                              dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal)
!                if (depth_geothermal > mlayer_PEM(ilay) .or. depth_geothermal < mlayer_PEM(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay - 1),depth_geothermal,mlayer_PEM(ilay)
              index_tmp = ilay
           endif
           ! otherwise leave depth_geothermal unchanged
           exit
        endif
        massfillabove = massfillafter
    enddo
    if (index_tmp > 0) index_geothermal = index_tmp
end if

END SUBROUTINE find_layering_icetable
!=======================================================================

!=======================================================================
FUNCTION constriction(porefill) RESULT(eta)
!-----------------------------------------------------------------------
! NAME
!     constriction
!
! DESCRIPTION
!     Compute the constriction of vapor flux by pore ice.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

! ARGUMENTS
! ---------
real, intent(in) :: porefill ! pore filling fraction

! LOCAL VARIABLES
! ---------------
real :: eta ! constriction

! CODE
! ----
if (porefill <= 0.) then
    eta = 1.
else if (0 < porefill .and. porefill < 1.) then
    eta = (1-porefill)**2 ! Hudson et al., JGR, 2009
else
    eta = 0.
endif

END FUNCTION constriction
!=======================================================================

!=======================================================================
SUBROUTINE compute_Tice_pem(nsoil,ptsoil,ptsurf,ice_depth,Tice)
!-----------------------------------------------------------------------
! NAME
!     compute_Tice_pem
!
! DESCRIPTION
!     Compute subsurface ice temperature by interpolating the
!     temperatures between the two adjacent cells.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 12/12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use soil,     only: mlayer_PEM
use stop_pem, only: stop_clean

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

! ARGUMENTS
! ---------
integer,                intent(in)  :: nsoil     ! Number of soil layers
real, dimension(nsoil), intent(in)  :: ptsoil    ! Soil temperature (K)
real,                   intent(in)  :: ptsurf    ! Surface temperature (K)
real,                   intent(in)  :: ice_depth ! Ice depth (m)
real,                   intent(out) :: Tice      ! Ice temperatures (K)

! LOCAL VARIABLES
! ---------------
integer :: ik       ! Loop variables
integer :: indexice ! Index of the ice

! CODE
! ----
indexice = -1
if (ice_depth >= mlayer_PEM(nsoil - 1)) then
    Tice = ptsoil(nsoil)
else
    if(ice_depth < mlayer_PEM(0)) then
        indexice = 0.
    else
        do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations
            if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then
                indexice = ik + 1
                exit
            endif
        enddo
    endif
    if (indexice < 0) then
        call stop_clean("compute_Tice_pem","Subsurface ice is below the last soil layer!",1)
    else
        if(indexice >= 1) then ! Linear inteprolation between soil temperature
            Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1)
        else ! Linear inteprolation between the 1st soil temperature and the surface temperature
            Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1)
        endif ! index ice >= 1
    endif !indexice < 0
endif ! icedepth > mlayer_PEM(nsoil - 1)

END SUBROUTINE compute_Tice_pem
!=======================================================================

!=======================================================================
FUNCTION rho_ice(T,ice) RESULT(rho)
!-----------------------------------------------------------------------
! NAME
!     rho_ice
!
! DESCRIPTION
!     Compute subsurface ice density.
!
! AUTHORS & DATE
!     JB Clement, 12/12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
real,         intent(in) :: T
character(3), intent(in) :: ice

! LOCAL VARIABLES
! ---------------
real :: rho

! CODE
! ----
select case (trim(adjustl(ice)))
    case('h2o')
        rho = -3.5353e-4*T**2 + 0.0351*T + 933.5030 ! Rottgers 2012
    case('co2')
        rho = (1.72391 - 2.53e-4*T-2.87*1e-7*T**2)*1.e3 ! Mangan et al. 2017
    case default
        call stop_clean("rho_ice","Type of ice not known!",1)
end select

END FUNCTION rho_ice
!=======================================================================

END MODULE ice_table