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