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