| 1 | module ice_table_mod |
|---|
| 2 | |
|---|
| 3 | implicit none |
|---|
| 4 | |
|---|
| 5 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 6 | !!! |
|---|
| 7 | !!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static) |
|---|
| 8 | !!! Author: LL, 02/2023 |
|---|
| 9 | !!! |
|---|
| 10 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 11 | |
|---|
| 12 | LOGICAL,SAVE :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium |
|---|
| 13 | LOGICAL,SAVE :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) |
|---|
| 14 | real,save,allocatable :: porefillingice_depth(:,:) ! ngrid x nslope: Depth of the ice table [m] |
|---|
| 15 | real,save,allocatable :: porefillingice_thickness(:,:) ! ngrid x nslope: Thickness of the ice table [m] |
|---|
| 16 | |
|---|
| 17 | contains |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | subroutine ini_ice_table_porefilling(ngrid,nslope) |
|---|
| 21 | |
|---|
| 22 | implicit none |
|---|
| 23 | integer,intent(in) :: ngrid ! number of atmospheric columns |
|---|
| 24 | integer,intent(in) :: nslope ! number of slope within a mesh |
|---|
| 25 | |
|---|
| 26 | allocate(porefillingice_depth(ngrid,nslope)) |
|---|
| 27 | allocate(porefillingice_thickness(ngrid,nslope)) |
|---|
| 28 | |
|---|
| 29 | end subroutine ini_ice_table_porefilling |
|---|
| 30 | |
|---|
| 31 | subroutine end_ice_table_porefilling |
|---|
| 32 | |
|---|
| 33 | implicit none |
|---|
| 34 | if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) |
|---|
| 35 | if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) |
|---|
| 36 | |
|---|
| 37 | end subroutine end_ice_table_porefilling |
|---|
| 38 | |
|---|
| 39 | !!! -------------------------------------- |
|---|
| 40 | |
|---|
| 41 | SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) |
|---|
| 42 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 43 | !!! |
|---|
| 44 | !!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water |
|---|
| 45 | !!! density at the surface and at depth. |
|---|
| 46 | !!! Computations are made following the methods in Schorgofer et al., 2005 |
|---|
| 47 | !!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere |
|---|
| 48 | !!! |
|---|
| 49 | !!! Author: LL |
|---|
| 50 | !!! |
|---|
| 51 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 52 | use math_mod,only: findroot |
|---|
| 53 | USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM ! Depth of the vertical grid |
|---|
| 54 | USE soil_thermalproperties_mod, only: ice_thermal_properties |
|---|
| 55 | implicit none |
|---|
| 56 | ! inputs |
|---|
| 57 | ! ----------- |
|---|
| 58 | integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM |
|---|
| 59 | logical,intent(in) :: watercaptag(ngrid) ! Boolean to check the presence of a perennial glacier |
|---|
| 60 | real,intent(in) :: rhowatersurf_ave(ngrid,nslope) ! Water density at the surface, yearly averaged [kg/m^3] |
|---|
| 61 | 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] |
|---|
| 62 | real,intent(in) :: regolith_inertia(ngrid,nslope) ! Thermal inertia of the regolith layer [SI] |
|---|
| 63 | ! Ouputs |
|---|
| 64 | ! ----------- |
|---|
| 65 | real,intent(out) :: ice_table_beg(ngrid,nslope) ! ice table depth [m] |
|---|
| 66 | real,intent(out) :: ice_table_thickness(ngrid,nslope) ! ice table thickness [m] |
|---|
| 67 | ! Local |
|---|
| 68 | ! ----------- |
|---|
| 69 | real :: z1,z2 ! intermediate variables used when doing a linear interpolation between two depths to find the root |
|---|
| 70 | integer ig, islope,isoil,isoilend ! loop variables |
|---|
| 71 | real :: diff_rho(nsoil_PEM) ! difference of water vapor density between the surface and at depth [kg/m^3] |
|---|
| 72 | real :: ice_table_end ! depth of the end of the ice table [m] |
|---|
| 73 | real :: previous_icetable_depth(ngrid,nslope) ! Ice table computed at previous ice depth [m] |
|---|
| 74 | real :: stretch ! stretch factor to improve the convergence of the ice table |
|---|
| 75 | real :: wice_inertia ! Water Ice thermal Inertia [USI] |
|---|
| 76 | ! Code |
|---|
| 77 | ! ----------- |
|---|
| 78 | |
|---|
| 79 | previous_icetable_depth(:,:) = ice_table_beg(:,:) |
|---|
| 80 | do ig = 1,ngrid |
|---|
| 81 | if(watercaptag(ig)) then |
|---|
| 82 | ice_table_beg(ig,:) = 0. |
|---|
| 83 | ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) |
|---|
| 84 | else |
|---|
| 85 | do islope = 1,nslope |
|---|
| 86 | ice_table_beg(ig,islope) = -1. |
|---|
| 87 | ice_table_thickness(ig,islope) = 0. |
|---|
| 88 | do isoil = 1,nsoil_PEM |
|---|
| 89 | diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) |
|---|
| 90 | enddo |
|---|
| 91 | if(diff_rho(1) > 0) then ! ice is at the surface |
|---|
| 92 | ice_table_beg(ig,islope) = 0. |
|---|
| 93 | do isoilend = 2,nsoil_PEM-1 |
|---|
| 94 | if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then |
|---|
| 95 | call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end) |
|---|
| 96 | ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) |
|---|
| 97 | exit |
|---|
| 98 | endif |
|---|
| 99 | enddo |
|---|
| 100 | else |
|---|
| 101 | 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 |
|---|
| 102 | if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then |
|---|
| 103 | call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope)) |
|---|
| 104 | ! Now let's find the end of the ice table |
|---|
| 105 | ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) |
|---|
| 106 | do isoilend = isoil+1,nsoil_PEM-1 |
|---|
| 107 | if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then |
|---|
| 108 | call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end) |
|---|
| 109 | ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) |
|---|
| 110 | exit |
|---|
| 111 | endif |
|---|
| 112 | enddo |
|---|
| 113 | endif !diff_rho(z) <0 & diff_rho(z+1) > 0 |
|---|
| 114 | enddo |
|---|
| 115 | endif ! diff_rho(1) > 0 |
|---|
| 116 | enddo |
|---|
| 117 | endif ! watercaptag |
|---|
| 118 | enddo |
|---|
| 119 | |
|---|
| 120 | ! Small trick to speed up the convergence, Oded's idea. |
|---|
| 121 | do islope = 1,nslope |
|---|
| 122 | do ig = 1,ngrid |
|---|
| 123 | if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then |
|---|
| 124 | call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia) |
|---|
| 125 | stretch = (regolith_inertia(ig,islope)/wice_inertia)**2 |
|---|
| 126 | |
|---|
| 127 | ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - & |
|---|
| 128 | previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) |
|---|
| 129 | ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch |
|---|
| 130 | endif |
|---|
| 131 | enddo |
|---|
| 132 | enddo |
|---|
| 133 | |
|---|
| 134 | RETURN |
|---|
| 135 | END |
|---|
| 136 | |
|---|
| 137 | !!! -------------------------------------- |
|---|
| 138 | |
|---|
| 139 | 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) |
|---|
| 140 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 141 | !!! |
|---|
| 142 | !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM |
|---|
| 143 | !!! |
|---|
| 144 | !!! Author: LL |
|---|
| 145 | !!! |
|---|
| 146 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 147 | use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM |
|---|
| 148 | use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple |
|---|
| 149 | implicit none |
|---|
| 150 | ! inputs |
|---|
| 151 | ! ------ |
|---|
| 152 | |
|---|
| 153 | real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice |
|---|
| 154 | real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa] |
|---|
| 155 | real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] |
|---|
| 156 | real,intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] |
|---|
| 157 | real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] |
|---|
| 158 | real,intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) |
|---|
| 159 | integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] |
|---|
| 160 | real, intent(inout) :: depth_filling ! depth where pore filling begins [m] |
|---|
| 161 | |
|---|
| 162 | ! outputs |
|---|
| 163 | ! ------- |
|---|
| 164 | integer,intent(out) :: index_filling ! index where the pore filling begins [1] |
|---|
| 165 | integer, intent(out) :: index_geothermal ! index where the ice table stops [1] |
|---|
| 166 | real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] |
|---|
| 167 | real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase |
|---|
| 168 | |
|---|
| 169 | ! local |
|---|
| 170 | |
|---|
| 171 | real :: eta(nsoilmx_PEM) ! constriction |
|---|
| 172 | integer :: ilay ! index for loop |
|---|
| 173 | real :: old_depth_filling ! depth_filling saved |
|---|
| 174 | real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn |
|---|
| 175 | integer :: index_tmp ! for loop |
|---|
| 176 | real :: Jdry ! flux trought the dry layer |
|---|
| 177 | real :: Jsat ! flux trought the ice layer |
|---|
| 178 | real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer |
|---|
| 179 | integer :: index_firstice ! first index where ice appears (i.e., f > 0) |
|---|
| 180 | real :: dz_eta(nsoilmx_PEM) ! \partial z \eta |
|---|
| 181 | real :: dz_eta_low ! same but evaluated at the interface for ice |
|---|
| 182 | real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat |
|---|
| 183 | real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal |
|---|
| 184 | |
|---|
| 185 | ! constant |
|---|
| 186 | real :: pvap2rho = 18.e-3/8.314 |
|---|
| 187 | ! |
|---|
| 188 | |
|---|
| 189 | ! 0. Compute constriction over the layer |
|---|
| 190 | ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 |
|---|
| 191 | if (index_IS.lt.0) then |
|---|
| 192 | index_tmp = nsoilmx_PEM |
|---|
| 193 | do ilay = 1,nsoilmx_PEM |
|---|
| 194 | call constriction(porefill(ilay),eta(ilay)) |
|---|
| 195 | enddo |
|---|
| 196 | else |
|---|
| 197 | index_tmp = index_IS |
|---|
| 198 | do ilay = 1,index_IS-1 |
|---|
| 199 | call constriction(porefill(ilay),eta(ilay)) |
|---|
| 200 | enddo |
|---|
| 201 | do ilay = index_IS,nsoilmx_PEM |
|---|
| 202 | eta(ilay) = 0. |
|---|
| 203 | enddo |
|---|
| 204 | endif |
|---|
| 205 | |
|---|
| 206 | ! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| 207 | |
|---|
| 208 | old_depth_filling = depth_filling |
|---|
| 209 | |
|---|
| 210 | call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) |
|---|
| 211 | |
|---|
| 212 | do ilay = 1,index_tmp |
|---|
| 213 | Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| 214 | Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| 215 | if((Jdry - Jsat).le.0) then |
|---|
| 216 | index_filling = ilay |
|---|
| 217 | exit |
|---|
| 218 | endif |
|---|
| 219 | enddo |
|---|
| 220 | |
|---|
| 221 | if(index_filling.eq.1) depth_filling = mlayer_PEM(1) |
|---|
| 222 | |
|---|
| 223 | if(index_filling.gt.1) then |
|---|
| 224 | Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1) |
|---|
| 225 | Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1) |
|---|
| 226 | call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling) |
|---|
| 227 | endif |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010)) |
|---|
| 231 | |
|---|
| 232 | |
|---|
| 233 | ! 2.0 preliminary: depth to shallowest ice (discontinuity at interface) |
|---|
| 234 | |
|---|
| 235 | index_firstice = -1 |
|---|
| 236 | do ilay = 1,nsoilmx_PEM |
|---|
| 237 | if (porefill(ilay).le.0.) then |
|---|
| 238 | index_firstice = ilay ! first point with ice |
|---|
| 239 | exit |
|---|
| 240 | endif |
|---|
| 241 | enddo |
|---|
| 242 | |
|---|
| 243 | ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) |
|---|
| 244 | |
|---|
| 245 | call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta) |
|---|
| 246 | |
|---|
| 247 | if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then |
|---|
| 248 | call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) |
|---|
| 249 | endif |
|---|
| 250 | |
|---|
| 251 | call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) |
|---|
| 252 | dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:)) |
|---|
| 253 | |
|---|
| 254 | ! 3. Ice table boundary due to geothermal heating |
|---|
| 255 | |
|---|
| 256 | if(index_IS.gt.0) index_geothermal = -1 |
|---|
| 257 | if(index_geothermal.lt.0) depth_geothermal = -1. |
|---|
| 258 | if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010) |
|---|
| 259 | index_geothermal = -1 |
|---|
| 260 | do ilay=2,nsoilmx_PEM |
|---|
| 261 | if (dz_psat(ilay).gt.0.) then ! first point with reversed flux |
|---|
| 262 | index_geothermal=ilay |
|---|
| 263 | call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) |
|---|
| 264 | exit |
|---|
| 265 | endif |
|---|
| 266 | enddo |
|---|
| 267 | else |
|---|
| 268 | index_geothermal = -1 |
|---|
| 269 | endif |
|---|
| 270 | if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 24 from Schorgofer, Icarus (2010) |
|---|
| 271 | call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) |
|---|
| 272 | index_tmp = -1 |
|---|
| 273 | do ilay=index_geothermal,nsoilmx_PEM |
|---|
| 274 | if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full |
|---|
| 275 | call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) |
|---|
| 276 | if (massfillafter<dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG |
|---|
| 277 | if (ilay>index_geothermal) then |
|---|
| 278 | ! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal |
|---|
| 279 | call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & |
|---|
| 280 | dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) |
|---|
| 281 | ! if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*) |
|---|
| 282 | ! '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay) |
|---|
| 283 | index_tmp=ilay |
|---|
| 284 | endif |
|---|
| 285 | ! otherwise leave depth_geothermal unchanged |
|---|
| 286 | exit |
|---|
| 287 | endif |
|---|
| 288 | massfillabove = massfillafter |
|---|
| 289 | enddo |
|---|
| 290 | if (index_tmp.gt.0) index_geothermal = index_tmp |
|---|
| 291 | end if |
|---|
| 292 | return |
|---|
| 293 | end subroutine |
|---|
| 294 | |
|---|
| 295 | !!! -------------------------------------- |
|---|
| 296 | |
|---|
| 297 | SUBROUTINE constriction(porefill,eta) |
|---|
| 298 | |
|---|
| 299 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 300 | !!! |
|---|
| 301 | !!! Purpose: Compute the constriction of vapor flux by pore ice |
|---|
| 302 | !!! |
|---|
| 303 | !!! Author: LL |
|---|
| 304 | !!! |
|---|
| 305 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 306 | implicit none |
|---|
| 307 | real,intent(in) :: porefill ! pore filling fraction |
|---|
| 308 | real,intent(out) :: eta ! constriction |
|---|
| 309 | |
|---|
| 310 | !!! |
|---|
| 311 | |
|---|
| 312 | if (porefill.le.0.) eta = 1. |
|---|
| 313 | if ((porefill.gt.0.) .and.(porefill.lt.1.)) then |
|---|
| 314 | eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 |
|---|
| 315 | endif |
|---|
| 316 | if (porefill.le.1.) eta = 0. |
|---|
| 317 | return |
|---|
| 318 | end subroutine |
|---|
| 319 | |
|---|
| 320 | end module |
|---|
