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,SAVE :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium LOGICAL,SAVE :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) real,save,allocatable :: porefillingice_depth(:,:) ! ngrid x nslope: Depth of the ice table [m] real,save,allocatable :: porefillingice_thickness(:,:) ! ngrid x nslope: Thickness of the ice table [m] contains subroutine ini_ice_table_porefilling(ngrid,nslope) implicit none integer,intent(in) :: ngrid ! number of atmospheric columns integer,intent(in) :: nslope ! number of slope within a mesh allocate(porefillingice_depth(ngrid,nslope)) allocate(porefillingice_thickness(ngrid,nslope)) end subroutine ini_ice_table_porefilling subroutine end_ice_table_porefilling implicit none if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) end subroutine end_ice_table_porefilling SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,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 #ifndef CPP_STD 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_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM logical,intent(in) :: watercaptag(ngrid) ! Boolean to check the presence of a perennial glacier real,intent(in) :: rhowatersurf_ave(ngrid,nslope) ! Water density at the surface, yearly averaged [kg/m^3] real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope) ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3] real,intent(in) :: regolith_inertia(ngrid,nslope) ! Thermal inertia of the regolith layer [SI] ! Ouputs ! ----------- real,intent(out) :: ice_table_beg(ngrid,nslope) ! ice table depth [m] real,intent(out) :: ice_table_thickness(ngrid,nslope) ! ice table thickness [m] ! Local ! ----------- real :: z1,z2 ! intermediate variables used when doing a linear interpolation between two depths to find the root integer ig, islope,isoil,isoilend ! loop variables real :: diff_rho(nsoil_PEM) ! 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 :: previous_icetable_depth(ngrid,nslope) ! 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_PEM) ! 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_PEM 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_PEM-1 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.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_PEM -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).lt.0).and.(diff_rho(isoil+1).gt.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_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) do isoilend = isoil+1,nsoil_PEM-1 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.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).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.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 !======================================================================= RETURN #endif END 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,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice real,intent(in) :: psat_soil(nsoilmx_PEM) ! 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, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase ! local real :: eta(nsoilmx_PEM) ! constriction integer :: ilay ! index for loop real :: old_depth_filling ! depth_filling saved real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn integer :: index_tmp ! for loop 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 :: dz_eta(nsoilmx_PEM) ! \partial z \eta real :: dz_eta_low ! same but evaluated at the interface for ice real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal ! constant real :: pvap2rho = 18.e-3/8.314 ! ! 0. Compute constriction over the layer ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 if (index_IS.lt.0) then index_tmp = nsoilmx_PEM do ilay = 1,nsoilmx_PEM call constriction(porefill(ilay),eta(ilay)) enddo else index_tmp = index_IS do ilay = 1,index_IS-1 call constriction(porefill(ilay),eta(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).le.0) then index_filling = ilay exit endif enddo if(index_filling.eq.1) depth_filling = mlayer_PEM(1) if(index_filling.gt.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).le.0.) then index_firstice = ilay ! first point with ice exit endif enddo ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta) if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) endif 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.gt.0) index_geothermal = -1 if(index_geothermal.lt.0) depth_geothermal = -1. if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010) index_geothermal = -1 do ilay=2,nsoilmx_PEM if (dz_psat(ilay).gt.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.gt.0).and.(index_IS.lt.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 (massfillafterindex_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.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.