Changeset 2902 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- Feb 22, 2023, 6:56:01 PM (22 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r2888 r2902 28 28 29 29 #ifndef CPP_STD 30 USE comsoil_h_PEM, only: layer_PEM ! Depth of the vertical grid30 USE comsoil_h_PEM, only: mlayer_PEM ! Depth of the vertical grid 31 31 implicit none 32 32 … … 46 46 else 47 47 do islope = 1,nslope 48 ice_table(ig,islope) = -1 0.48 ice_table(ig,islope) = -1. 49 49 50 50 do isoil = 1,nsoil_PEM … … 57 57 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 58 58 if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then 59 z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1)) 60 z2 = -layer_PEM(isoil+1)*z1 + diff_rho(isoil+1) 61 ice_table(ig,islope) = -z2/z1 59 call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope)) 62 60 exit 63 61 endif !diff_rho(z) <0 & diff_rho(z+1) > 0 … … 72 70 END 73 71 74 75 76 end module 72 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) 73 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 74 !!! 75 !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM 76 !!! 77 !!! Author: LL 78 !!! 79 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 80 use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM 81 implicit none 82 ! inputs 83 ! ------ 84 85 real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice 86 real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa] 87 real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] 88 real,intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] 89 real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] 90 real,intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) 91 integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] 92 real, intent(inout) :: depth_filling ! depth where pore filling begins [m] 93 94 ! outputs 95 ! ------- 96 integer,intent(out) :: index_filling ! index where the pore filling begins [1] 97 integer, intent(out) :: index_geothermal ! index where the ice table stops [1] 98 real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] 99 real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase 100 101 ! local 102 103 real :: eta(nsoilmx_PEM) ! constriction 104 integer :: ilay ! index for loop 105 real :: old_depth_filling ! depth_filling saved 106 real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn 107 integer :: index_tmp ! for loop 108 real :: Jdry ! flux trought the dry layer 109 real :: Jsat ! flux trought the ice layer 110 real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer 111 integer :: index_firstice ! first index where ice appears (i.e., f > 0) 112 real :: dz_eta(nsoilmx_PEM) ! \partial z \eta 113 real :: dz_eta_low ! same but evaluated at the interface for ice 114 real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat 115 real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal 116 117 ! constant 118 real :: pvap2rho = 18.e-3/8.314 119 ! 120 121 122 123 124 125 ! 0. Compute constriction over the layer 126 ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 127 if (index_IS.lt.0) then 128 index_tmp = nsoilmx_PEM 129 do ilay = 1,nsoilmx_PEM 130 call constriction(porefill(ilay),eta(ilay)) 131 enddo 132 else 133 index_tmp = index_IS 134 do ilay = 1,index_IS-1 135 call constriction(porefill(ilay),eta(ilay)) 136 enddo 137 do ilay = index_IS,nsoilmx_PEM 138 eta(ilay) = 0. 139 enddo 140 endif 141 142 ! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) 143 144 old_depth_filling = depth_filling 145 146 call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 147 148 do ilay = 1,index_tmp 149 Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) 150 Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) 151 if((Jdry - Jsat).le.0) then 152 index_filling = ilay 153 exit 154 endif 155 enddo 156 157 if(index_filling.eq.1) depth_filling = mlayer_PEM(1) 158 159 if(index_filling.gt.1) then 160 Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1) 161 Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1) 162 call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling) 163 endif 164 165 166 ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010)) 167 168 169 ! 2.0 preliminary: depth to shallowest ice (discontinuity at interface) 170 171 index_firstice = -1 172 do ilay = 1,nsoilmx_PEM 173 if (porefill(ilay).le.0.) then 174 index_firstice = ilay ! first point with ice 175 exit 176 endif 177 enddo 178 179 ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) 180 181 call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta) 182 183 if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then 184 call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) 185 endif 186 187 call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) 188 dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:)) 189 190 ! 3. Ice table boundary due to geothermal heating 191 192 if(index_IS.gt.0) index_geothermal = -1 193 if(index_geothermal.lt.0) depth_geothermal = -1. 194 if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010) 195 index_geothermal = -1 196 do ilay=2,nsoilmx_PEM 197 if (dz_psat(ilay).gt.0.) then ! first point with reversed flux 198 index_geothermal=ilay 199 call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) 200 exit 201 endif 202 enddo 203 else 204 index_geothermal = -1 205 endif 206 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 24 from Schorgofer, Icarus (2010) 207 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) 208 index_tmp = -1 209 do ilay=index_geothermal,nsoilmx_PEM 210 if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full 211 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) 212 if (massfillafter<dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG 213 if (ilay>index_geothermal) then 214 ! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal 215 call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 216 dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) 217 ! if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*) 218 ! '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay) 219 index_tmp=ilay 220 endif 221 ! otherwise leave depth_geothermal unchanged 222 exit 223 endif 224 massfillabove = massfillafter 225 enddo 226 if (index_tmp.gt.0) index_geothermal = index_tmp 227 end if 228 return 229 end subroutine 230 231 232 233 234 235 SUBROUTINE findroot(y1,y2,z1,z2,zr) 236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 237 !!! 238 !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2 239 !!! 240 !!! Author: LL 241 !!! 242 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 243 244 implicit none 245 real,intent(in) :: y1,y2 ! difference between surface water density and at depth [kg/m^3] 246 real,intent(in) :: z1,z2 ! depth [m] 247 real,intent(out) :: zr ! depth at which we have zero 248 zr = (y1*z2 - y2*z1)/(y1-y2) 249 RETURN 250 end 251 252 SUBROUTINE constriction(porefill,eta) 253 254 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 255 !!! 256 !!! Purpose: Compute the constriction of vapor flux by pore ice 257 !!! 258 !!! Author: LL 259 !!! 260 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 261 implicit none 262 real,intent(in) :: porefill ! pore filling fraction 263 real,intent(out) :: eta ! constriction 264 265 !!! 266 267 268 if (porefill.le.0.) eta = 1. 269 if ((porefill.gt.0.) .and.(porefill.lt.1.)) then 270 eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 271 endif 272 if (porefill.le.1.) eta = 0. 273 return 274 end 275 276 277 278 subroutine deriv1(z,nz,y,y0,ybot,dzY) 279 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 280 !!! 281 !!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 282 !!! 283 !!! Author: From N.S (N.S, Icarus 2010), impletented here by LL 284 !!! 285 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 286 implicit none 287 288 ! first derivative of a function y(z) on irregular grid 289 ! upper boundary conditions: y(0)=y0 290 ! lower boundary condition.: yp = ybottom 291 integer, intent(IN) :: nz ! number of layer 292 real, intent(IN) :: z(nz) ! depth layer 293 real, intent(IN) :: y(nz) ! function which needs to be derived 294 real, intent(IN) :: y0,ybot ! boundary conditions 295 real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth 296 ! local 297 integer :: j 298 real :: hm,hp,c1,c2,c3 299 300 hp = z(2)-z(1) 301 c1 = z(1)/(hp*z(2)) 302 c2 = 1/z(1) - 1/(z(2)-z(1)) 303 c3 = -hp/(z(1)*z(2)) 304 dzY(1) = c1*y(2) + c2*y(1) + c3*y0 305 do j=2,nz-1 306 hp = z(j+1)-z(j) 307 hm = z(j)-z(j-1) 308 c1 = +hm/(hp*(z(j+1)-z(j-1))) 309 c2 = 1/hm - 1/hp 310 c3 = -hp/(hm*(z(j+1)-z(j-1))) 311 dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) 312 enddo 313 dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1))) 314 return 315 end subroutine deriv1 316 317 318 319 subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2) 320 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 321 !!! 322 !!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 323 !!! 324 !!! Author: N.S (raw copy/paste from MSIM) 325 !!! 326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 327 328 ! second derivative y_zz on irregular grid 329 ! boundary conditions: y(0)=y0, y(nz+1)=yNp1 330 implicit none 331 integer, intent(IN) :: nz 332 real, intent(IN) :: z(nz),y(nz),y0,yNp1 333 real, intent(OUT) :: yp2(nz) 334 integer j 335 real hm,hp,c1,c2,c3 336 337 c1 = +2./((z(2)-z(1))*z(2)) 338 c2 = -2./((z(2)-z(1))*z(1)) 339 c3 = +2./(z(1)*z(2)) 340 yp2(1) = c1*y(2) + c2*y(1) + c3*y0 341 342 do j=2,nz-1 343 hp = z(j+1)-z(j) 344 hm = z(j)-z(j-1) 345 c1 = +2./(hp*(z(j+1)-z(j-1))) 346 c2 = -2./(hp*hm) 347 c3 = +2./(hm*(z(j+1)-z(j-1))) 348 yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) 349 enddo 350 yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2 351 return 352 end subroutine deriv2_simple 353 354 355 356 subroutine deriv1_onesided(j,z,nz,y,dy_zj) 357 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 358 !!! 359 !!! Purpose: First derivative of function y(z) at z(j) one-sided derivative on irregular grid 360 !!! 361 !!! Author: N.S (raw copy/paste from MSIM) 362 !!! 363 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 364 365 implicit none 366 integer, intent(IN) :: nz,j 367 real, intent(IN) :: z(nz),y(nz) 368 real, intent(out) :: dy_zj 369 real h1,h2,c1,c2,c3 370 if (j<1 .or. j>nz-2) then 371 dy_zj = -1. 372 else 373 h1 = z(j+1)-z(j) 374 h2 = z(j+2)-z(j+1) 375 c1 = -(2*h1+h2)/(h1*(h1+h2)) 376 c2 = (h1+h2)/(h1*h2) 377 c3 = -h1/(h2*(h1+h2)) 378 dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2) 379 endif 380 return 381 end subroutine deriv1_onesided 382 383 384 subroutine colint(y,z,nz,i1,i2,integral) 385 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 386 !!! 387 !!! Purpose: Column integrates y on irregular grid 388 !!! 389 !!! Author: N.S (raw copy/paste from MSIM) 390 !!! 391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 392 393 implicit none 394 integer, intent(IN) :: nz, i1, i2 395 real, intent(IN) :: y(nz), z(nz) 396 real,intent(out) :: integral 397 integer i 398 real dz(nz) 399 400 dz(1) = (z(2)-0.)/2 401 do i=2,nz-1 402 dz(i) = (z(i+1)-z(i-1))/2. 403 enddo 404 dz(nz) = z(nz)-z(nz-1) 405 integral = sum(y(i1:i2)*dz(i1:i2)) 406 end subroutine colint 407 408 409 410 411 412 413 end module
Note: See TracChangeset
for help on using the changeset viewer.
