| [3320] | 1 | MODULE ice_table_mod |
|---|
| [2888] | 2 | |
|---|
| [3320] | 3 | implicit none |
|---|
| [2888] | 4 | |
|---|
| 5 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 6 | !!! |
|---|
| [2961] | 7 | !!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static) |
|---|
| [2888] | 8 | !!! Author: LL, 02/2023 |
|---|
| 9 | !!! |
|---|
| 10 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 11 | |
|---|
| [3320] | 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, allocatable, dimension(:,:), save :: porefillingice_depth ! ngrid x nslope: Depth of the ice table [m] |
|---|
| 15 | real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m] |
|---|
| [2888] | 16 | |
|---|
| [3321] | 17 | !----------------------------------------------------------------------- |
|---|
| [3327] | 18 | contains |
|---|
| [3321] | 19 | !----------------------------------------------------------------------- |
|---|
| [3320] | 20 | SUBROUTINE ini_ice_table_porefilling(ngrid,nslope) |
|---|
| [2888] | 21 | |
|---|
| [3320] | 22 | implicit none |
|---|
| [2961] | 23 | |
|---|
| [3320] | 24 | integer, intent(in) :: ngrid ! number of atmospheric columns |
|---|
| 25 | integer, intent(in) :: nslope ! number of slope within a mesh |
|---|
| [2961] | 26 | |
|---|
| [3320] | 27 | allocate(porefillingice_depth(ngrid,nslope)) |
|---|
| 28 | allocate(porefillingice_thickness(ngrid,nslope)) |
|---|
| [2961] | 29 | |
|---|
| [3320] | 30 | END SUBROUTINE ini_ice_table_porefilling |
|---|
| [2961] | 31 | |
|---|
| [3321] | 32 | !----------------------------------------------------------------------- |
|---|
| [3320] | 33 | SUBROUTINE end_ice_table_porefilling |
|---|
| [2961] | 34 | |
|---|
| [3320] | 35 | implicit none |
|---|
| [2961] | 36 | |
|---|
| [3320] | 37 | if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) |
|---|
| 38 | if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) |
|---|
| [2961] | 39 | |
|---|
| [3321] | 40 | END SUBROUTINE end_ice_table_porefilling |
|---|
| [2961] | 41 | |
|---|
| [3321] | 42 | !------------------------------------------------------------------------------------------------------ |
|---|
| [3320] | 43 | SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) |
|---|
| [2888] | 44 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 45 | !!! |
|---|
| [2961] | 46 | !!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water |
|---|
| [2888] | 47 | !!! density at the surface and at depth. |
|---|
| 48 | !!! Computations are made following the methods in Schorgofer et al., 2005 |
|---|
| [3320] | 49 | !!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere |
|---|
| 50 | !!! |
|---|
| [2888] | 51 | !!! Author: LL |
|---|
| 52 | !!! |
|---|
| 53 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| [3320] | 54 | use math_mod, only: findroot |
|---|
| 55 | use comsoil_h_PEM, only: mlayer_PEM, layer_PEM ! Depth of the vertical grid |
|---|
| 56 | use soil_thermalproperties_mod, only: ice_thermal_properties |
|---|
| 57 | |
|---|
| 58 | implicit none |
|---|
| 59 | |
|---|
| [3327] | 60 | ! Inputs |
|---|
| 61 | ! ------ |
|---|
| [3320] | 62 | integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM |
|---|
| 63 | logical, dimension(ngrid), intent(in) :: watercaptag ! Boolean to check the presence of a perennial glacier |
|---|
| 64 | real, dimension(ngrid,nslope), intent(in) :: rhowatersurf_ave ! Water density at the surface, yearly averaged [kg/m^3] |
|---|
| 65 | real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3] |
|---|
| 66 | real, dimension(ngrid,nslope), intent(in) :: regolith_inertia ! Thermal inertia of the regolith layer [SI] |
|---|
| [3327] | 67 | ! Ouputs |
|---|
| 68 | ! ------ |
|---|
| [3320] | 69 | real, dimension(ngrid,nslope), intent(out) :: ice_table_beg ! ice table depth [m] |
|---|
| 70 | real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m] |
|---|
| [3327] | 71 | ! Locals |
|---|
| 72 | ! ------ |
|---|
| [3321] | 73 | integer :: ig, islope, isoil, isoilend ! loop variables |
|---|
| 74 | real, dimension(nsoil_PEM) :: diff_rho ! difference of water vapor density between the surface and at depth [kg/m^3] |
|---|
| 75 | real :: ice_table_end ! depth of the end of the ice table [m] |
|---|
| 76 | real, dimension(ngrid,nslope) :: previous_icetable_depth ! Ice table computed at previous ice depth [m] |
|---|
| 77 | real :: stretch ! stretch factor to improve the convergence of the ice table |
|---|
| 78 | real :: wice_inertia ! Water Ice thermal Inertia [USI] |
|---|
| [3327] | 79 | ! Code |
|---|
| 80 | ! ---- |
|---|
| [3321] | 81 | previous_icetable_depth = ice_table_beg |
|---|
| [3320] | 82 | do ig = 1,ngrid |
|---|
| 83 | if (watercaptag(ig)) then |
|---|
| 84 | ice_table_beg(ig,:) = 0. |
|---|
| 85 | ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) |
|---|
| 86 | else |
|---|
| [2888] | 87 | do islope = 1,nslope |
|---|
| [3320] | 88 | ice_table_beg(ig,islope) = -1. |
|---|
| 89 | ice_table_thickness(ig,islope) = 0. |
|---|
| 90 | do isoil = 1,nsoil_PEM |
|---|
| 91 | diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) |
|---|
| 92 | enddo |
|---|
| 93 | if (diff_rho(1) > 0) then ! ice is at the surface |
|---|
| 94 | ice_table_beg(ig,islope) = 0. |
|---|
| [3321] | 95 | do isoilend = 2,nsoil_PEM - 1 |
|---|
| 96 | if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then |
|---|
| [3320] | 97 | call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) |
|---|
| 98 | ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) |
|---|
| 99 | exit |
|---|
| 100 | endif |
|---|
| 101 | enddo |
|---|
| 102 | else |
|---|
| 103 | 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 |
|---|
| 104 | if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then |
|---|
| 105 | call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope)) |
|---|
| 106 | ! Now let's find the end of the ice table |
|---|
| 107 | ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) |
|---|
| [3321] | 108 | do isoilend = isoil + 1,nsoil_PEM - 1 |
|---|
| [3320] | 109 | if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then |
|---|
| 110 | call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) |
|---|
| 111 | ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) |
|---|
| 112 | exit |
|---|
| 113 | endif |
|---|
| 114 | enddo |
|---|
| [3321] | 115 | endif ! diff_rho(z) <0 & diff_rho(z+1) > 0 |
|---|
| [3320] | 116 | enddo |
|---|
| 117 | endif ! diff_rho(1) > 0 |
|---|
| [2888] | 118 | enddo |
|---|
| [3320] | 119 | endif ! watercaptag |
|---|
| 120 | enddo |
|---|
| [2961] | 121 | |
|---|
| 122 | ! Small trick to speed up the convergence, Oded's idea. |
|---|
| [3320] | 123 | do islope = 1,nslope |
|---|
| 124 | do ig = 1,ngrid |
|---|
| 125 | if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then |
|---|
| [2961] | 126 | call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia) |
|---|
| 127 | stretch = (regolith_inertia(ig,islope)/wice_inertia)**2 |
|---|
| 128 | ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - & |
|---|
| [3320] | 129 | previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) |
|---|
| 130 | ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch |
|---|
| 131 | endif |
|---|
| 132 | enddo |
|---|
| 133 | enddo |
|---|
| [2961] | 134 | |
|---|
| [3321] | 135 | END SUBROUTINE computeice_table_equilibrium |
|---|
| [2888] | 136 | |
|---|
| [3321] | 137 | !----------------------------------------------------------------------- |
|---|
| [3320] | 138 | SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o) |
|---|
| [3031] | 139 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 140 | !!! |
|---|
| [3320] | 141 | !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed |
|---|
| 142 | !!! |
|---|
| [3031] | 143 | !!! Author: LL |
|---|
| 144 | !!! |
|---|
| 145 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| [3320] | 146 | use comsoil_h_PEM, only: mlayer_PEM |
|---|
| 147 | use comslope_mod, only: subslope_dist, def_slope_mean |
|---|
| 148 | use constants_marspem_mod, only: porosity |
|---|
| [3082] | 149 | #ifndef CPP_STD |
|---|
| 150 | use comcstfi_h, only: pi |
|---|
| 151 | #else |
|---|
| 152 | use comcstfi_mod, only: pi |
|---|
| 153 | #endif |
|---|
| 154 | |
|---|
| [3320] | 155 | implicit none |
|---|
| [3031] | 156 | |
|---|
| [3327] | 157 | ! Inputs |
|---|
| 158 | ! ------ |
|---|
| [3320] | 159 | integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM |
|---|
| 160 | real, dimension(ngrid,nslope), intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m] |
|---|
| 161 | real, dimension(ngrid,nslope), intent(in) :: new_ice_table_thickness ! ice table thickness at the current iteration [m] |
|---|
| 162 | real, dimension(ngrid,nslope), intent(in) :: ice_table_depth ! ice table depth [m] |
|---|
| 163 | real, dimension(ngrid,nslope), intent(in) :: tsurf ! Surface temperature [K] |
|---|
| 164 | real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil ! Soil temperature [K] |
|---|
| [3327] | 165 | ! Outputs |
|---|
| 166 | ! ------- |
|---|
| [3320] | 167 | 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] |
|---|
| [3327] | 168 | ! Locals |
|---|
| 169 | !------- |
|---|
| [3321] | 170 | integer :: ig, islope, ilay, iref ! loop index |
|---|
| 171 | real, dimension(ngrid,nslope) :: rho ! density of water ice [kg/m^3] |
|---|
| 172 | real, dimension(ngrid,nslope) :: Tice ! ice temperature [k] |
|---|
| [3327] | 173 | ! Code |
|---|
| 174 | ! ---- |
|---|
| [3320] | 175 | rho = 0. |
|---|
| 176 | Tice = 0. |
|---|
| [3031] | 177 | !1. First let's compute Tice using a linear interpolation between the mlayer level |
|---|
| [3320] | 178 | do ig = 1,ngrid |
|---|
| 179 | do islope = 1,nslope |
|---|
| [3327] | 180 | call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope)) |
|---|
| [3320] | 181 | rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 |
|---|
| 182 | enddo |
|---|
| 183 | enddo |
|---|
| [3031] | 184 | |
|---|
| 185 | !2. Let's compute the amount of ice that has sublimated in each subslope |
|---|
| [3320] | 186 | do ig = 1,ngrid |
|---|
| 187 | do islope = 1,nslope |
|---|
| [3321] | 188 | delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses |
|---|
| [3320] | 189 | *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) |
|---|
| 190 | enddo |
|---|
| 191 | enddo |
|---|
| [3031] | 192 | |
|---|
| [3321] | 193 | END SUBROUTINE compute_massh2o_exchange_ssi |
|---|
| [3031] | 194 | |
|---|
| [3321] | 195 | !----------------------------------------------------------------------- |
|---|
| [3320] | 196 | 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) |
|---|
| [2902] | 197 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 198 | !!! |
|---|
| [3320] | 199 | !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM |
|---|
| 200 | !!! |
|---|
| [2902] | 201 | !!! Author: LL |
|---|
| 202 | !!! |
|---|
| 203 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| [3320] | 204 | use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM |
|---|
| 205 | use math_mod, only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple |
|---|
| 206 | |
|---|
| [2902] | 207 | implicit none |
|---|
| [3320] | 208 | |
|---|
| 209 | ! Inputs |
|---|
| [2902] | 210 | ! ------ |
|---|
| [3320] | 211 | real, dimension(nsoilmx_PEM), intent(in) :: porefill ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice |
|---|
| 212 | real, dimension(nsoilmx_PEM), intent(in) :: psat_soil ! Soil water pressure at saturation, yearly averaged [Pa] |
|---|
| 213 | real, intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] |
|---|
| 214 | real, intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] |
|---|
| 215 | real, intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] |
|---|
| 216 | real, intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) |
|---|
| 217 | integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] |
|---|
| [2902] | 218 | real, intent(inout) :: depth_filling ! depth where pore filling begins [m] |
|---|
| [3320] | 219 | ! Outputs |
|---|
| [2902] | 220 | ! ------- |
|---|
| [3320] | 221 | integer, intent(out) :: index_filling ! index where the pore filling begins [1] |
|---|
| 222 | integer, intent(out) :: index_geothermal ! index where the ice table stops [1] |
|---|
| 223 | real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] |
|---|
| 224 | 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 |
|---|
| [3327] | 225 | ! Locals |
|---|
| 226 | !------- |
|---|
| [3320] | 227 | real, dimension(nsoilmx_PEM) :: eta ! constriction |
|---|
| 228 | real, dimension(nsoilmx_PEM) :: dz_psat ! first derivative of the vapor pressure at saturation |
|---|
| 229 | real, dimension(nsoilmx_PEM) :: dz_eta ! \partial z \eta |
|---|
| 230 | real, dimension(nsoilmx_PEM) :: dzz_psat ! \partial \partial psat |
|---|
| 231 | integer :: ilay, index_tmp ! index for loop |
|---|
| 232 | real :: old_depth_filling ! depth_filling saved |
|---|
| 233 | real :: Jdry ! flux trought the dry layer |
|---|
| 234 | real :: Jsat ! flux trought the ice layer |
|---|
| 235 | real :: Jdry_prevlay, Jsat_prevlay ! same but for the previous ice layer |
|---|
| 236 | integer :: index_firstice ! first index where ice appears (i.e., f > 0) |
|---|
| 237 | real :: massfillabove, massfillafter ! h2O mass above and after index_geothermal |
|---|
| [3327] | 238 | ! Constant |
|---|
| 239 | !--------- |
|---|
| [3320] | 240 | real, parameter :: pvap2rho = 18.e-3/8.314 |
|---|
| [3327] | 241 | ! Code |
|---|
| 242 | !----- |
|---|
| [2902] | 243 | ! 0. Compute constriction over the layer |
|---|
| 244 | ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 |
|---|
| [3321] | 245 | if (index_IS < 0) then |
|---|
| [3320] | 246 | index_tmp = nsoilmx_PEM |
|---|
| 247 | do ilay = 1,nsoilmx_PEM |
|---|
| 248 | call constriction(porefill(ilay),eta(ilay)) |
|---|
| 249 | enddo |
|---|
| [2902] | 250 | else |
|---|
| [3320] | 251 | index_tmp = index_IS |
|---|
| [3321] | 252 | do ilay = 1,index_IS - 1 |
|---|
| [3320] | 253 | call constriction(porefill(ilay),eta(ilay)) |
|---|
| 254 | enddo |
|---|
| 255 | do ilay = index_IS,nsoilmx_PEM |
|---|
| 256 | eta(ilay) = 0. |
|---|
| 257 | enddo |
|---|
| [2902] | 258 | endif |
|---|
| 259 | |
|---|
| [3320] | 260 | ! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| [2902] | 261 | old_depth_filling = depth_filling |
|---|
| 262 | |
|---|
| [3320] | 263 | call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) |
|---|
| [2902] | 264 | |
|---|
| 265 | do ilay = 1,index_tmp |
|---|
| [3320] | 266 | Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| 267 | Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) |
|---|
| 268 | if (Jdry - Jsat <= 0) then |
|---|
| 269 | index_filling = ilay |
|---|
| 270 | exit |
|---|
| 271 | endif |
|---|
| [2902] | 272 | enddo |
|---|
| 273 | |
|---|
| [3320] | 274 | if (index_filling == 1) then |
|---|
| 275 | depth_filling = mlayer_PEM(1) |
|---|
| 276 | else if (index_filling > 1) then |
|---|
| 277 | Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1) |
|---|
| 278 | Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1) |
|---|
| 279 | call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling) |
|---|
| [2902] | 280 | endif |
|---|
| 281 | |
|---|
| 282 | ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010)) |
|---|
| 283 | ! 2.0 preliminary: depth to shallowest ice (discontinuity at interface) |
|---|
| 284 | index_firstice = -1 |
|---|
| 285 | do ilay = 1,nsoilmx_PEM |
|---|
| [3320] | 286 | if (porefill(ilay) <= 0.) then |
|---|
| [2902] | 287 | index_firstice = ilay ! first point with ice |
|---|
| 288 | exit |
|---|
| 289 | endif |
|---|
| 290 | enddo |
|---|
| 291 | |
|---|
| 292 | ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) |
|---|
| [3320] | 293 | call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta) |
|---|
| [2902] | 294 | |
|---|
| [3320] | 295 | if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) |
|---|
| [2902] | 296 | |
|---|
| 297 | call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) |
|---|
| [3320] | 298 | dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat) |
|---|
| [2902] | 299 | |
|---|
| 300 | ! 3. Ice table boundary due to geothermal heating |
|---|
| [3320] | 301 | if (index_IS > 0) index_geothermal = -1 |
|---|
| 302 | if (index_geothermal < 0) depth_geothermal = -1. |
|---|
| 303 | if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010) |
|---|
| 304 | index_geothermal = -1 |
|---|
| 305 | do ilay = 2,nsoilmx_PEM |
|---|
| 306 | if (dz_psat(ilay) > 0.) then ! first point with reversed flux |
|---|
| 307 | index_geothermal = ilay |
|---|
| 308 | call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) |
|---|
| 309 | exit |
|---|
| [2902] | 310 | endif |
|---|
| [3320] | 311 | enddo |
|---|
| 312 | else |
|---|
| 313 | index_geothermal = -1 |
|---|
| 314 | endif |
|---|
| 315 | if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010) |
|---|
| [3321] | 316 | call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,index_geothermal - 1,nsoilmx_PEM,massfillabove) |
|---|
| [3320] | 317 | index_tmp = -1 |
|---|
| 318 | do ilay = index_geothermal,nsoilmx_PEM |
|---|
| 319 | if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full |
|---|
| [3321] | 320 | call colint(porefill/eta,mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) |
|---|
| [3320] | 321 | if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG |
|---|
| 322 | if (ilay > index_geothermal) then |
|---|
| 323 | ! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal |
|---|
| 324 | call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, & |
|---|
| 325 | dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) |
|---|
| 326 | ! 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) |
|---|
| 327 | index_tmp = ilay |
|---|
| [2902] | 328 | endif |
|---|
| 329 | ! otherwise leave depth_geothermal unchanged |
|---|
| 330 | exit |
|---|
| 331 | endif |
|---|
| 332 | massfillabove = massfillafter |
|---|
| [3320] | 333 | enddo |
|---|
| 334 | if (index_tmp > 0) index_geothermal = index_tmp |
|---|
| 335 | end if |
|---|
| [2902] | 336 | |
|---|
| [3321] | 337 | END SUBROUTINE find_layering_icetable |
|---|
| [2902] | 338 | |
|---|
| [3321] | 339 | !----------------------------------------------------------------------- |
|---|
| [2902] | 340 | SUBROUTINE constriction(porefill,eta) |
|---|
| 341 | |
|---|
| [3320] | 342 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 343 | !!! |
|---|
| 344 | !!! Purpose: Compute the constriction of vapor flux by pore ice |
|---|
| 345 | !!! |
|---|
| 346 | !!! Author: LL |
|---|
| 347 | !!! |
|---|
| 348 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| [2902] | 349 | |
|---|
| [3320] | 350 | implicit none |
|---|
| [2902] | 351 | |
|---|
| [3320] | 352 | real, intent(in) :: porefill ! pore filling fraction |
|---|
| 353 | real, intent(out) :: eta ! constriction |
|---|
| [2902] | 354 | |
|---|
| [3320] | 355 | if (porefill <= 0.) then |
|---|
| 356 | eta = 1. |
|---|
| 357 | else if (0 < porefill .and. porefill < 1.) then |
|---|
| 358 | eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 |
|---|
| 359 | else |
|---|
| 360 | eta = 0. |
|---|
| 361 | endif |
|---|
| 362 | |
|---|
| [3321] | 363 | END SUBROUTINE constriction |
|---|
| [3320] | 364 | |
|---|
| [3327] | 365 | !----------------------------------------------------------------------- |
|---|
| 366 | SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice) |
|---|
| 367 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 368 | !!! |
|---|
| 369 | !!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells. |
|---|
| 370 | !!! |
|---|
| 371 | !!! Author: LL |
|---|
| 372 | !!! |
|---|
| 373 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 374 | |
|---|
| 375 | use comsoil_h_PEM, only: layer_PEM, mlayer_PEM |
|---|
| 376 | |
|---|
| 377 | implicit none |
|---|
| 378 | |
|---|
| 379 | ! Inputs |
|---|
| 380 | ! ------ |
|---|
| 381 | integer, intent(in) :: nsoil ! Number of soil layers |
|---|
| 382 | real, dimension(nsoil), intent(in) :: ptsoil ! Soil temperature (K) |
|---|
| 383 | real, intent(in) :: ptsurf ! Soil temperature (K) |
|---|
| 384 | real, intent(in) :: ice_depth ! Ice depth (m) |
|---|
| 385 | |
|---|
| 386 | ! Outputs |
|---|
| 387 | ! ------ |
|---|
| 388 | real, intent(out) :: Tice ! Ice temperatures (K) |
|---|
| 389 | |
|---|
| 390 | ! Locals |
|---|
| 391 | ! ------ |
|---|
| 392 | integer :: ik ! Loop variables |
|---|
| 393 | integer :: indexice ! Index of the ice |
|---|
| 394 | |
|---|
| 395 | ! Code |
|---|
| 396 | !----- |
|---|
| 397 | indexice = -1 |
|---|
| 398 | if (ice_depth >= mlayer_PEM(nsoil - 1)) then |
|---|
| 399 | Tice = ptsoil(nsoil) |
|---|
| 400 | else |
|---|
| 401 | if(ice_depth < mlayer_PEM(0)) then |
|---|
| 402 | indexice = 0. |
|---|
| 403 | else |
|---|
| 404 | do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations |
|---|
| 405 | if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then |
|---|
| 406 | indexice = ik + 1 |
|---|
| 407 | exit |
|---|
| 408 | endif |
|---|
| 409 | enddo |
|---|
| 410 | endif |
|---|
| 411 | if (indexice < 0) then |
|---|
| 412 | call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1) |
|---|
| 413 | else |
|---|
| 414 | if(indexice >= 1) then ! Linear inteprolation between soil temperature |
|---|
| 415 | Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1) |
|---|
| 416 | else ! Linear inteprolation between the 1st soil temperature and the surface temperature |
|---|
| 417 | Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1) |
|---|
| 418 | endif ! index ice >= 1 |
|---|
| 419 | endif !indexice < 0 |
|---|
| 420 | endif ! icedepth > mlayer_PEM(nsoil - 1) |
|---|
| 421 | |
|---|
| 422 | END SUBROUTINE compute_Tice_pem |
|---|
| 423 | |
|---|
| [3321] | 424 | END MODULE ice_table_mod |
|---|
