source: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 @ 3563

MODULE ice_table_mod

implicit none

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static)
!!! Author:  LL, 02/2023
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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       ! ngrid x nslope: Depth of the ice table [m]
real, allocatable, dimension(:,:)   :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
real, allocatable, dimension(:,:,:) :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]

!-----------------------------------------------------------------------
contains

!-----------------------------------------------------------------------
SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)

implicit none

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

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

implicit none

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_beg,ice_table_thickness)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
!!! density at the surface and at depth.
!!! Computations are made following the methods in Schorgofer et al., 2005
!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
!!!
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use math_mod,                   only: findroot
use comsoil_h_PEM,              only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
use soil_thermalproperties_mod, only: ice_thermal_properties

implicit none

! Inputs
! ------
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]
! Ouputs
! ------
real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
! Locals
! ------
integer                       :: ig, islope, isoil, isoilend ! loop variables
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_beg
do ig = 1,ngrid
    if (watercaptag(ig)) then
        ice_table_beg(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_beg(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_beg(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),mlayer_PEM(isoilend + 1),ice_table_end)
                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(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_beg(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),mlayer_PEM(isoilend + 1),ice_table_end)
                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(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_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
            call ice_thermal_properties(.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_beg(ig,islope) - &
                                             previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
!!!
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use comsoil_h_PEM,         only: mlayer_PEM
use comslope_mod,          only: subslope_dist, def_slope_mean
use constants_marspem_mod, only: porosity
#ifndef CPP_STD
    use comcstfi_h,   only: pi
#else
    use comcstfi_mod, only: pi
#endif

implicit none

! Inputs
! ------
integer,                             intent(in) :: ngrid, nslope, nsoil   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
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]
! Outputs
! -------
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]
! Locals
!-------
integer :: ig, islope, isoil ! loop index
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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
!!!
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM
use math_mod,      only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple

implicit none

! Inputs
! ------
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]
! Outputs
! -------
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
! Locals
!-------
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
! Constant
!---------
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)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the constriction of vapor flux by pore ice
!!!
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

implicit none

real, intent(in) :: porefill ! pore filling fraction
real :: eta ! constriction

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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
!!!
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
use abort_pem_mod, only: abort_pem

implicit none

! Inputs
! ------
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)

! Outputs
! -------
real, intent(out) :: Tice ! Ice temperatures (K)

! Locals
! ------
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 abort_pem("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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute subsurface ice density
!!!
!!! Author: JBC
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use abort_pem_mod, only: abort_pem

implicit none

real,         intent(in) :: T
character(3), intent(in) :: ice
real :: rho

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 abort_pem("rho_ice","Type of ice not known!",1)
end select

END FUNCTION rho_ice

END MODULE ice_table_mod