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