1 !----------------------------------- 2 subroutine sflx & 3 !................................... 4 ! --- inputs: 5 & ( nsoil, couple, icein, ffrozp, dt, zlvl, sldpth, & 6 & swdn, swnet, lwdn, sfcems, sfcprs, sfctmp, & 7 & sfcspd, prcp, q2, q2sat, dqsdt2, th2, & 8 & vegtyp, soiltyp, slopetyp, shdmin, alb, snoalb, & 9 ! --- input/outputs: 10 & tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm, & 11 ! --- outputs: 12 & nroot, shdfac, snowh, albedo, eta, sheat, ec, & 13 & edir, et, ett, esnow, drip, dew, beta, etp, ssoil, & 14 & flx1, flx2, flx3, runoff1, runoff2, runoff3, & 15 & snomlt, sncovr, rc, pc, rsmin, xlai, rcs, rct, rcq, & 16 & rcsoil, soilw, soilm, smcwlt, smcdry, smcref, smcmax & 17 & ) 18 19 ! ===================================================================== ! 20 ! description: ! 21 ! ! 22 ! subroutine sflx - version 2.7: ! 23 ! sub-driver for "noah/osu lsm" family of physics subroutines for a ! 24 ! soil/veg/snowpack land-surface model to update soil moisture, soil ! 25 ! ice, soil temperature, skin temperature, snowpack water content, ! 26 ! snowdepth, and all terms of the surface energy balance and surface ! 27 ! water balance (excluding input atmospheric forcings of downward ! 28 ! radiation and precip) ! 29 ! ! 30 ! usage: ! 31 ! ! 32 ! call sflx ! 33 ! --- inputs: ! 34 ! ( nsoil, couple, icein, ffrozp, dt, zlvl, sldpth, ! 35 ! swdn, swnet, lwdn, sfcems, sfcprs, sfctmp, ! 36 ! sfcspd, prcp, q2, q2sat, dqsdt2, th2, ! 37 ! vegtyp, soiltyp, slopetyp, shdmin, alb, snoalb, ! 38 ! --- input/outputs: ! 39 ! tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm, ! 40 ! --- outputs: ! 41 ! nroot, shdfac, snowh, albedo, eta, sheat, ec, ! 42 ! edir, et, ett, esnow, drip, dew, beta, etp, ssoil, ! 43 ! flx1, flx2, flx3, runoff1, runoff2, runoff3, ! 44 ! snomlt, sncovr, rc, pc, rsmin, xlai, rcs, rct, rcq, ! 45 ! rcsoil, soilw, soilm, smcwlt, smcdry, smcref, smcmax ) ! 46 ! ! 47 ! ! 48 ! subprograms called: redprm, snow_new, csnow, snfrac, alcalc, ! 49 ! tdfcnd, snowz0, sfcdif, penman, canres, nopac, snopac. ! 50 ! ! 51 ! ! 52 ! program history log: ! 53 ! jun 2003 -- K. Mitchell et. al -- created version 2.7 ! 54 ! 200x -- sarah lu modified the code including: ! 55 ! added passing argument, couple; replaced soldn ! 56 ! and solnet by radflx; call sfcdif if couple=0; ! 57 ! apply time filter to stc and tskin; and the ! 58 ! way of namelist inport. ! 59 ! feb 2004 -- m. ek noah v2.7.1 non-linear weighting of snow vs ! 60 ! non-snow covered portions of gridbox ! 61 ! apr 2009 -- y.-t. hou added lw surface emissivity effect, ! 62 ! streamlined and reformatted the code, and ! 63 ! consolidated constents/parameters by using ! 64 ! module physcons, and added program documentation! ! 65 ! sep 2009 -- s. moorthi minor fixes 66 ! ! 67 ! ==================== defination of variables ==================== ! 68 ! ! 69 ! inputs: size ! 70 ! nsoil - integer, number of soil layers (>=2 but <=nsold) 1 ! 71 ! couple - integer, =0:uncoupled (land model only) 1 ! 72 ! =1:coupled with parent atmos model ! 73 ! icein - integer, sea-ice flag (=1: sea-ice, =0: land) 1 ! 74 ! ffrozp - real, 1 ! 75 ! dt - real, time step (<3600 sec) 1 ! 76 ! zlvl - real, height abv atmos ground forcing vars (m) 1 ! 77 ! sldpth - real, thickness of each soil layer (m) nsoil ! 78 ! swdn - real, downward sw radiation flux (w/m**2) 1 ! 79 ! swnet - real, downward sw net (dn-up) flux (w/m**2) 1 ! 80 ! lwdn - real, downward lw radiation flux (w/m**2) 1 ! 81 ! sfcems - real, sfc lw emissivity (fractional) 1 ! 82 ! sfcprs - real, pressure at height zlvl abv ground(pascals) 1 ! 83 ! sfctmp - real, air temp at height zlvl abv ground (k) 1 ! 84 ! sfcspd - real, wind speed at height zlvl abv ground (m/s) 1 ! 85 ! prcp - real, precip rate (kg m-2 s-1) 1 ! 86 ! q2 - real, mixing ratio at hght zlvl abv grnd (kg/kg) 1 ! 87 ! q2sat - real, sat mixing ratio at zlvl abv grnd (kg/kg) 1 ! 88 ! dqsdt2 - real, slope of sat specific humidity curve at 1 ! 89 ! t=sfctmp (kg kg-1 k-1) ! 90 ! th2 - real, air potential temp at zlvl abv grnd (k) 1 ! 91 ! vegtyp - integer, vegetation type (integer index) 1 ! 92 ! soiltyp - integer, soil type (integer index) 1 ! 93 ! slopetyp - integer, class of sfc slope (integer index) 1 ! 94 ! shdmin - real, min areal coverage of green veg (fraction) 1 ! 95 ! alb - real, bkground snow-free sfc albedo (fraction) 1 ! 96 ! snoalb - real, max albedo over deep snow (fraction) 1 ! 97 ! ! 98 ! input/outputs: ! 99 ! tbot - real, bottom soil temp (k) 1 ! 100 ! (local yearly-mean sfc air temp) ! 101 ! cmc - real, canopy moisture content (m) 1 ! 102 ! t1 - real, ground/canopy/snowpack eff skin temp (k) 1 ! 103 ! stc - real, soil temp (k) nsoil ! 104 ! smc - real, total soil moisture (vol fraction) nsoil ! 105 ! sh2o - real, unfrozen soil moisture (vol fraction) nsoil ! 106 ! note: frozen part = smc-sh2o ! 107 ! sneqv - real, water-equivalent snow depth (m) 1 ! 108 ! note: snow density = snwqv/snowh ! 109 ! ch - real, sfc exchange coeff for heat & moisture (m/s)1 ! 110 ! note: conductance since it's been mult by wind ! 111 ! cm - real, sfc exchange coeff for momentum (m/s) 1 ! 112 ! note: conductance since it's been mult by wind ! 113 ! ! 114 ! outputs: ! 115 ! nroot - integer, number of root layers 1 ! 116 ! shdfac - real, aeral coverage of green veg (fraction) 1 ! 117 ! snowh - real, snow depth (m) 1 ! 118 ! albedo - real, sfc albedo incl snow effect (fraction) 1 ! 119 ! eta - real, downward latent heat flux (w/m2) 1 ! 120 ! sheat - real, downward sensible heat flux (w/m2) 1 ! 121 ! ec - real, canopy water evaporation (w/m2) 1 ! 122 ! edir - real, direct soil evaporation (w/m2) 1 ! 123 ! et - real, plant transpiration (w/m2) nsoil ! 124 ! ett - real, total plant transpiration (w/m2) 1 ! 125 ! esnow - real, sublimation from snowpack (w/m2) 1 ! 126 ! drip - real, through-fall of precip and/or dew in excess 1 ! 127 ! of canopy water-holding capacity (m) ! 128 ! dew - real, dewfall (or frostfall for t<273.15) (m) 1 ! 129 ! beta - real, ratio of actual/potential evap 1 ! 130 ! etp - real, potential evaporation (w/m2) 1 ! 131 ! ssoil - real, upward soil heat flux (w/m2) 1 ! 132 ! flx1 - real, precip-snow sfc flux (w/m2) 1 ! 133 ! flx2 - real, freezing rain latent heat flux (w/m2) 1 ! 134 ! flx3 - real, phase-change heat flux from snowmelt (w/m2) 1 ! 135 ! snomlt - real, snow melt (m) (water equivalent) 1 ! 136 ! sncovr - real, fractional snow cover 1 ! 137 ! runoff1 - real, surface runoff (m/s) not infiltrating sfc 1 ! 138 ! runoff2 - real, sub sfc runoff (m/s) (baseflow) 1 ! 139 ! runoff3 - real, excess of porosity for a given soil layer 1 ! 140 ! rc - real, canopy resistance (s/m) 1 ! 141 ! pc - real, plant coeff (fraction) where pc*etp=transpi 1 ! 142 ! rsmin - real, minimum canopy resistance (s/m) 1 ! 143 ! xlai - real, leaf area index (dimensionless) 1 ! 144 ! rcs - real, incoming solar rc factor (dimensionless) 1 ! 145 ! rct - real, air temp rc factor (dimensionless) 1 ! 146 ! rcq - real, atoms vapor press deficit rc factor 1 ! 147 ! rcsoil - real, soil moisture rc factor (dimensionless) 1 ! 148 ! soilw - real, available soil mois in root zone 1 ! 149 ! soilm - real, total soil column mois (frozen+unfrozen) (m)1 ! 150 ! smcwlt - real, wilting point (volumetric) 1 ! 151 ! smcdry - real, dry soil mois threshold (volumetric) 1 ! 152 ! smcref - real, soil mois threshold (volumetric) 1 ! 153 ! smcmax - real, porosity (sat val of soil mois) 1 ! 154 ! ! 155 ! ==================== end of description ===================== ! 156 ! 157 use machine , only : kind_phys 158 ! 159 use physcons, only : con_cp, con_rd, con_t0c, con_g, con_pi, & 160 & con_cliq, con_csol, con_hvap, con_hfus, & 161 & con_sbc 162 ! 163 implicit none 164 165 ! --- constant parameters: 166 ! *** note: some of the constants are different in subprograms and need to 167 ! be consolidated with the standard def in module physcons at sometime 168 ! at the present time, those diverse values are kept temperately to 169 ! provide the same result as the original codes. -- y.t.h. may09 170 171 integer, parameter :: nsold = 4 ! max soil layers 172 173 real (kind=kind_phys), parameter :: gs = con_g ! con_g =9.80665 174 real (kind=kind_phys), parameter :: gs1 = 9.8 ! con_g in sfcdif 175 real (kind=kind_phys), parameter :: gs2 = 9.81 ! con_g in snowpack, frh2o 176 real (kind=kind_phys), parameter :: tfreez = con_t0c ! con_t0c =275.15 177 real (kind=kind_phys), parameter :: lsubc = 2.501e+6 ! con_hvap=2.5000e+6 178 real (kind=kind_phys), parameter :: lsubf = 3.335e5 ! con_hfus=3.3358e+5 179 real (kind=kind_phys), parameter :: lsubs = 2.83e+6 ! ? in sflx, snopac 180 real (kind=kind_phys), parameter :: elcp = 2.4888e+3 ! ? in penman 181 real (kind=kind_phys), parameter :: rd = con_rd ! con_rd =287.05 182 real (kind=kind_phys), parameter :: rd1 = 287.04 ! con_rd in sflx, penman, canres 183 real (kind=kind_phys), parameter :: cp = con_cp ! con_cp =1004.6 184 real (kind=kind_phys), parameter :: cp1 = 1004.5 ! con_cp in sflx, canres 185 real (kind=kind_phys), parameter :: cp2 = 1004.0 ! con_cp in htr 186 real (kind=kind_phys), parameter :: cph2o = con_cliq ! con_cliq=4.1855e+3 187 real (kind=kind_phys), parameter :: cph2o1 = 4.218e+3 ! con_cliq in penman, snopac 188 real (kind=kind_phys), parameter :: cph2o2 = 4.2e6 ! con_cliq in hrt *unit diff! 189 real (kind=kind_phys), parameter :: cpice = con_csol ! con_csol=2.106e+3 190 real (kind=kind_phys), parameter :: cpice1 = 2.106e6 ! con_csol in hrt *unit diff! 191 real (kind=kind_phys), parameter :: sigma = con_sbc ! con_sbc=5.6704e-8 192 real (kind=kind_phys), parameter :: sigma1 = 5.67e-8 ! con_sbc in penman, nopac, snopac 193 194 ! --- inputs: 195 integer, intent(in) :: nsoil, couple, icein, vegtyp, soiltyp, & 196 & slopetyp 197 198 real (kind=kind_phys), intent(in) :: ffrozp, dt, zlvl, lwdn, & 199 & sldpth(nsoil), swdn, swnet, sfcems, sfcprs, sfctmp, & 200 & sfcspd, prcp, q2, q2sat, dqsdt2, th2, shdmin, alb, snoalb 201 202 ! --- input/outputs: 203 real (kind=kind_phys), intent(inout) :: tbot, cmc, t1, sneqv, & 204 & stc(nsoil), smc(nsoil), sh2o(nsoil), ch, cm 205 206 ! --- outputs: 207 integer, intent(out) :: nroot 208 209 real (kind=kind_phys), intent(out) :: shdfac, snowh, albedo, & 210 & eta, sheat, ec, edir, et(nsoil), ett, esnow, drip, dew, & 211 & beta, etp, ssoil, flx1, flx2, flx3, snomlt, sncovr, & 212 & runoff1, runoff2, runoff3, rc, pc, rsmin, xlai, rcs, & 213 & rct, rcq, rcsoil, soilw, soilm, smcwlt, smcdry, smcref, & 214 & smcmax 215 216 ! --- locals: 217 real (kind=kind_phys) :: bexp, cfactr, cmcmax, csoil, czil, & 218 & df1, df1h, df1a, dksat, dwsat, dsoil, dtot, frcsno, & 219 & frcsoi, epsca, fdown, f1, fxexp, frzx, hs, kdt, prcp1, & 220 & psisat, quartz, rch, refkdt, rr, rgl, rsmax, sndens, & 221 & sncond, sbeta, sn_new, slope, snup, salp, soilwm, soilww, & 222 & t1v, t24, t2v, th2v, topt, tsnow, zbot, z0 223 224 real (kind=kind_phys), dimension(nsold) :: rtdis, zsoil 225 226 logical :: frzgra, snowng 227 228 integer :: ice, k, kz 229 230 ! 231 !===> ... begin here 232 ! 233 ! --- ... initialization 234 235 runoff1 = 0.0 236 runoff2 = 0.0 237 runoff3 = 0.0 238 snomlt = 0.0 239 240 ! --- ... define local variable ice to achieve: 241 ! sea-ice case, ice = 1 242 ! non-glacial land, ice = 0 243 ! glacial-ice land, ice = -1 244 ! if vegtype=13 (glacial-ice), re-set ice flag = -1 (glacial-ice) 245 ! note - for open-sea, sflx should *not* have been called. set green 246 ! vegetation fraction (shdfac) = 0. 247 248 ice = icein 249 250 if (vegtyp == 13) then 251 ice = -1 252 shdfac = 0.0 253 endif 254 255 if (ice == 1) then 256 257 shdfac = 0.0 258 259 ! --- ... set green vegetation fraction (shdfac) = 0. 260 ! set sea-ice layers of equal thickness and sum to 3 meters 261 262 do kz = 1, nsoil 263 zsoil(kz) = -3.0 * float(kz) / float(nsoil) 264 enddo 265 266 else 267 268 ! --- ... calculate depth (negative) below ground from top skin sfc to 269 ! bottom of each soil layer. 270 ! note - sign of zsoil is negative (denoting below ground) 271 272 zsoil(1) = -sldpth(1) 273 do kz = 2, nsoil 274 zsoil(kz) = -sldpth(kz) + zsoil(kz-1) 275 end do 276 277 endif ! end if_ice_block 278 279 ! --- ... next is crucial call to set the land-surface parameters, 280 ! including soil-type and veg-type dependent parameters. 281 ! set shdfac=0.0 for bare soil surfaces 282 283 call redprm 284 ! --- inputs: ! 285 ! ( nsoil, vegtyp, soiltyp, slopetyp, sldpth, zsoil, ! 286 ! --- outputs: ! 287 ! cfactr, cmcmax, rsmin, rsmax, topt, refkdt, kdt, ! 288 ! sbeta, shdfac, rgl, hs, zbot, frzx, psisat, slope, ! 289 ! snup, salp, bexp, dksat, dwsat, smcmax, smcwlt, ! 290 ! smcref, smcdry, f1, quartz, fxexp, rtdis, nroot, ! 291 ! z0, czil, xlai, csoil ) ! 292 293 ! --- ... initialize precipitation logicals. 294 295 snowng = .false. 296 frzgra = .false. 297 298 ! --- ... over sea-ice or glacial-ice, if s.w.e. (sneqv) below threshold 299 ! lower bound (0.01 m for sea-ice, 0.10 m for glacial-ice), then 300 ! set at lower bound and store the source increment in subsurface 301 ! runoff/baseflow (runoff2). 302 ! note - runoff2 is then a negative value (as a flag) over sea-ice or 303 ! glacial-ice, in order to achieve water balance. 304 305 if (ice == 1) then 306 307 if (sneqv < 0.01) then 308 sneqv = 0.01 309 snowh = 0.10 310 ! snowh = sneqv / sndens 311 endif 312 313 elseif (ice == -1) then 314 315 if (sneqv < 0.10) then 316 ! sndens = sneqv / snowh 317 ! runoff2 = -(0.10 - sneqv) / dt 318 sneqv = 0.10 319 snowh = 1.00 320 ! snowh = sneqv / sndens 321 endif 322 323 endif ! end if_ice_block 324 325 ! --- ... for sea-ice and glacial-ice cases, set smc and sh2o values = 1 326 ! as a flag for non-soil medium 327 328 if (ice /= 0) then 329 do kz = 1, nsoil 330 smc (kz) = 1.0 331 sh2o(kz) = 1.0 332 enddo 333 endif 334 335 ! --- ... if input snowpack is nonzero, then compute snow density "sndens" 336 ! and snow thermal conductivity "sncond" (note that csnow is a 337 ! function subroutine) 338 339 if (sneqv .eq. 0.0) then 340 sndens = 0.0 341 snowh = 0.0 342 sncond = 1.0 343 else 344 sndens = sneqv / snowh 345 sndens = max( 0.0, min( 1.0, sndens )) ! added by moorthi 346 347 call csnow 348 ! --- inputs: ! 349 ! ( sndens, ! 350 ! --- outputs: ! 351 ! sncond ) ! 352 353 endif 354 355 ! --- ... determine if it's precipitating and what kind of precip it is. 356 ! if it's prcping and the air temp is colder than 0 c, it's snowing! 357 ! if it's prcping and the air temp is warmer than 0 c, but the grnd 358 ! temp is colder than 0 c, freezing rain is presumed to be falling. 359 360 if (prcp > 0.0) then 361 if (ffrozp > 0.5) then 362 snowng = .true. 363 else 364 if (t1 <= tfreez) frzgra = .true. 365 endif 366 endif 367 368 ! --- ... if either prcp flag is set, determine new snowfall (converting 369 ! prcp rate from kg m-2 s-1 to a liquid equiv snow depth in meters) 370 ! and add it to the existing snowpack. 371 ! note - that since all precip is added to snowpack, no precip infiltrates 372 ! into the soil so that prcp1 is set to zero. 373 374 if (snowng .or. frzgra) then 375 376 sn_new = prcp * dt * 0.001 377 sneqv = sneqv + sn_new 378 prcp1 = 0.0 379 380 ! --- ... update snow density based on new snowfall, using old and new 381 ! snow. update snow thermal conductivity 382 383 call snow_new 384 ! --- inputs: ! 385 ! ( sfctmp, sn_new, ! 386 ! --- input/outputs: ! 387 ! snowh, sndens ) ! 388 389 call csnow 390 ! --- inputs: ! 391 ! ( sndens, ! 392 ! --- outputs: ! 393 ! sncond ) ! 394 395 else 396 397 ! --- ... precip is liquid (rain), hence save in the precip variable 398 ! that later can wholely or partially infiltrate the soil (along 399 ! with any canopy "drip" added to this later) 400 401 prcp1 = prcp 402 403 endif ! end if_snowng_block 404 405 ! --- ... determine snowcover fraction and albedo fraction over land. 406 407 if (ice /= 0) then 408 409 ! --- ... snow cover, albedo over sea-ice, glacial-ice 410 411 sncovr = 1.0 412 albedo = 0.65 413 414 else 415 416 ! --- ... non-glacial land 417 ! if snow depth=0, set snowcover fraction=0, albedo=snow free albedo. 418 419 if (sneqv == 0.0) then 420 421 sncovr = 0.0 422 albedo = alb 423 424 else 425 426 ! --- ... determine snow fraction cover. 427 ! determine surface albedo modification due to snowdepth state. 428 429 call snfrac 430 ! --- inputs: ! 431 ! ( sneqv, snup, salp, snowh, ! 432 ! --- outputs: ! 433 ! sncovr ) ! 434 435 call alcalc 436 ! --- inputs: ! 437 ! ( alb, snoalb, shdfac, shdmin, sncovr, tsnow, ! 438 ! --- outputs: ! 439 ! albedo ) ! 440 441 endif ! end if_sneqv_block 442 443 endif ! end if_ice_block 444 445 ! --- ... thermal conductivity for sea-ice case, glacial-ice case 446 447 if (ice /= 0) then 448 449 df1 = 2.2 450 451 else 452 453 ! --- ... next calculate the subsurface heat flux, which first requires 454 ! calculation of the thermal diffusivity. treatment of the 455 ! latter follows that on pages 148-149 from "heat transfer in 456 ! cold climates", by v. j. lunardini (published in 1981 457 ! by van nostrand reinhold co.) i.e. treatment of two contiguous 458 ! "plane parallel" mediums (namely here the first soil layer 459 ! and the snowpack layer, if any). this diffusivity treatment 460 ! behaves well for both zero and nonzero snowpack, including the 461 ! limit of very thin snowpack. this treatment also eliminates 462 ! the need to impose an arbitrary upper bound on subsurface 463 ! heat flux when the snowpack becomes extremely thin. 464 465 ! --- ... first calculate thermal diffusivity of top soil layer, using 466 ! both the frozen and liquid soil moisture, following the 467 ! soil thermal diffusivity function of peters-lidard et al. 468 ! (1998,jas, vol 55, 1209-1224), which requires the specifying 469 ! the quartz content of the given soil class (see routine redprm) 470 471 call tdfcnd & 472 ! --- inputs: 473 & ( smc(1), quartz, smcmax, sh2o(1), & 474 ! --- outputs: 475 & df1 & 476 & ) 477 478 ! --- ... next add subsurface heat flux reduction effect from the 479 ! overlying green canopy, adapted from section 2.1.2 of 480 ! peters-lidard et al. (1997, jgr, vol 102(d4)) 481 482 df1 = df1 * exp( sbeta*shdfac ) 483 484 endif ! end if_ice_block 485 486 ! --- ... finally "plane parallel" snowpack effect following 487 ! v.j. linardini reference cited above. note that dtot is 488 ! combined depth of snowdepth and thickness of first soil layer 489 490 dsoil = -0.5 * zsoil(1) 491 492 if (sneqv == 0.0) then 493 494 ssoil = df1 * (t1 - stc(1)) / dsoil 495 496 else 497 498 dtot = snowh + dsoil 499 frcsno = snowh / dtot 500 frcsoi = dsoil / dtot 501 502 ! --- ... 1. harmonic mean (series flow) 503 504 ! df1 = (sncond*df1) / (frcsoi*sncond + frcsno*df1) 505 ! df1h = (sncond*df1) / (frcsoi*sncond + frcsno*df1) 506 507 ! --- ... 2. arithmetic mean (parallel flow) 508 509 ! df1 = frcsno*sncond + frcsoi*df1 510 df1a = frcsno*sncond + frcsoi*df1 511 512 ! --- ... 3. geometric mean (intermediate between harmonic and arithmetic mean) 513 514 ! df1 = (sncond**frcsno) * (df1**frcsoi) 515 ! df1 = df1h*sncovr + df1a*(1.0-sncovr) 516 ! df1 = df1h*sncovr + df1 *(1.0-sncovr) 517 df1 = df1a*sncovr + df1 *(1.0-sncovr) 518 519 ! --- ... calculate subsurface heat flux, ssoil, from final thermal 520 ! diffusivity of surface mediums, df1 above, and skin 521 ! temperature and top mid-layer soil temperature 522 523 ssoil = df1 * (t1 - stc(1)) / dtot 524 525 endif ! end if_sneqv_block 526 527 ! --- ... determine surface roughness over snowpack using snow condition 528 ! from the previous timestep. 529 530 if (couple == 0) then ! uncoupled mode 531 if (sncovr > 0.0) then 532 533 call snowz0 534 ! --- inputs: ! 535 ! ( sncovr, ! 536 ! --- input/outputs: ! 537 ! z0 ) ! 538 539 endif 540 endif 541 542 ! --- ... calc virtual temps and virtual potential temps needed by 543 ! subroutines sfcdif and penman. 544 545 t2v = sfctmp * (1.0 + 0.61*q2) 546 547 ! --- ... next call routine sfcdif to calculate the sfc exchange coef (ch) 548 ! for heat and moisture. 549 ! note - comment out call sfcdif, if sfcdif already called in calling 550 ! program (such as in coupled atmospheric model). 551 ! - do not call sfcdif until after above call to redprm, in case 552 ! alternative values of roughness length (z0) and zilintinkevich 553 ! coef (czil) are set there via namelist i/o. 554 ! - routine sfcdif returns a ch that represents the wind spd times 555 ! the "original" nondimensional "ch" typical in literature. hence 556 ! the ch returned from sfcdif has units of m/s. the important 557 ! companion coefficient of ch, carried here as "rch", is the ch 558 ! from sfcdif times air density and parameter "cp". "rch" is 559 ! computed in "call penman". rch rather than ch is the coeff 560 ! usually invoked later in eqns. 561 ! - sfcdif also returns the surface exchange coefficient for momentum, 562 ! cm, also known as the surface drage coefficient, but cm is not 563 ! used here. 564 565 ! --- ... key required radiation term is the total downward radiation 566 ! (fdown) = net solar (swnet) + downward longwave (lwdn), 567 ! for use in penman ep calculation (penman) and other surface 568 ! energy budget calcuations. also need downward solar (swdn) 569 ! for canopy resistance routine (canres). 570 ! note - fdown, swdn are derived differently in the uncoupled and 571 ! coupled modes. 572 573 if (couple == 0) then !......uncoupled mode 574 575 ! --- ... uncoupled mode: 576 ! compute surface exchange coefficients 577 578 t1v = t1 * (1.0 + 0.61 * q2) 579 th2v = th2 * (1.0 + 0.61 * q2) 580 581 call sfcdif 582 ! --- inputs: ! 583 ! ( zlvl, z0, t1v, th2v, sfcspd, czil, ! 584 ! --- input/outputs: ! 585 ! cm, ch ) ! 586 587 ! swnet = net solar radiation into the ground (w/m2; dn-up) from input 588 ! fdown = net solar + downward lw flux at sfc (w/m2) 589 590 fdown = swnet + lwdn 591 592 else !......coupled mode 593 594 ! --- ... coupled mode (couple .ne. 0): 595 ! surface exchange coefficients computed externally and passed in, 596 ! hence subroutine sfcdif not called. 597 598 ! swnet = net solar radiation into the ground (w/m2; dn-up) from input 599 ! fdown = net solar + downward lw flux at sfc (w/m2) 600 601 fdown = swnet + lwdn 602 603 endif ! end if_couple_block 604 605 ! --- ... call penman subroutine to calculate potential evaporation (etp), 606 ! and other partial products and sums save in common/rite for later 607 ! calculations. 608 609 call penman 610 ! --- inputs: ! 611 ! ( sfctmp, sfcprs, sfcems, ch, t2v, th2, prcp, fdown, ! 612 ! ssoil, q2, q2sat, dqsdt2, snowng, frzgra, ! 613 ! --- outputs: ! 614 ! t24, etp, rch, epsca, rr, flx2 ) ! 615 616 ! --- ... call canres to calculate the canopy resistance and convert it 617 ! into pc if nonzero greenness fraction 618 619 if (shdfac > 0.) then 620 621 ! --- ... frozen ground extension: total soil water "smc" was replaced 622 ! by unfrozen soil water "sh2o" in call to canres below 623 624 call canres 625 ! --- inputs: ! 626 ! ( nsoil, nroot, swdn, ch, q2, q2sat, dqsdt2, sfctmp, ! 627 ! sfcprs, sfcems, sh2o, smcwlt, smcref, zsoil, rsmin, ! 628 ! rsmax, topt, rgl, hs, xlai, ! 629 ! --- outputs: ! 630 ! rc, pc, rcs, rct, rcq, rcsoil ) ! 631 632 endif 633 634 ! --- ... now decide major pathway branch to take depending on whether 635 ! snowpack exists or not: 636 637 esnow = 0.0 638 639 if (sneqv .eq. 0.0) then 640 641 call nopac 642 ! --- inputs: ! 643 ! ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, ! 644 ! smcdry, cmcmax, dt, shdfac, sbeta, sfctmp, sfcems, ! 645 ! t24, th2, fdown, epsca, bexp, pc, rch, rr, cfactr, ! 646 ! slope, kdt, frzx, psisat, zsoil, dksat, dwsat, ! 647 ! zbot, ice, rtdis, quartz, fxexp, csoil, ! 648 ! --- input/outputs: ! 649 ! cmc, t1, stc, sh2o, tbot, ! 650 ! --- outputs: ! 651 ! eta, smc, ssoil, runoff1, runoff2, runoff3, edir, ! 652 ! ec, et, ett, beta, drip, dew, flx1, flx3 ) ! 653 654 else 655 656 call snopac 657 ! --- inputs: ! 658 ! ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, smcdry, ! 659 ! cmcmax, dt, df1, sfcems, sfctmp, t24, th2, fdown, epsca, ! 660 ! bexp, pc, rch, rr, cfactr, slope, kdt, frzx, psisat, ! 661 ! zsoil, dwsat, dksat, zbot, shdfac, ice, rtdis, quartz, ! 662 ! fxexp, csoil, flx2, snowng, ! 663 ! --- input/outputs: ! 664 ! prcp1, cmc, t1, stc, sncovr, sneqv, sndens, snowh, ! 665 ! sh2o, tbot, beta, ! 666 ! --- outputs: ! 667 ! smc, ssoil, runoff1, runoff2, runoff3, edir, ec, et, ! 668 ! ett, snomlt, drip, dew, flx1, flx3, esnow ) ! 669 670 endif 671 672 ! --- ... prepare sensible heat (h) for return to parent model 673 674 sheat = -(ch*cp1*sfcprs) / (rd1*t2v) * (th2 - t1) 675 676 ! --- ... convert units and/or sign of total evap (eta), potential evap (etp), 677 ! subsurface heat flux (s), and runoffs for what parent model expects 678 ! convert eta from kg m-2 s-1 to w m-2 679 ! eta = eta * lsubc 680 ! etp = etp * lsubc 681 682 edir = edir * lsubc 683 ec = ec * lsubc 684 685 do k = 1, 4 686 et(k) = et(k) * lsubc 687 enddo 688 689 ett = ett * lsubc 690 esnow = esnow * lsubs 691 etp = etp * ((1.0 - sncovr)*lsubc + sncovr*lsubs) 692 693 if (etp > 0.) then 694 eta = edir + ec + ett + esnow 695 else 696 eta = etp 697 endif 698 699 beta = eta / etp 700 701 ! --- ... convert the sign of soil heat flux so that: 702 ! ssoil>0: warm the surface (night time) 703 ! ssoil<0: cool the surface (day time) 704 705 ssoil = -1.0 * ssoil 706 707 if (ice == 0) then 708 709 ! --- ... for the case of land (but not glacial-ice): 710 ! convert runoff3 (internal layer runoff from supersat) from m 711 ! to m s-1 and add to subsurface runoff/baseflow (runoff2). 712 ! runoff2 is already a rate at this point. 713 714 runoff3 = runoff3 / dt 715 runoff2 = runoff2 + runoff3 716 717 else 718 719 ! --- ... for the case of sea-ice (ice=1) or glacial-ice (ice=-1), add any 720 ! snowmelt directly to surface runoff (runoff1) since there is no 721 ! soil medium, and thus no call to subroutine smflx (for soil 722 ! moisture tendency). 723 724 runoff1 = snomlt / dt 725 endif 726 727 ! --- ... total column soil moisture in meters (soilm) and root-zone 728 ! soil moisture availability (fraction) relative to porosity/saturation 729 730 soilm = -1.0 * smc(1) * zsoil(1) 731 do k = 2, nsoil 732 soilm = soilm + smc(k)*(zsoil(k-1) - zsoil(k)) 733 enddo 734 735 soilwm = -1.0 * (smcmax - smcwlt) * zsoil(1) 736 soilww = -1.0 * (smc(1) - smcwlt) * zsoil(1) 737 do k = 2, nroot 738 soilwm = soilwm + (smcmax - smcwlt) * (zsoil(k-1) - zsoil(k)) 739 soilww = soilww + (smc(k) - smcwlt) * (zsoil(k-1) - zsoil(k)) 740 enddo 741 742 soilw = soilww / soilwm 743 ! 744 return 745 746 747 ! ================= 748 contains 749 ! ================= 750 751 !*************************************! 752 ! section-1 1st level subprograms ! 753 !*************************************! 754 755 !----------------------------------- 756 subroutine alcalc 757 !................................... 758 ! --- inputs: 759 ! & ( alb, snoalb, shdfac, shdmin, sncovr, tsnow, & 760 ! --- outputs: 761 ! & albedo & 762 ! & ) 763 764 ! ===================================================================== ! 765 ! description: ! 766 ! ! 767 ! subroutine alcalc calculates albedo including snow effect (0 -> 1) ! 768 ! ! 769 ! subprogram called: none ! 770 ! ! 771 ! ! 772 ! ==================== defination of variables ==================== ! 773 ! ! 774 ! inputs from calling program: size ! 775 ! alb - real, snowfree albedo 1 ! 776 ! snoalb - real, maximum (deep) snow albedo 1 ! 777 ! shdfac - real, areal fractional coverage of green veg. 1 ! 778 ! shdmin - real, minimum areal coverage of green veg. 1 ! 779 ! sncovr - real, fractional snow cover 1 ! 780 ! tsnow - real, snow surface temperature (k) 1 ! 781 ! ! 782 ! outputs to calling program: ! 783 ! albedo - real, surface albedo including snow effect 1 ! 784 ! ! 785 ! ==================== end of description ===================== ! 786 ! 787 ! --- inputs: 788 ! real (kind=kind_phys), intent(in) :: alb, snoalb, shdfac, & 789 ! & shdmin, sncovr, tsnow 790 791 ! --- outputs: 792 ! real (kind=kind_phys), intent(out) :: albedo 793 794 ! --- locals: (none) 795 796 ! 797 !===> ... begin here 798 ! 799 ! --- ... snoalb is argument representing maximum albedo over deep snow, 800 ! as passed into sflx, and adapted from the satellite-based 801 ! maximum snow albedo fields provided by d. robinson and g. kukla 802 ! (1985, jcam, vol 24, 402-411) 803 804 ! albedo = alb + (1.0-(shdfac-shdmin))*sncovr*(snoalb-alb) 805 albedo = alb + sncovr*(snoalb - alb) 806 807 if (albedo > snoalb) albedo = snoalb 808 809 ! --- ... base formulation (dickinson et al., 1986, cogley et al., 1990) 810 811 ! if (tsnow <= 263.16) then 812 ! albedo = snoalb 813 ! else 814 ! if (tsnow < 273.16) then 815 ! tm = 0.1 * (tsnow - 263.16) 816 ! albedo = 0.5 * ((0.9 - 0.2*(tm**3)) + (0.8 - 0.16*(tm**3))) 817 ! else 818 ! albedo = 0.67 819 ! endif 820 ! endif 821 822 ! --- ... isba formulation (verseghy, 1991; baker et al., 1990) 823 824 ! if (tsnow < 273.16) then 825 ! albedo = snoalb - 0.008*dt/86400 826 ! else 827 ! albedo = (snoalb - 0.5) * exp( -0.24*dt/86400 ) + 0.5 828 ! endif 829 830 ! 831 return 832 !................................... 833 end subroutine alcalc 834 !----------------------------------- 835 836 837 !----------------------------------- 838 subroutine canres 839 !................................... 840 ! --- inputs: 841 ! & ( nsoil, nroot, swdn, ch, q2, q2sat, dqsdt2, sfctmp, & 842 ! & sfcprs, sfcems, sh2o, smcwlt, smcref, zsoil, rsmin, & 843 ! & rsmax, topt, rgl, hs, xlai, & 844 ! --- outputs: 845 ! & rc, pc, rcs, rct, rcq, rcsoil & 846 ! & ) 847 848 ! ===================================================================== ! 849 ! description: ! 850 ! ! 851 ! subroutine canres calculates canopy resistance which depends on ! 852 ! incoming solar radiation, air temperature, atmospheric water vapor ! 853 ! pressure deficit at the lowest model level, and soil moisture ! 854 ! (preferably unfrozen soil moisture rather than total) ! 855 ! ! 856 ! source: jarvis (1976), noilhan and planton (1989, mwr), jacquemin ! 857 ! and noilhan (1990, blm) ! 858 ! see also: chen et al (1996, jgr, vol 101(d3), 7251-7268), eqns ! 859 ! 12-14 and table 2 of sec. 3.1.2 ! 860 ! ! 861 ! subprogram called: none ! 862 ! ! 863 ! ! 864 ! ==================== defination of variables ==================== ! 865 ! ! 866 ! inputs from calling program: size ! 867 ! nsoil - integer, no. of soil layers 1 ! 868 ! nroot - integer, no. of soil layers in root zone (<nsoil) 1 ! 869 ! swdn - real, incoming solar radiation 1 ! 870 ! ch - real, sfc exchange coeff for heat and moisture 1 ! 871 ! q2 - real, air humidity at 1st level above ground 1 ! 872 ! q2sat - real, sat. air humidity at 1st level abv ground 1 ! 873 ! dqsdt2 - real, slope of sat. humidity function wrt temp 1 ! 874 ! sfctmp - real, sfc temperature at 1st level above ground 1 ! 875 ! sfcprs - real, sfc pressure 1 ! 876 ! sfcems - real, sfc emissivity for lw radiation 1 ! 877 ! sh2o - real, volumetric soil moisture nsoil ! 878 ! smcwlt - real, wilting point 1 ! 879 ! smcref - real, reference soil moisture 1 ! 880 ! zsoil - real, soil depth (negative sign, as below grd) nsoil ! 881 ! rsmin - real, mimimum stomatal resistance 1 ! 882 ! rsmax - real, maximum stomatal resistance 1 ! 883 ! topt - real, optimum transpiration air temperature 1 ! 884 ! rgl - real, canopy resistance func (in solar rad term) 1 ! 885 ! hs - real, canopy resistance func (vapor deficit term) 1 ! 886 ! xlai - real, leaf area index 1 ! 887 ! ! 888 ! outputs to calling program: ! 889 ! rc - real, canopy resistance 1 ! 890 ! pc - real, plant coefficient 1 ! 891 ! rcs - real, incoming solar rc factor 1 ! 892 ! rct - real, air temp rc factor 1 ! 893 ! rcq - real, atoms vapor press deficit rc factor 1 ! 894 ! rcsoil - real, soil moisture rc factor 1 ! 895 ! ! 896 ! ==================== end of description ===================== ! 897 ! 898 ! --- inputs: 899 ! integer, intent(in) :: nsoil, nroot 900 901 ! real (kind=kind_phys), intent(in) :: swdn, ch, q2, q2sat, & 902 ! & dqsdt2, sfctmp, sfcprs, sfcems, smcwlt, smcref, rsmin, & 903 ! & rsmax, topt, rgl, hs, xlai, sh2o(nsoil), zsoil(nsoil) 904 905 ! --- outputs: 906 ! real (kind=kind_phys), intent(out) :: rc, pc, rcs, rct, rcq, & 907 ! & rcsoil 908 909 ! --- locals: 910 real (kind=kind_phys) :: delta, ff, gx, rr, part(nsold) 911 912 integer :: k 913 914 ! 915 !===> ... begin here 916 ! 917 ! --- ... initialize canopy resistance multiplier terms. 918 919 rcs = 0.0 920 rct = 0.0 921 rcq = 0.0 922 rcsoil = 0.0 923 rc = 0.0 924 925 ! --- ... contribution due to incoming solar radiation 926 927 ff = 0.55 * 2.0 * swdn / (rgl*xlai) 928 rcs = (ff + rsmin/rsmax) / (1.0 + ff) 929 rcs = max( rcs, 0.0001 ) 930 931 ! --- ... contribution due to air temperature at first model level above ground 932 ! rct expression from noilhan and planton (1989, mwr). 933 934 rct = 1.0 - 0.0016 * (topt - sfctmp)**2.0 935 rct = max( rct, 0.0001 ) 936 937 ! --- ... contribution due to vapor pressure deficit at first model level. 938 ! rcq expression from ssib 939 940 rcq = 1.0 / (1.0 + hs*(q2sat-q2)) 941 rcq = max( rcq, 0.01 ) 942 943 ! --- ... contribution due to soil moisture availability. 944 ! determine contribution from each soil layer, then add them up. 945 946 gx = (sh2o(1) - smcwlt) / (smcref - smcwlt) 947 gx = max( 0.0, min( 1.0, gx ) ) 948 949 ! --- ... use soil depth as weighting factor 950 part(1) = (zsoil(1)/zsoil(nroot)) * gx 951 952 ! --- ... use root distribution as weighting factor 953 ! part(1) = rtdis(1) * gx 954 955 do k = 2, nroot 956 957 gx = (sh2o(k) - smcwlt) / (smcref - smcwlt) 958 gx = max( 0.0, min( 1.0, gx ) ) 959 960 ! --- ... use soil depth as weighting factor 961 part(k) = ((zsoil(k) - zsoil(k-1)) / zsoil(nroot)) * gx 962 963 ! --- ... use root distribution as weighting factor 964 ! part(k) = rtdis(k) * gx 965 966 enddo 967 968 do k = 1, nroot 969 rcsoil = rcsoil + part(k) 970 enddo 971 rcsoil = max( rcsoil, 0.0001 ) 972 973 ! --- ... determine canopy resistance due to all factors. convert canopy 974 ! resistance (rc) to plant coefficient (pc) to be used with 975 ! potential evap in determining actual evap. pc is determined by: 976 ! pc * linerized penman potential evap = penman-monteith actual 977 ! evaporation (containing rc term). 978 979 rc = rsmin / (xlai*rcs*rct*rcq*rcsoil) 980 rr = (4.0*sfcems*sigma1*rd1/cp1) * (sfctmp**4.0)/(sfcprs*ch) + 1.0 981 delta = (lsubc/cp1) * dqsdt2 982 983 pc = (rr + delta) / (rr*(1.0 + rc*ch) + delta) 984 ! 985 return 986 !................................... 987 end subroutine canres 988 !----------------------------------- 989 990 991 !----------------------------------- 992 subroutine csnow 993 !................................... 994 ! --- inputs: 995 ! & ( sndens, & 996 ! --- outputs: 997 ! & sncond & 998 ! & ) 999 1000 ! ===================================================================== ! 1001 ! description: ! 1002 ! ! 1003 ! subroutine csnow calculates snow termal conductivity ! 1004 ! ! 1005 ! subprogram called: none ! 1006 ! ! 1007 ! ==================== defination of variables ==================== ! 1008 ! ! 1009 ! inputs from the calling program: size ! 1010 ! sndens - real, snow density 1 ! 1011 ! ! 1012 ! outputs to the calling program: ! 1013 ! sncond - real, snow termal conductivity 1 ! 1014 ! ! 1015 ! ==================== end of description ===================== ! 1016 ! 1017 ! --- constant parameters: 1018 real (kind=kind_phys), parameter :: unit = 0.11631 1019 1020 ! --- inputs: 1021 ! real (kind=kind_phys), intent(in) :: sndens 1022 1023 ! --- outputs: 1024 ! real (kind=kind_phys), intent(out) :: sncond 1025 1026 ! --- locals: 1027 real (kind=kind_phys) :: c 1028 1029 ! 1030 !===> ... begin here 1031 ! 1032 ! --- ... sncond in units of cal/(cm*hr*c), returned in w/(m*c) 1033 ! basic version is dyachkova equation (1960), for range 0.1-0.4 1034 1035 c = 0.328 * 10**(2.25*sndens) 1036 sncond = unit * c 1037 1038 ! --- ... de vaux equation (1933), in range 0.1-0.6 1039 1040 ! sncond = 0.0293 * (1.0 + 100.0*sndens**2) 1041 1042 ! --- ... e. andersen from flerchinger 1043 1044 ! sncond = 0.021 + 2.51 * sndens**2 1045 ! 1046 return 1047 !................................... 1048 end subroutine csnow 1049 !----------------------------------- 1050 1051 1052 !----------------------------------- 1053 subroutine nopac 1054 !................................... 1055 ! --- inputs: 1056 ! & ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, & 1057 ! & smcdry, cmcmax, dt, shdfac, sbeta, sfctmp, sfcems, & 1058 ! & t24, th2, fdown, epsca, bexp, pc, rch, rr, cfactr, & 1059 ! & slope, kdt, frzx, psisat, zsoil, dksat, dwsat, & 1060 ! & zbot, ice, rtdis, quartz, fxexp, csoil, & 1061 ! --- input/outputs: 1062 ! & cmc, t1, stc, sh2o, tbot, & 1063 ! --- outputs: 1064 ! & eta, smc, ssoil, runoff1, runoff2, runoff3, edir, & 1065 ! & ec, et, ett, beta, drip, dew, flx1, flx3 & 1066 ! & ) 1067 1068 ! ===================================================================== ! 1069 ! description: ! 1070 ! ! 1071 ! subroutine nopac calculates soil moisture and heat flux values and ! 1072 ! update soil moisture content and soil heat content values for the ! 1073 ! case when no snow pack is present. ! 1074 ! ! 1075 ! ! 1076 ! subprograms called: evapo, smflx, tdfcnd, shflx ! 1077 ! ! 1078 ! ==================== defination of variables ==================== ! 1079 ! ! 1080 ! inputs from calling program: size ! 1081 ! nsoil - integer, number of soil layers 1 ! 1082 ! nroot - integer, number of root layers 1 ! 1083 ! etp - real, potential evaporation 1 ! 1084 ! prcp - real, precip rate 1 ! 1085 ! smcmax - real, porosity (sat val of soil mois) 1 ! 1086 ! smcwlt - real, wilting point 1 ! 1087 ! smcref - real, soil mois threshold 1 ! 1088 ! smcdry - real, dry soil mois threshold 1 ! 1089 ! cmcmax - real, maximum canopy water parameters 1 ! 1090 ! dt - real, time step 1 ! 1091 ! shdfac - real, aeral coverage of green veg 1 ! 1092 ! sbeta - real, param to cal veg effect on soil heat flux 1 ! 1093 ! sfctmp - real, air temp at height zlvl abv ground 1 ! 1094 ! sfcems - real, sfc lw emissivity 1 ! 1095 ! t24 - real, sfctmp**4 1 ! 1096 ! th2 - real, air potential temp at zlvl abv grnd 1 ! 1097 ! fdown - real, net solar + downward lw flux at sfc 1 ! 1098 ! epsca - real, 1 ! 1099 ! bexp - real, soil type "b" parameter 1 ! 1100 ! pc - real, plant coeff 1 ! 1101 ! rch - real, companion coefficient of ch 1 ! 1102 ! rr - real, 1 ! 1103 ! cfactr - real, canopy water parameters 1 ! 1104 ! slope - real, linear reservoir coefficient 1 ! 1105 ! kdt - real, 1 ! 1106 ! frzx - real, frozen ground parameter 1 ! 1107 ! psisat - real, saturated soil potential 1 ! 1108 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 1109 ! dksat - real, saturated soil hydraulic conductivity 1 ! 1110 ! dwsat - real, saturated soil diffusivity 1 ! 1111 ! zbot - real, specify depth of lower bd soil 1 ! 1112 ! ice - integer, sea-ice flag (=1: sea-ice, =0: land) 1 ! 1113 ! rtdis - real, root distribution nsoil ! 1114 ! quartz - real, soil quartz content 1 ! 1115 ! fxexp - real, bare soil evaporation exponent 1 ! 1116 ! csoil - real, soil heat capacity 1 ! 1117 ! ! 1118 ! input/outputs from and to the calling program: ! 1119 ! cmc - real, canopy moisture content 1 ! 1120 ! t1 - real, ground/canopy/snowpack eff skin temp 1 ! 1121 ! stc - real, soil temp nsoil ! 1122 ! sh2o - real, unfrozen soil moisture nsoil ! 1123 ! tbot - real, bottom soil temp 1 ! 1124 ! ! 1125 ! outputs to the calling program: ! 1126 ! eta - real, downward latent heat flux 1 ! 1127 ! smc - real, total soil moisture nsoil ! 1128 ! ssoil - real, upward soil heat flux 1 ! 1129 ! runoff1 - real, surface runoff not infiltrating sfc 1 ! 1130 ! runoff2 - real, sub surface runoff (baseflow) 1 ! 1131 ! runoff3 - real, excess of porosity 1 ! 1132 ! edir - real, direct soil evaporation 1 ! 1133 ! ec - real, canopy water evaporation 1 ! 1134 ! et - real, plant transpiration nsoil ! 1135 ! ett - real, total plant transpiration 1 ! 1136 ! beta - real, ratio of actual/potential evap 1 ! 1137 ! drip - real, through-fall of precip and/or dew 1 ! 1138 ! dew - real, dewfall (or frostfall) 1 ! 1139 ! flx1 - real, precip-snow sfc flux 1 ! 1140 ! flx3 - real, phase-change heat flux from snowmelt 1 ! 1141 ! ! 1142 ! ==================== end of description ===================== ! 1143 ! 1144 ! --- inputs: 1145 ! integer, intent(in) :: nsoil, nroot, ice 1146 1147 ! real (kind=kind_phys), intent(in) :: etp, prcp, smcmax, & 1148 ! & smcwlt, smcref, smcdry, cmcmax, dt, shdfac, sbeta, & 1149 ! & sfctmp, sfcems, t24, th2, fdown, epsca, bexp, pc, & 1150 ! & rch, rr, cfactr, slope, kdt, frzx, psisat, & 1151 ! & zsoil(nsoil), dksat, dwsat, zbot, rtdis(nsoil), & 1152 ! & quartz, fxexp, csoil 1153 1154 ! --- input/outputs: 1155 ! real (kind=kind_phys), intent(inout) :: cmc, t1, stc(nsoil), & 1156 ! & sh2o(nsoil), tbot 1157 1158 ! --- outputs: 1159 ! real (kind=kind_phys), intent(out) :: eta, smc(nsoil), ssoil, & 1160 ! & runoff1, runoff2, runoff3, edir, ec, et(nsoil), ett, & 1161 ! & beta, drip, dew, flx1, flx3 1162 1163 ! --- locals: 1164 real (kind=kind_phys) :: df1, eta1, etp1, prcp1, yy, yynum, & 1165 & zz1, ec1, edir1, et1(nsoil), ett1 1166 1167 integer :: k 1168 1169 ! 1170 !===> ... begin here 1171 ! 1172 ! --- ... convert etp from kg m-2 s-1 to ms-1 and initialize dew. 1173 1174 prcp1= prcp * 0.001 1175 etp1 = etp * 0.001 1176 dew = 0.0 1177 edir = 0.0 1178 edir1= 0.0 1179 ec = 0.0 1180 ec1 = 0.0 1181 1182 do k = 1, nsoil 1183 et (k) = 0.0 1184 et1(k) = 0.0 1185 enddo 1186 1187 ett = 0.0 1188 ett1 = 0.0 1189 1190 if (etp > 0.0) then 1191 1192 ! --- ... convert prcp from 'kg m-2 s-1' to 'm s-1'. 1193 1194 call evapo & 1195 ! --- inputs: 1196 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, & 1197 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, & 1198 & shdfac, cfactr, rtdis, fxexp, & 1199 ! --- outputs: 1200 & eta1, edir1, ec1, et1, ett1 & 1201 & ) 1202 1203 call smflx & 1204 ! --- inputs: 1205 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, & 1206 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, & 1207 & edir1, ec1, et1, & 1208 ! --- input/outputs: 1209 & cmc, sh2o, & 1210 ! --- outputs: 1211 & smc, runoff1, runoff2, runoff3, drip & 1212 & ) 1213 1214 else 1215 1216 ! --- ... if etp < 0, assume dew forms (transform etp1 into dew and 1217 ! reinitialize etp1 to zero). 1218 1219 eta1 = 0.0 1220 dew = -etp1 1221 1222 ! --- ... convert prcp from 'kg m-2 s-1' to 'm s-1' and add dew amount. 1223 1224 prcp1 = prcp1 + dew 1225 1226 call smflx & 1227 ! --- inputs: 1228 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, & 1229 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, & 1230 & edir1, ec1, et1, & 1231 ! --- input/outputs: 1232 & cmc, sh2o, & 1233 ! --- outputs: 1234 & smc, runoff1, runoff2, runoff3, drip & 1235 & ) 1236 1237 endif ! end if_etp_block 1238 1239 ! --- ... convert modeled evapotranspiration fm m s-1 to kg m-2 s-1 1240 1241 eta = eta1 * 1000.0 1242 edir = edir1 * 1000.0 1243 ec = ec1 * 1000.0 1244 1245 do k = 1, nsoil 1246 et(k) = et1(k) * 1000.0 1247 enddo 1248 1249 ett = ett1 * 1000.0 1250 1251 ! --- ... based on etp and e values, determine beta 1252 1253 if ( etp <= 0.0 ) then 1254 beta = 0.0 1255 if ( etp < 0.0 ) then 1256 beta = 1.0 1257 endif 1258 else 1259 beta = eta / etp 1260 endif 1261 1262 ! --- ... get soil thermal diffuxivity/conductivity for top soil lyr, 1263 ! calc. adjusted top lyr soil temp and adjusted soil flux, then 1264 ! call shflx to compute/update soil heat flux and soil temps. 1265 1266 call tdfcnd & 1267 ! --- inputs: 1268 & ( smc(1), quartz, smcmax, sh2o(1), & 1269 ! --- outputs: 1270 & df1 & 1271 & ) 1272 1273 ! --- ... vegetation greenness fraction reduction in subsurface heat 1274 ! flux via reduction factor, which is convenient to apply here 1275 ! to thermal diffusivity that is later used in hrt to compute 1276 ! sub sfc heat flux (see additional comments on veg effect 1277 ! sub-sfc heat flx in routine sflx) 1278 1279 df1 = df1 * exp( sbeta*shdfac ) 1280 1281 ! --- ... compute intermediate terms passed to routine hrt (via routine 1282 ! shflx below) for use in computing subsurface heat flux in hrt 1283 1284 yynum = fdown - sfcems*sigma1*t24 1285 yy = sfctmp + (yynum/rch + th2 - sfctmp - beta*epsca)/rr 1286 zz1 = df1/(-0.5*zsoil(1)*rch*rr) + 1.0 1287 1288 call shflx & 1289 ! --- inputs: 1290 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, & 1291 & psisat, bexp, df1, ice, quartz, csoil, & 1292 ! --- input/outputs: 1293 & stc, t1, tbot, sh2o, & 1294 ! --- outputs: 1295 & ssoil & 1296 & ) 1297 1298 ! --- ... set flx1 and flx3 (snopack phase change heat fluxes) to zero since 1299 ! they are not used here in snopac. flx2 (freezing rain heat flux) 1300 ! was similarly initialized in the penman routine. 1301 1302 flx1 = 0.0 1303 flx3 = 0.0 1304 ! 1305 return 1306 !................................... 1307 end subroutine nopac 1308 !----------------------------------- 1309 1310 1311 !----------------------------------- 1312 subroutine penman 1313 !................................... 1314 ! --- inputs: 1315 ! & ( sfctmp, sfcprs, sfcems, ch, t2v, th2, prcp, fdown, & 1316 ! & ssoil, q2, q2sat, dqsdt2, snowng, frzgra, & 1317 ! --- outputs: 1318 ! & t24, etp, rch, epsca, rr, flx2 & 1319 ! & ) 1320 1321 ! ===================================================================== ! 1322 ! description: ! 1323 ! ! 1324 ! subroutine penman calculates potential evaporation for the current ! 1325 ! point. various partial sums/products are also calculated and passed ! 1326 ! back to the calling routine for later use. ! 1327 ! ! 1328 ! ! 1329 ! subprogram called: none ! 1330 ! ! 1331 ! ==================== defination of variables ==================== ! 1332 ! ! 1333 ! inputs: size ! 1334 ! sfctmp - real, sfc temperature at 1st level above ground 1 ! 1335 ! sfcprs - real, sfc pressure 1 ! 1336 ! sfcems - real, sfc emissivity for lw radiation 1 ! 1337 ! ch - real, sfc exchange coeff for heat & moisture 1 ! 1338 ! t2v - real, sfc virtual temperature 1 ! 1339 ! th2 - real, air potential temp at zlvl abv grnd 1 ! 1340 ! prcp - real, precip rate 1 ! 1341 ! fdown - real, net solar + downward lw flux at sfc 1 ! 1342 ! ssoil - real, upward soil heat flux 1 ! 1343 ! q2 - real, mixing ratio at hght zlvl abv ground 1 ! 1344 ! q2sat - real, sat mixing ratio at zlvl abv ground 1 ! 1345 ! dqsdt2 - real, slope of sat specific humidity curve 1 ! 1346 ! snowng - logical, snow flag 1 ! 1347 ! frzgra - logical, freezing rain flag 1 ! 1348 ! ! 1349 ! outputs: ! 1350 ! t24 - real, sfctmp**4 1 ! 1351 ! etp - real, potential evaporation 1 ! 1352 ! rch - real, companion coefficient of ch 1 ! 1353 ! epsca - real, 1 ! 1354 ! rr - real, 1 ! 1355 ! flx2 - real, freezing rain latent heat flux 1 ! 1356 ! ! 1357 ! ==================== end of description ===================== ! 1358 ! 1359 ! --- inputs: 1360 ! real (kind=kind_phys), intent(in) :: sfctmp, sfcprs, sfcems, & 1361 ! & ch, t2v, th2, prcp, fdown, ssoil, q2, q2sat, dqsdt2 1362 1363 ! logical, intent(in) :: snowng, frzgra 1364 1365 ! --- outputs: 1366 ! real (kind=kind_phys), intent(out) :: t24, etp, rch, epsca, & 1367 ! & rr, flx2 1368 1369 ! --- locals: 1370 real (kind=kind_phys) :: a, delta, fnet, rad, rho 1371 1372 ! 1373 !===> ... begin here 1374 ! 1375 flx2 = 0.0 1376 1377 ! --- ... prepare partial quantities for penman equation. 1378 1379 delta = elcp * dqsdt2 1380 t24 = sfctmp * sfctmp * sfctmp * sfctmp 1381 rr = t24 * 6.48e-8 / (sfcprs*ch) + 1.0 1382 rho = sfcprs / (rd1*t2v) 1383 rch = rho * cp * ch 1384 1385 ! --- ... adjust the partial sums / products with the latent heat 1386 ! effects caused by falling precipitation. 1387 1388 if (.not. snowng) then 1389 if (prcp > 0.0) rr = rr + cph2o1*prcp/rch 1390 else 1391 rr = rr + cpice*prcp/rch 1392 endif 1393 1394 fnet = fdown - sfcems*sigma1*t24 - ssoil 1395 1396 ! --- ... include the latent heat effects of frzng rain converting to ice 1397 ! on impact in the calculation of flx2 and fnet. 1398 1399 if (frzgra) then 1400 flx2 = -lsubf * prcp 1401 fnet = fnet - flx2 1402 endif 1403 1404 ! --- ... finish penman equation calculations. 1405 1406 rad = fnet/rch + th2 - sfctmp 1407 a = elcp * (q2sat - q2) 1408 epsca = (a*rr + rad*delta) / (delta + rr) 1409 etp = epsca * rch / lsubc 1410 ! 1411 return 1412 !................................... 1413 end subroutine penman 1414 !----------------------------------- 1415 1416 1417 !----------------------------------- 1418 subroutine redprm 1419 !................................... 1420 ! --- inputs: 1421 ! & ( nsoil, vegtyp, soiltyp, slopetyp, sldpth, zsoil, & 1422 ! --- outputs: 1423 ! & cfactr, cmcmax, rsmin, rsmax, topt, refkdt, kdt, & 1424 ! & sbeta, shdfac, rgl, hs, zbot, frzx, psisat, slope, & 1425 ! & snup, salp, bexp, dksat, dwsat, smcmax, smcwlt, & 1426 ! & smcref, smcdry, f1, quartz, fxexp, rtdis, nroot, & 1427 ! & z0, czil, xlai, csoil & 1428 ! & ) 1429 1430 ! ===================================================================== ! 1431 ! description: ! 1432 ! ! 1433 ! subroutine redprm internally sets(default valuess), or optionally ! 1434 ! read-in via namelist i/o, all soil and vegetation parameters ! 1435 ! required for the execusion of the noah lsm. ! 1436 ! ! 1437 ! optional non-default parameters can be read in, accommodating up to ! 1438 ! 30 soil, veg, or slope classes, if the default max number of soil, ! 1439 ! veg, and/or slope types is reset. ! 1440 ! ! 1441 ! future upgrades of routine redprm must expand to incorporate some ! 1442 ! of the empirical parameters of the frozen soil and snowpack physics ! 1443 ! (such as in routines frh2o, snowpack, and snow_new) not yet set in ! 1444 ! this redprm routine, but rather set in lower level subroutines. ! 1445 ! ! 1446 ! all soil, veg, slope, and universal parameters values are defined ! 1447 ! externally (in subroutine "set_soilveg.f") and then accessed via ! 1448 ! "use namelist_soilveg" (below) and then set here. ! 1449 ! ! 1450 ! soil types zobler (1986) cosby et al (1984) (quartz cont.(1)) ! 1451 ! 1 coarse loamy sand (0.82) ! 1452 ! 2 medium silty clay loam (0.10) ! 1453 ! 3 fine light clay (0.25) ! 1454 ! 4 coarse-medium sandy loam (0.60) ! 1455 ! 5 coarse-fine sandy clay (0.52) ! 1456 ! 6 medium-fine clay loam (0.35) ! 1457 ! 7 coarse-med-fine sandy clay loam (0.60) ! 1458 ! 8 organic loam (0.40) ! 1459 ! 9 glacial land ice loamy sand (na using 0.82)! 1460 ! 13: <old>- glacial land ice -<old> ! 1461 ! 13: glacial-ice (no longer use these parameters), now ! 1462 ! treated as ice-only surface and sub-surface ! 1463 ! (in subroutine hrtice) ! 1464 ! ! 1465 ! ssib vegetation types (dorman and sellers, 1989; jam) ! 1466 ! 1: broadleaf-evergreen trees (tropical forest) ! 1467 ! 2: broadleaf-deciduous trees ! 1468 ! 3: broadleaf and needleleaf trees (mixed forest) ! 1469 ! 4: needleleaf-evergreen trees ! 1470 ! 5: needleleaf-deciduous trees (larch) ! 1471 ! 6: broadleaf trees with groundcover (savanna) ! 1472 ! 7: groundcover only (perennial) ! 1473 ! 8: broadleaf shrubs with perennial groundcover ! 1474 ! 9: broadleaf shrubs with bare soil ! 1475 ! 10: dwarf trees and shrubs with groundcover (tundra) ! 1476 ! 11: bare soil ! 1477 ! 12: cultivations (the same parameters as for type 7) ! 1478 ! 13: <old>- glacial (the same parameters as for type 11) -<old> ! 1479 ! 13: glacial-ice (no longer use these parameters), now treated as ! 1480 ! ice-only surface and sub-surface (in subroutine hrtice) ! 1481 ! ! 1482 ! slopetyp is to estimate linear reservoir coefficient slope to the ! 1483 ! baseflow runoff out of the bottom layer. lowest class (slopetyp=0) ! 1484 ! means highest slope parameter = 1. ! 1485 ! ! 1486 ! slope class percent slope ! 1487 ! 1 0-8 ! 1488 ! 2 8-30 ! 1489 ! 3 > 30 ! 1490 ! 4 0-30 ! 1491 ! 5 0-8 & > 30 ! 1492 ! 6 8-30 & > 30 ! 1493 ! 7 0-8, 8-30, > 30 ! 1494 ! 9 glacial ice ! 1495 ! blank ocean/sea ! 1496 ! ! 1497 ! note: class 9 from zobler file should be replaced by 8 and 'blank' 9 ! 1498 ! ! 1499 ! ! 1500 ! subprogram called: none ! 1501 ! ! 1502 ! ==================== defination of variables ==================== ! 1503 ! ! 1504 ! inputs from calling program: size ! 1505 ! nsoil - integer, number of soil layers 1 ! 1506 ! vegtyp - integer, vegetation type (integer index) 1 ! 1507 ! soiltyp - integer, soil type (integer index) 1 ! 1508 ! slopetyp - integer, class of sfc slope (integer index) 1 ! 1509 ! sldpth - integer, thickness of each soil layer (m) nsoil ! 1510 ! zsoil - integer, soil depth (negative sign) (m) nsoil ! 1511 ! ! 1512 ! outputs to the calling program: ! 1513 ! cfactr - real, canopy water parameters 1 ! 1514 ! cmcmax - real, maximum canopy water parameters 1 ! 1515 ! rsmin - real, mimimum stomatal resistance 1 ! 1516 ! rsmax - real, maximum stomatal resistance 1 ! 1517 ! topt - real, optimum transpiration air temperature 1 ! 1518 ! refkdt - real, =2.e-6 the sat. dk. val for soil type 2 1 ! 1519 ! kdt - real, 1 ! 1520 ! sbeta - real, param to cal veg effect on soil heat flux 1 ! 1521 ! shdfac - real, vegetation greenness fraction 1 ! 1522 ! rgl - real, canopy resistance func (in solar rad term) 1 ! 1523 ! hs - real, canopy resistance func (vapor deficit term) 1 ! 1524 ! zbot - real, specify depth of lower bd soil temp (m) 1 ! 1525 ! frzx - real, frozen ground parameter, ice content 1 ! 1526 ! threshold above which frozen soil is impermeable ! 1527 ! psisat - real, saturated soil potential 1 ! 1528 ! slope - real, linear reservoir coefficient 1 ! 1529 ! snup - real, threshold snow depth (water equi m) 1 ! 1530 ! salp - real, snow cover shape parameter 1 ! 1531 ! from anderson's hydro-17 best fit salp = 2.6 ! 1532 ! bexp - real, the 'b' parameter 1 ! 1533 ! dksat - real, saturated soil hydraulic conductivity 1 ! 1534 ! dwsat - real, saturated soil diffusivity 1 ! 1535 ! smcmax - real, max soil moisture content (porosity) 1 ! 1536 ! smcwlt - real, wilting pt soil moisture contents 1 ! 1537 ! smcref - real, reference soil moisture (onset stress) 1 ! 1538 ! smcdry - real, air dry soil moist content limits 1 ! 1539 ! f1 - real, used to comp soil diffusivity/conductivity 1 ! 1540 ! quartz - real, soil quartz content 1 ! 1541 ! fxexp - real, bare soil evaporation exponent 1 ! 1542 ! rtdis - real, root distribution nsoil ! 1543 ! nroot - integer, number of root layers 1 ! 1544 ! z0 - real, roughness length (m) 1 ! 1545 ! czil - real, param to cal roughness length of heat 1 ! 1546 ! xlai - real, leaf area index 1 ! 1547 ! csoil - real, soil heat capacity (j m-3 k-1) 1 ! 1548 ! ! 1549 ! ==================== end of description ===================== ! 1550 ! 1551 use namelist_soilveg 1552 1553 ! --- input: 1554 ! integer, intent(in) :: nsoil, vegtyp, soiltyp, slopetyp 1555 1556 ! real (kind=kind_phys), intent(in) :: sldpth(nsoil), zsoil(nsoil) 1557 1558 ! --- outputs: 1559 ! real (kind=kind_phys), intent(out) :: cfactr, cmcmax, rsmin, & 1560 ! & rsmax, topt, refkdt, kdt, sbeta, shdfac, rgl, hs, zbot, & 1561 ! & frzx, psisat, slope, snup, salp, bexp, dksat, dwsat, & 1562 ! & smcmax, smcwlt, smcref, smcdry, f1, quartz, fxexp, z0, & 1563 ! & czil, xlai, csoil, rtdis(nsoil) 1564 1565 ! integer, intent(out) :: nroot 1566 1567 ! --- locals: 1568 real (kind=kind_phys) :: frzfact, frzk, refdk 1569 1570 integer :: i 1571 1572 ! 1573 !===> ... begin here 1574 ! 1575 if (soiltyp > defined_soil) then 1576 write(*,*) 'warning: too many soil types' 1577 stop 333 1578 endif 1579 1580 if (vegtyp > defined_veg) then 1581 write(*,*) 'warning: too many veg types' 1582 stop 333 1583 endif 1584 1585 if (slopetyp > defined_slope) then 1586 write(*,*) 'warning: too many slope types' 1587 stop 333 1588 endif 1589 1590 ! --- ... set-up universal parameters (not dependent on soiltyp, vegtyp 1591 ! or slopetyp) 1592 1593 zbot = zbot_data 1594 salp = salp_data 1595 cfactr = cfactr_data 1596 cmcmax = cmcmax_data 1597 sbeta = sbeta_data 1598 rsmax = rsmax_data 1599 topt = topt_data 1600 refdk = refdk_data 1601 frzk = frzk_data 1602 fxexp = fxexp_data 1603 refkdt = refkdt_data 1604 czil = czil_data 1605 csoil = csoil_data 1606 1607 ! --- ... set-up soil parameters 1608 1609 bexp = bb (soiltyp) 1610 dksat = satdk(soiltyp) 1611 dwsat = satdw(soiltyp) 1612 f1 = f11 (soiltyp) 1613 kdt = refkdt * dksat / refdk 1614 1615 psisat = satpsi(soiltyp) 1616 quartz = qtz (soiltyp) 1617 smcdry = drysmc(soiltyp) 1618 smcmax = maxsmc(soiltyp) 1619 smcref = refsmc(soiltyp) 1620 smcwlt = wltsmc(soiltyp) 1621 1622 frzfact = (smcmax / smcref) * (0.412 / 0.468) 1623 1624 ! --- ... to adjust frzk parameter to actual soil type: frzk * frzfact 1625 1626 frzx = frzk * frzfact 1627 1628 ! --- ... set-up vegetation parameters 1629 1630 nroot = nroot_data(vegtyp) 1631 snup = snupx(vegtyp) 1632 rsmin = rsmtbl(vegtyp) 1633 1634 rgl = rgltbl(vegtyp) 1635 hs = hstbl(vegtyp) 1636 z0 = z0_data(vegtyp) 1637 xlai= lai_data(vegtyp) 1638 1639 if (vegtyp == bare) shdfac = 0.0 1640 1641 if (nroot > nsoil) then 1642 write(*,*) 'warning: too many root layers' 1643 stop 333 1644 endif 1645 1646 ! --- ... calculate root distribution. present version assumes uniform 1647 ! distribution based on soil layer depths. 1648 1649 do i = 1, nroot 1650 rtdis(i) = -sldpth(i) / zsoil(nroot) 1651 enddo 1652 1653 ! --- ... set-up slope parameter 1654 1655 slope = slope_data(slopetyp) 1656 ! 1657 return 1658 !................................... 1659 end subroutine redprm 1660 !----------------------------------- 1661 1662 1663 !----------------------------------- 1664 subroutine sfcdif 1665 !................................... 1666 ! --- inputs: 1667 ! & ( zlvl, z0, t1v, th2v, sfcspd, czil, & 1668 ! --- input/outputs: 1669 ! & cm, ch & 1670 ! & ) 1671 1672 ! ===================================================================== ! 1673 ! description: ! 1674 ! ! 1675 ! subroutine sfcdif calculates surface layer exchange coefficients ! 1676 ! via iterative process. see chen et al (1997, blm) ! 1677 ! ! 1678 ! subprogram called: none ! 1679 ! ! 1680 ! ! 1681 ! ==================== defination of variables ==================== ! 1682 ! ! 1683 ! inputs from the calling program: size ! 1684 ! zlvl - real, height abv atmos ground forcing vars (m) 1 ! 1685 ! z0 - real, roughness length (m) 1 ! 1686 ! t1v - real, surface exchange coefficient 1 ! 1687 ! th2v - real, surface exchange coefficient 1 ! 1688 ! sfcspd - real, wind speed at height zlvl abv ground (m/s) 1 ! 1689 ! czil - real, param to cal roughness length of heat 1 ! 1690 ! ! 1691 ! input/outputs from and to the calling program: ! 1692 ! cm - real, sfc exchange coeff for momentum (m/s) 1 ! 1693 ! ch - real, sfc exchange coeff for heat & moisture (m/s)1 ! 1694 ! ! 1695 ! ==================== end of description ===================== ! 1696 ! 1697 ! --- constant parameters: 1698 integer, parameter :: itrmx = 5 1699 real (kind=kind_phys), parameter :: wwst = 1.2 1700 real (kind=kind_phys), parameter :: wwst2 = wwst*wwst 1701 real (kind=kind_phys), parameter :: vkrm = 0.40 1702 real (kind=kind_phys), parameter :: excm = 0.001 1703 real (kind=kind_phys), parameter :: beta = 1.0/270.0 1704 real (kind=kind_phys), parameter :: btg = beta*gs1 1705 real (kind=kind_phys), parameter :: elfc = vkrm*btg 1706 real (kind=kind_phys), parameter :: wold = 0.15 1707 real (kind=kind_phys), parameter :: wnew = 1.0-wold 1708 real (kind=kind_phys), parameter :: pihf = 3.14159265/2.0 ! con_pi/2.0 1709 1710 real (kind=kind_phys), parameter :: epsu2 = 1.e-4 1711 real (kind=kind_phys), parameter :: epsust = 0.07 1712 real (kind=kind_phys), parameter :: ztmin = -5.0 1713 real (kind=kind_phys), parameter :: ztmax = 1.0 1714 real (kind=kind_phys), parameter :: hpbl = 1000.0 1715 real (kind=kind_phys), parameter :: sqvisc = 258.2 1716 1717 real (kind=kind_phys), parameter :: ric = 0.183 1718 real (kind=kind_phys), parameter :: rric = 1.0/ric 1719 real (kind=kind_phys), parameter :: fhneu = 0.8 1720 real (kind=kind_phys), parameter :: rfc = 0.191 1721 real (kind=kind_phys), parameter :: rfac = ric/(fhneu*rfc*rfc) 1722 1723 ! --- inputs: 1724 ! real (kind=kind_phys), intent(in) :: zlvl, z0, t1v, th2v, & 1725 ! & sfcspd, czil 1726 1727 ! --- input/outputs: 1728 ! real (kind=kind_phys), intent(inout) :: cm, ch 1729 1730 ! --- locals: 1731 real (kind=kind_phys) :: zilfc, zu, zt, rdz, cxch, dthv, du2, & 1732 & btgh, wstar2, ustar, zslu, zslt, rlogu, rlogt, rlmo, & 1733 & zetalt, zetalu, zetau, zetat, xlu4, xlt4, xu4, xt4, & 1734 & xlu, xlt, xu, xt, psmz, simm, pshz, simh, ustark, & 1735 & rlmn, rlma 1736 1737 integer :: ilech, itr 1738 1739 ! --- define local in-line functions: 1740 1741 real (kind=kind_phys) :: pslmu, pslms, pslhu, pslhs, zz 1742 real (kind=kind_phys) :: pspmu, pspms, psphu, psphs, xx, yy 1743 1744 ! ... 1) lech's surface functions 1745 1746 pslmu( zz ) = -0.96 * log( 1.0-4.5*zz ) 1747 pslms( zz ) = zz*rric - 2.076*(1.0 - 1.0/(zz + 1.0)) 1748 pslhu( zz ) = -0.96 * log( 1.0-4.5*zz ) 1749 pslhs( zz ) = zz*rfac - 2.076*(1.0 - 1.0/(zz + 1.0)) 1750 1751 ! ... 2) paulson's surface functions 1752 1753 pspmu( xx ) = -2.0 * log( (xx + 1.0)*0.5 ) & 1754 & - log( (xx*xx + 1.0)*0.5 ) + 2.0*atan(xx) - pihf 1755 pspms( yy ) = 5.0 * yy 1756 psphu( xx ) = -2.0 * log( (xx*xx + 1.0)*0.5 ) 1757 psphs( yy ) = 5.0 * yy 1758 1759 ! 1760 !===> ... begin here 1761 ! 1762 ! --- ... this routine sfcdif can handle both over open water (sea, ocean) and 1763 ! over solid surface (land, sea-ice). 1764 1765 ilech = 0 1766 1767 ! --- ... ztfc: ratio of zoh/zom less or equal than 1 1768 ! czil: constant c in zilitinkevich, s. s.1995,:note about zt 1769 1770 zilfc = -czil * vkrm * sqvisc 1771 1772 zu = z0 1773 1774 rdz = 1.0 / zlvl 1775 cxch = excm * rdz 1776 dthv = th2v - t1v 1777 du2 = max( sfcspd*sfcspd, epsu2 ) 1778 1779 ! --- ... beljars correction of ustar 1780 1781 btgh = btg * hpbl 1782 1783 ! --- ... if statements to avoid tangent linear problems near zero 1784 if (btgh*ch*dthv /= 0.0) then 1785 wstar2 = wwst2 * abs( btgh*ch*dthv )**(2.0/3.0) 1786 else 1787 wstar2 = 0.0 1788 endif 1789 1790 ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust ) 1791 1792 ! --- ... zilitinkevitch approach for zt 1793 1794 zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0 1795 1796 zslu = zlvl + zu 1797 zslt = zlvl + zt 1798 1799 ! print*,'zslt=',zslt 1800 ! print*,'zlvl=',zvll 1801 ! print*,'zt=',zt 1802 1803 rlogu = log( zslu/zu ) 1804 rlogt = log( zslt/zt ) 1805 1806 rlmo = elfc*ch*dthv / ustar**3 1807 1808 ! print*,'rlmo=',rlmo 1809 ! print*,'elfc=',elfc 1810 ! print*,'ch=',ch 1811 ! print*,'dthv=',dthv 1812 ! print*,'ustar=',ustar 1813 1814 do itr = 1, itrmx 1815 1816 ! --- ... 1./ monin-obukkhov length-scale 1817 1818 zetalt = max( zslt*rlmo, ztmin ) 1819 rlmo = zetalt / zslt 1820 zetalu = zslu * rlmo 1821 zetau = zu * rlmo 1822 zetat = zt * rlmo 1823 1824 if (ilech == 0) then 1825 1826 if (rlmo < 0.0) then 1827 xlu4 = 1.0 - 16.0 * zetalu 1828 xlt4 = 1.0 - 16.0 * zetalt 1829 xu4 = 1.0 - 16.0 * zetau 1830 xt4 = 1.0 - 16.0* zetat 1831 1832 xlu = sqrt( sqrt( xlu4 ) ) 1833 xlt = sqrt( sqrt( xlt4 ) ) 1834 xu = sqrt( sqrt( xu4 ) ) 1835 xt = sqrt( sqrt( xt4 ) ) 1836 1837 psmz = pspmu(xu) 1838 1839 ! print*,'-----------1------------' 1840 ! print*,'psmz=',psmz 1841 ! print*,'pspmu(zetau)=',pspmu( zetau ) 1842 ! print*,'xu=',xu 1843 ! print*,'------------------------' 1844 1845 simm = pspmu( xlu ) - psmz + rlogu 1846 pshz = psphu( xt ) 1847 simh = psphu( xlt ) - pshz + rlogt 1848 else 1849 zetalu = min( zetalu, ztmax ) 1850 zetalt = min( zetalt, ztmax ) 1851 psmz = pspms( zetau ) 1852 1853 ! print*,'-----------2------------' 1854 ! print*,'psmz=',psmz 1855 ! print*,'pspms(zetau)=',pspms( zetau ) 1856 ! print*,'zetau=',zetau 1857 ! print*,'------------------------' 1858 1859 simm = pspms( zetalu ) - psmz + rlogu 1860 pshz = psphs( zetat ) 1861 simh = psphs( zetalt ) - pshz + rlogt 1862 endif ! end if_rlmo_block 1863 1864 else 1865 1866 ! --- ... lech's functions 1867 1868 if (rlmo < 0.0) then 1869 psmz = pslmu( zetau ) 1870 1871 ! print*,'-----------3------------' 1872 ! print*,'psmz=',psmz 1873 ! print*,'pslmu(zetau)=',pslmu( zetau ) 1874 ! print*,'zetau=',zetau 1875 ! print*,'------------------------' 1876 1877 simm = pslmu( zetalu ) - psmz + rlogu 1878 pshz = pslhu( zetat ) 1879 simh = pslhu( zetalt ) - pshz + rlogt 1880 else 1881 zetalu = min( zetalu, ztmax ) 1882 zetalt = min( zetalt, ztmax ) 1883 1884 psmz = pslms( zetau ) 1885 1886 ! print*,'-----------4------------' 1887 ! print*,'psmz=',psmz 1888 ! print*,'pslms(zetau)=',pslms( zetau ) 1889 ! print*,'zetau=',zetau 1890 ! print*,'------------------------' 1891 1892 simm = pslms( zetalu ) - psmz + rlogu 1893 pshz = pslhs( zetat ) 1894 simh = pslhs( zetalt ) - pshz + rlogt 1895 endif ! end if_rlmo_block 1896 1897 endif ! end if_ilech_block 1898 1899 ! --- ... beljaars correction for ustar 1900 1901 ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust ) 1902 1903 ! --- ... zilitinkevitch fix for zt 1904 1905 zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0 1906 1907 zslt = zlvl + zt 1908 rlogt = log( zslt/zt ) 1909 1910 ustark = ustar * vkrm 1911 cm = max( ustark/simm, cxch ) 1912 ch = max( ustark/simh, cxch ) 1913 1914 ! --- ... if statements to avoid tangent linear problems near zero 1915 1916 if (btgh*ch*dthv /= 0.0) then 1917 wstar2 = wwst2 * abs(btgh*ch*dthv) ** (2.0/3.0) 1918 else 1919 wstar2 = 0.0 1920 endif 1921 1922 rlmn = elfc*ch*dthv / ustar**3 1923 rlma = rlmo*wold + rlmn*wnew 1924 1925 rlmo = rlma 1926 1927 enddo ! end do_itr_loop 1928 1929 ! print*,'----------------------------' 1930 ! print*,'sfcdif output ! ! ! ! ! ! ! ! ! ! ! !' 1931 ! 1932 ! print*,'zlvl=',zlvl 1933 ! print*,'z0=',z0 1934 ! print*,'t1v=',t1v 1935 ! print*,'th2v=',th2v 1936 ! print*,'sfcspd=',sfcspd 1937 ! print*,'czil=',czil 1938 ! print*,'cm=',cm 1939 ! print*,'ch=',ch 1940 ! print*,'----------------------------' 1941 ! 1942 return 1943 !................................... 1944 end subroutine sfcdif 1945 !----------------------------------- 1946 1947 1948 !----------------------------------- 1949 subroutine snfrac 1950 !................................... 1951 ! --- inputs: 1952 ! & ( sneqv, snup, salp, snowh, & 1953 ! --- outputs: 1954 ! & sncovr & 1955 ! & ) 1956 1957 ! ===================================================================== ! 1958 ! description: ! 1959 ! ! 1960 ! subroutine snfrac calculatexsnow fraction (0 -> 1) ! 1961 ! ! 1962 ! subprogram called: none ! 1963 ! ! 1964 ! ! 1965 ! ==================== defination of variables ==================== ! 1966 ! ! 1967 ! inputs from the calling program: size ! 1968 ! sneqv - real, snow water equivalent (m) 1 ! 1969 ! snup - real, threshold sneqv depth above which sncovr=1 1 ! 1970 ! salp - real, tuning parameter 1 ! 1971 ! snowh - real, snow depth (m) 1 ! 1972 ! ! 1973 ! outputs to the calling program: ! 1974 ! sncovr - real, fractional snow cover 1 ! 1975 ! ! 1976 ! ==================== end of description ===================== ! 1977 ! 1978 ! --- inputs: 1979 ! real (kind=kind_phys), intent(in) :: sneqv, snup, salp, snowh 1980 1981 ! --- outputs: 1982 ! real (kind=kind_phys), intent(out) :: sncovr 1983 1984 ! --- locals: 1985 real (kind=kind_phys) :: rsnow, z0n 1986 1987 ! 1988 !===> ... begin here 1989 ! 1990 ! --- ... snup is veg-class dependent snowdepth threshhold (set in routine 1991 ! redprm) above which snocvr=1. 1992 1993 if (sneqv < snup) then 1994 rsnow = sneqv / snup 1995 sncovr = 1.0 - (exp(-salp*rsnow) - rsnow*exp(-salp)) 1996 else 1997 sncovr = 1.0 1998 endif 1999 2000 z0n = 0.035 2001 2002 ! --- ... formulation of dickinson et al. 1986 2003 2004 ! sncovr = snowh / (snowh + 5.0*z0n) 2005 2006 ! --- ... formulation of marshall et al. 1994 2007 2008 ! sncovr = sneqv / (sneqv + 2.0*z0n) 2009 2010 ! 2011 return 2012 !................................... 2013 end subroutine snfrac 2014 !----------------------------------- 2015 2016 2017 !----------------------------------- 2018 subroutine snopac 2019 !................................... 2020 ! --- inputs: 2021 ! & ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, smcdry, & 2022 ! & cmcmax, dt, df1, sfcems, sfctmp, t24, th2, fdown, epsca, & 2023 ! & bexp, pc, rch, rr, cfactr, slope, kdt, frzx, psisat, & 2024 ! & zsoil, dwsat, dksat, zbot, shdfac, ice, rtdis, quartz, & 2025 ! & fxexp, csoil, flx2, snowng, & 2026 ! --- input/outputs: 2027 ! & prcp1, cmc, t1, stc, sncovr, sneqv, sndens, snowh, & 2028 ! & sh2o, tbot, beta, & 2029 ! --- outputs: 2030 ! & smc, ssoil, runoff1, runoff2, runoff3, edir, ec, et, & 2031 ! & ett, snomlt, drip, dew, flx1, flx3, esnow & 2032 ! & ) 2033 2034 ! ===================================================================== ! 2035 ! description: ! 2036 ! ! 2037 ! subroutine snopac calculates soil moisture and heat flux values and ! 2038 ! update soil moisture content and soil heat content values for the ! 2039 ! case when a snow pack is present. ! 2040 ! ! 2041 ! ! 2042 ! subprograms called: evapo, smflx, shflx, snowpack 2043 ! ! 2044 ! ==================== defination of variables ==================== ! 2045 ! ! 2046 ! inputs from the calling program: size ! 2047 ! nsoil - integer, number of soil layers 1 ! 2048 ! nroot - integer, number of root layers 1 ! 2049 ! etp - real, potential evaporation 1 ! 2050 ! prcp - real, precip rate 1 ! 2051 ! smcmax - real, porosity 1 ! 2052 ! smcwlt - real, wilting point 1 ! 2053 ! smcref - real, soil mois threshold 1 ! 2054 ! smcdry - real, dry soil mois threshold 1 ! 2055 ! cmcmax - real, maximum canopy water parameters 1 ! 2056 ! dt - real, time step 1 ! 2057 ! df1 - real, thermal diffusivity m ! 2058 ! sfcems - real, lw surface emissivity 1 ! 2059 ! sfctmp - real, sfc temperature 1 ! 2060 ! t24 - real, sfctmp**4 1 ! 2061 ! th2 - real, sfc air potential temperature 1 ! 2062 ! fdown - real, net solar + downward lw flux at sfc 1 ! 2063 ! epsca - real, 1 ! 2064 ! bexp - real, soil type "b" parameter 1 ! 2065 ! pc - real, plant coeff 1 ! 2066 ! rch - real, companion coefficient of ch 1 ! 2067 ! rr - real, 1 ! 2068 ! cfactr - real, canopy water parameters 1 ! 2069 ! slope - real, linear reservoir coefficient 1 ! 2070 ! kdt - real, 1 ! 2071 ! frzx - real, frozen ground parameter 1 ! 2072 ! psisat - real, saturated soil potential 1 ! 2073 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 2074 ! dwsat - real, saturated soil diffusivity 1 ! 2075 ! dksat - real, saturated soil hydraulic conductivity 1 ! 2076 ! zbot - real, specify depth of lower bd soil 1 ! 2077 ! shdfac - real, aeral coverage of green vegetation 1 ! 2078 ! ice - integer, sea-ice flag (=1: sea-ice, =0: land) 1 ! 2079 ! rtdis - real, root distribution nsoil ! 2080 ! quartz - real, soil quartz content 1 ! 2081 ! fxexp - real, bare soil evaporation exponent 1 ! 2082 ! csoil - real, soil heat capacity 1 ! 2083 ! flx2 - real, freezing rain latent heat flux 1 ! 2084 ! snowng - logical, snow flag 1 ! 2085 ! ! 2086 ! input/outputs from and to the calling program: ! 2087 ! prcp1 - real, effective precip 1 ! 2088 ! cmc - real, canopy moisture content 1 ! 2089 ! t1 - real, ground/canopy/snowpack eff skin temp 1 ! 2090 ! stc - real, soil temperature nsoil ! 2091 ! sncovr - real, snow cover 1 ! 2092 ! sneqv - real, water-equivalent snow depth 1 ! 2093 ! sndens - real, snow density 1 ! 2094 ! snowh - real, snow depth 1 ! 2095 ! sh2o - real, unfrozen soil moisture nsoil ! 2096 ! tbot - real, bottom soil temperature 1 ! 2097 ! beta - real, ratio of actual/potential evap 1 ! 2098 ! ! 2099 ! outputs to the calling program: ! 2100 ! smc - real, total soil moisture nsoil ! 2101 ! ssoil - real, upward soil heat flux 1 ! 2102 ! runoff1 - real, surface runoff not infiltrating sfc 1 ! 2103 ! runoff2 - real, sub surface runoff 1 ! 2104 ! runoff3 - real, excess of porosity for a given soil layer 1 ! 2105 ! edir - real, direct soil evaporation 1 ! 2106 ! ec - real, canopy water evaporation 1 ! 2107 ! et - real, plant transpiration nsoil ! 2108 ! ett - real, total plant transpiration 1 ! 2109 ! snomlt - real, snow melt water equivalent 1 ! 2110 ! drip - real, through-fall of precip 1 ! 2111 ! dew - real, dewfall (or frostfall) 1 ! 2112 ! flx1 - real, precip-snow sfc flux 1 ! 2113 ! flx3 - real, phase-change heat flux from snowmelt 1 ! 2114 ! esnow - real, sublimation from snowpack 1 ! 2115 ! ! 2116 ! ==================== end of description ===================== ! 2117 ! 2118 ! --- constant parameters: 2119 real, parameter :: esdmin = 1.e-6 2120 2121 ! --- inputs: 2122 ! integer, intent(in) :: nsoil, nroot, ice 2123 2124 ! real (kind=kind_phys), intent(in) :: etp, prcp, smcmax, smcref, & 2125 ! & smcwlt, smcdry, cmcmax, dt, df1, sfcems, sfctmp, t24, & 2126 ! & th2, fdown, epsca, bexp, pc, rch, rr, cfactr, slope, kdt, & 2127 ! & frzx, psisat, dwsat, dksat, zbot, shdfac, quartz, & 2128 ! & csoil, fxexp, flx2, zsoil(nsoil), rtdis(nsoil) 2129 2130 ! logical, intent(in) :: snowng 2131 2132 ! --- input/outputs: 2133 ! real (kind=kind_phys), intent(inout) :: prcp1, t1, sncovr, sneqv, & 2134 ! & sndens, snowh, cmc, tbot, beta, sh2o(nsoil), stc(nsoil) 2135 2136 ! --- outputs: 2137 ! real (kind=kind_phys), intent(out) :: ssoil, runoff1, runoff2, & 2138 ! & runoff3, edir, ec, et(nsoil), ett, snomlt, drip, dew, & 2139 ! & flx1, flx3, esnow, smc(nsoil) 2140 2141 ! --- locals: 2142 real (kind=kind_phys):: denom, dsoil, dtot, eta, etp1, ssoil1, & 2143 & snoexp, ex, t11, t12, t12a, t12b, yy, zz1, seh, t14, & 2144 & ec1, edir1, ett1, etns, etns1, esnow1, esnow2, etanrg, & 2145 & et1(nsoil) 2146 2147 integer k 2148 2149 ! data snoexp /1.0/ !!! <----- for noah v2.7 2150 data snoexp /2.0/ !!! <----- for noah v2.7.1 2151 2152 ! --- ... convert potential evap (etp) from kg m-2 s-1 to m s-1 and then to an 2153 ! amount (m) given timestep (dt) and call it an effective snowpack 2154 ! reduction amount, esnow2 (m) for a snowcover fraction = 1.0. this is 2155 ! the amount the snowpack would be reduced due to sublimation from the 2156 ! snow sfc during the timestep. sublimation will proceed at the 2157 ! potential rate unless the snow depth is less than the expected 2158 ! snowpack reduction. for snowcover fraction = 1.0, 0=edir=et=ec, and 2159 ! hence total evap = esnow = sublimation (potential evap rate) 2160 2161 ! --- ... if sea-ice (ice=1) or glacial-ice (ice=-1), snowcover fraction = 1.0, 2162 ! and sublimation is at the potential rate. 2163 ! for non-glacial land (ice=0), if snowcover fraction < 1.0, total 2164 ! evaporation < potential due to non-potential contribution from 2165 ! non-snow covered fraction. 2166 2167 prcp1 = prcp1 * 0.001 2168 2169 edir = 0.0 2170 edir1 = 0.0 2171 2172 ec = 0.0 2173 ec1 = 0.0 2174 2175 do k = 1, nsoil 2176 et (k) = 0.0 2177 et1(k) = 0.0 2178 enddo 2179 2180 ett = 0.0 2181 ett1 = 0.0 2182 etns = 0.0 2183 etns1 = 0.0 2184 esnow = 0.0 2185 esnow1= 0.0 2186 esnow2= 0.0 2187 2188 dew = 0.0 2189 etp1 = etp * 0.001 2190 2191 if (etp < 0.0) then 2192 2193 ! --- ... if etp<0 (downward) then dewfall (=frostfall in this case). 2194 2195 dew = -etp1 2196 esnow2 = etp1 * dt 2197 etanrg = etp * ((1.0-sncovr)*lsubc + sncovr*lsubs) 2198 2199 else 2200 2201 ! --- ... etp >= 0, upward moisture flux 2202 2203 if (ice /= 0) then ! for sea-ice and glacial-ice case 2204 2205 esnow = etp 2206 esnow1 = esnow * 0.001 2207 esnow2 = esnow1 * dt 2208 etanrg = esnow * lsubs 2209 2210 else ! for non-glacial land case 2211 2212 if (sncovr < 1.0) then 2213 2214 call evapo & 2215 ! --- inputs: 2216 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, & 2217 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, & 2218 & shdfac, cfactr, rtdis, fxexp, & 2219 ! --- outputs: 2220 & etns1, edir1, ec1, et1, ett1 & 2221 & ) 2222 2223 edir1 = edir1 * (1.0 - sncovr) 2224 ec1 = ec1 * (1.0 - sncovr) 2225 2226 do k = 1, nsoil 2227 et1(k) = et1(k) * (1.0 - sncovr) 2228 enddo 2229 2230 ett1 = ett1 * (1.0 - sncovr) 2231 etns1 = etns1 * (1.0 - sncovr) 2232 2233 edir = edir1 * 1000.0 2234 ec = ec1 * 1000.0 2235 2236 do k = 1, nsoil 2237 et(k) = et1(k) * 1000.0 2238 enddo 2239 2240 ett = ett1 * 1000.0 2241 etns = etns1 * 1000.0 2242 2243 endif ! end if_sncovr_block 2244 2245 esnow = etp * sncovr 2246 ! esnow1 = etp * 0.001 2247 esnow1 = esnow * 0.001 2248 esnow2 = esnow1 * dt 2249 etanrg = esnow*lsubs + etns*lsubc 2250 2251 endif ! end if_ice_block 2252 2253 endif ! end if_etp_block 2254 2255 ! --- ... if precip is falling, calculate heat flux from snow sfc to newly 2256 ! accumulating precip. note that this reflects the flux appropriate for 2257 ! the not-yet-updated skin temperature (t1). assumes temperature of the 2258 ! snowfall striking the gound is =sfctmp (lowest model level air temp). 2259 2260 flx1 = 0.0 2261 if ( snowng ) then 2262 flx1 = cpice * prcp * (t1 - sfctmp) 2263 else 2264 if (prcp > 0.0) flx1 = cph2o1 * prcp * (t1 - sfctmp) 2265 endif 2266 2267 ! --- ... calculate an 'effective snow-grnd sfc temp' (t12) based on heat 2268 ! fluxes between the snow pack and the soil and on net radiation. 2269 ! include flx1 (precip-snow sfc) and flx2 (freezing rain latent 2270 ! heat) fluxes. 2271 ! flx2 reflects freezing rain latent heat flux using t1 calculated 2272 ! in penman. 2273 2274 dsoil = -0.5 * zsoil(1) 2275 dtot = snowh + dsoil 2276 denom = 1.0 + df1 / (dtot * rr * rch) 2277 2278 ! t12a = ( (fdown - flx1 - flx2 - sigma1*t24) / rch & 2279 ! & + th2 - sfctmp - beta*epsca ) / rr 2280 t12a = ( (fdown - flx1 - flx2 - sfcems*sigma1*t24) / rch & 2281 & + th2 - sfctmp - etanrg/rch ) / rr 2282 2283 t12b = df1 * stc(1) / (dtot * rr * rch) 2284 t12 = (sfctmp + t12a + t12b) / denom 2285 2286 ! --- ... if the 'effective snow-grnd sfc temp' is at or below freezing, no snow 2287 ! melt will occur. set the skin temp to this effective temp. reduce 2288 ! (by sublimination ) or increase (by frost) the depth of the snowpack, 2289 ! depending on sign of etp. 2290 ! update soil heat flux (ssoil) using new skin temperature (t1) 2291 ! since no snowmelt, set accumulated snowmelt to zero, set 'effective' 2292 ! precip from snowmelt to zero, set phase-change heat flux from snowmelt 2293 ! to zero. 2294 2295 if (t12 <= tfreez) then 2296 2297 t1 = t12 2298 ssoil = df1 * (t1 - stc(1)) / dtot 2299 sneqv = max(0.0, sneqv-esnow2) 2300 flx3 = 0.0 2301 ex = 0.0 2302 snomlt = 0.0 2303 2304 else 2305 2306 ! --- ... if the 'effective snow-grnd sfc temp' is above freezing, snow melt 2307 ! will occur. call the snow melt rate,ex and amt, snomlt. revise the 2308 ! effective snow depth. revise the skin temp because it would have chgd 2309 ! due to the latent heat released by the melting. calc the latent heat 2310 ! released, flx3. set the effective precip, prcp1 to the snow melt rate, 2311 ! ex for use in smflx. adjustment to t1 to account for snow patches. 2312 ! calculate qsat valid at freezing point. note that esat (saturation 2313 ! vapor pressure) value of 6.11e+2 used here is that valid at frzzing 2314 ! point. note that etp from call penman in sflx is ignored here in 2315 ! favor of bulk etp over 'open water' at freezing temp. 2316 ! update soil heat flux (s) using new skin temperature (t1) 2317 2318 ! --- ... noah v2.7.1 mek feb2004 2319 ! non-linear weighting of snow vs non-snow covered portions of gridbox 2320 ! so with snoexp = 2.0 (>1), surface skin temperature is higher than 2321 ! for the linear case (snoexp = 1). 2322 2323 t1 = tfreez * sncovr**snoexp + t12 * (1.0 - sncovr**snoexp) 2324 2325 beta = 1.0 2326 ssoil = df1 * (t1 - stc(1)) / dtot 2327 2328 ! --- ... if potential evap (sublimation) greater than depth of snowpack. 2329 ! beta<1 2330 ! snowpack has sublimated away, set depth to zero. 2331 2332 if (sneqv-esnow2 <= esdmin) then 2333 2334 sneqv = 0.0 2335 ex = 0.0 2336 snomlt = 0.0 2337 flx3 = 0.0 2338 2339 else 2340 2341 ! --- ... potential evap (sublimation) less than depth of snowpack, retain 2342 ! beta=1. 2343 2344 sneqv = sneqv - esnow2 2345 seh = rch * (t1 - th2) 2346 2347 t14 = t1 * t1 2348 t14 = t14 * t14 2349 2350 flx3 = fdown - flx1 - flx2 - sfcems*sigma1*t14 & 2351 & - ssoil - seh - etanrg 2352 if (flx3 <= 0.0) flx3 = 0.0 2353 2354 ex = flx3 * 0.001 / lsubf 2355 2356 ! --- ... snowmelt reduction depending on snow cover 2357 ! if snow cover less than 5% no snowmelt reduction 2358 ! note: does 'if' below fail to match the melt water with the melt 2359 ! energy? 2360 2361 ! if (sncovr > 0.05) ex = ex * sncovr 2362 snomlt = ex * dt 2363 2364 ! --- ... esdmin represents a snowpack depth threshold value below which we 2365 ! choose not to retain any snowpack, and instead include it in snowmelt. 2366 2367 if (sneqv-snomlt >= esdmin) then 2368 2369 sneqv = sneqv - snomlt 2370 2371 else 2372 2373 ! --- ... snowmelt exceeds snow depth 2374 2375 ex = sneqv / dt 2376 flx3 = ex * 1000.0 * lsubf 2377 snomlt = sneqv 2378 sneqv = 0.0 2379 2380 endif ! end if_sneqv-snomlt_block 2381 2382 endif ! end if_sneqv-esnow2_block 2383 2384 ! prcp1 = prcp1 + ex 2385 2386 ! --- ... if non-glacial land, add snowmelt rate (ex) to precip rate to be used 2387 ! in subroutine smflx (soil moisture evolution) via infiltration. 2388 2389 ! --- ... for sea-ice and glacial-ice, the snowmelt will be added to subsurface 2390 ! runoff/baseflow later near the end of sflx (after return from call to 2391 ! subroutine snopac) 2392 2393 if (ice == 0) prcp1 = prcp1 + ex 2394 2395 endif ! end if_t12<=tfreez_block 2396 2397 ! --- ... final beta now in hand, so compute evaporation. evap equals etp 2398 ! unless beta<1. 2399 2400 ! eta = beta * etp 2401 2402 ! --- ... smflx returns updated soil moisture values for non-glacial land. 2403 ! if sea-ice (ice=1) or glacial-ice (ice=-1), skip call to smflx, since 2404 ! no soil medium for sea-ice or glacial-ice 2405 2406 if (ice == 0) then 2407 2408 call smflx & 2409 ! --- inputs: 2410 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, & 2411 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, & 2412 & edir1, ec1, et1, & 2413 ! --- input/outputs: 2414 & cmc, sh2o, & 2415 ! --- outputs: 2416 & smc, runoff1, runoff2, runoff3, drip & 2417 & ) 2418 2419 endif 2420 2421 ! --- ... before call shflx in this snowpack case, set zz1 and yy arguments to 2422 ! special values that ensure that ground heat flux calculated in shflx 2423 ! matches that already computer for below the snowpack, thus the sfc 2424 ! heat flux to be computed in shflx will effectively be the flux at the 2425 ! snow top surface. t11 is a dummy arguement so we will not use the 2426 ! skin temp value as revised by shflx. 2427 2428 zz1 = 1.0 2429 yy = stc(1) - 0.5*ssoil*zsoil(1)*zz1 / df1 2430 t11 = t1 2431 2432 ! --- ... shflx will calc/update the soil temps. note: the sub-sfc heat flux 2433 ! (ssoil1) and the skin temp (t11) output from this shflx call are not 2434 ! used in any subsequent calculations. rather, they are dummy variables 2435 ! here in the snopac case, since the skin temp and sub-sfc heat flux are 2436 ! updated instead near the beginning of the call to snopac. 2437 2438 call shflx & 2439 ! --- inputs: 2440 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, & 2441 & psisat, bexp, df1, ice, quartz, csoil, & 2442 ! --- input/outputs: 2443 & stc, t11, tbot, sh2o, & 2444 ! --- outputs: 2445 & ssoil1 & 2446 & ) 2447 2448 ! --- ... snow depth and density adjustment based on snow compaction. yy is 2449 ! assumed to be the soil temperture at the top of the soil column. 2450 2451 if (ice == 0) then ! for non-glacial land 2452 2453 if (sneqv > 0.0) then 2454 2455 call snowpack & 2456 ! --- inputs: 2457 & ( sneqv, dt, t1, yy, & 2458 ! --- input/outputs: 2459 & snowh, sndens & 2460 & ) 2461 2462 else 2463 2464 sneqv = 0.0 2465 snowh = 0.0 2466 sndens = 0.0 2467 ! sncond = 1.0 2468 sncovr = 0.0 2469 2470 endif ! end if_sneqv_block 2471 2472 ! --- ... over sea-ice or glacial-ice, if s.w.e. (sneqv) below threshold lower 2473 ! bound (0.01 m for sea-ice, 0.10 m for glacial-ice), then set at 2474 ! lower bound and store the source increment in subsurface runoff/ 2475 ! baseflow (runoff2). note: runoff2 is then a negative value (as 2476 ! a flag) over sea-ice or glacial-ice, in order to achieve water balance. 2477 2478 elseif (ice == 1) then ! for sea-ice 2479 2480 if (sneqv >= 0.01) then 2481 2482 call snowpack & 2483 ! --- inputs: 2484 & ( sneqv, dt, t1, yy, & 2485 ! --- input/outputs: 2486 & snowh, sndens & 2487 & ) 2488 2489 else 2490 2491 ! sndens = sneqv / snowh 2492 ! runoff2 = -(0.01 - sneqv) / dt 2493 sneqv = 0.01 2494 snowh = 0.05 2495 sncovr = 1.0 2496 ! snowh = sneqv / sndens 2497 2498 endif ! end if_sneqv_block 2499 2500 else ! for glacial-ice 2501 2502 if (sneqv >= 0.10) then 2503 2504 call snowpack & 2505 ! --- inputs: 2506 & ( sneqv, dt, t1, yy, & 2507 ! --- input/outputs: 2508 & snowh, sndens & 2509 & ) 2510 2511 else 2512 2513 ! sndens = sneqv / snowh 2514 ! runoff2 = -(0.10 - sneqv) / dt 2515 sneqv = 0.10 2516 snowh = 0.50 2517 sncovr = 1.0 2518 ! snowh = sneqv / sndens 2519 2520 endif ! end if_sneqv_block 2521 2522 endif ! end if_ice_block 2523 2524 ! 2525 return 2526 !................................... 2527 end subroutine snopac 2528 !----------------------------------- 2529 2530 2531 !----------------------------------- 2532 subroutine snow_new 2533 !................................... 2534 ! --- inputs: 2535 ! & ( sfctmp, sn_new, & 2536 ! --- input/outputs: 2537 ! & snowh, sndens & 2538 ! & ) 2539 2540 ! ===================================================================== ! 2541 ! description: ! 2542 ! ! 2543 ! subroutine snow_new calculates snow depth and densitity to account ! 2544 ! for the new snowfall. new values of snow depth & density returned. ! 2545 ! ! 2546 ! subprogram called: none ! 2547 ! ! 2548 ! ==================== defination of variables ==================== ! 2549 ! ! 2550 ! inputs from the calling program: size ! 2551 ! sfctmp - real, surface air temperature (k) 1 ! 2552 ! sn_new - real, new snowfall (m) 1 ! 2553 ! ! 2554 ! input/outputs from and to the calling program: ! 2555 ! snowh - real, snow depth (m) 1 ! 2556 ! sndens - real, snow density 1 ! 2557 ! (g/cm3=dimensionless fraction of h2o density) ! 2558 ! ! 2559 ! ==================== end of description ===================== ! 2560 ! 2561 ! --- inputs: 2562 ! real(kind=kind_phys), intent(in) :: sfctmp, sn_new 2563 2564 ! --- input/outputs: 2565 ! real(kind=kind_phys), intent(inout) :: snowh, sndens 2566 2567 ! --- locals: 2568 real(kind=kind_phys) :: dsnew, snowhc, hnewc, newsnc, tempc 2569 2570 ! 2571 !===> ... begin here 2572 ! 2573 ! --- ... conversion into simulation units 2574 2575 snowhc = snowh * 100.0 2576 newsnc = sn_new * 100.0 2577 tempc = sfctmp - tfreez 2578 2579 ! --- ... calculating new snowfall density depending on temperature 2580 ! equation from gottlib l. 'a general runoff model for 2581 ! snowcovered and glacierized basin', 6th nordic hydrological 2582 ! conference, vemadolen, sweden, 1980, 172-177pp. 2583 2584 if (tempc <= -15.0) then 2585 dsnew = 0.05 2586 else 2587 dsnew = 0.05 + 0.0017*(tempc + 15.0)**1.5 2588 endif 2589 2590 ! --- ... adjustment of snow density depending on new snowfall 2591 2592 hnewc = newsnc / dsnew 2593 sndens = (snowhc*sndens + hnewc*dsnew) / (snowhc + hnewc) 2594 snowhc = snowhc + hnewc 2595 snowh = snowhc * 0.01 2596 ! 2597 return 2598 !................................... 2599 end subroutine snow_new 2600 !----------------------------------- 2601 2602 2603 !----------------------------------- 2604 subroutine snowz0 2605 !................................... 2606 ! --- inputs: 2607 ! & ( sncovr, & 2608 ! --- input/outputs: 2609 ! & z0 & 2610 ! & ) 2611 2612 ! ===================================================================== ! 2613 ! description: ! 2614 ! ! 2615 ! subroutine snowz0 calculates total roughness length over snow ! 2616 ! ! 2617 ! subprogram called: none ! 2618 ! ! 2619 ! ==================== defination of variables ==================== ! 2620 ! ! 2621 ! inputs from the calling program: size ! 2622 ! sncovr - real, fractional snow cover 1 ! 2623 ! ! 2624 ! input/outputs from and to the calling program: ! 2625 ! z0 - real, roughness length (m) 1 ! 2626 ! ! 2627 ! ==================== end of description ===================== ! 2628 ! 2629 ! --- inputs: 2630 ! real(kind=kind_phys), intent(in) :: sncovr 2631 2632 ! --- input/outputs: 2633 ! real(kind=kind_phys), intent(inout) :: z0 2634 2635 ! --- locals: 2636 real(kind=kind_phys) :: z0s 2637 ! 2638 !===> ... begin here 2639 ! 2640 ! z0s = 0.001 ! snow roughness length:=0.001 (m) 2641 ! --- ... current noah lsm condition - mbek, 09-oct-2001 2642 z0s = z0 2643 2644 z0 = (1.0 - sncovr)*z0 + sncovr*z0s 2645 2646 ! 2647 return 2648 !................................... 2649 end subroutine snowz0 2650 !----------------------------------- 2651 2652 2653 !----------------------------------- 2654 subroutine tdfcnd & 2655 !................................... 2656 ! --- inputs: 2657 & ( smc, qz, smcmax, sh2o, & 2658 ! --- outputs: 2659 & df & 2660 & ) 2661 2662 ! ===================================================================== ! 2663 ! description: ! 2664 ! ! 2665 ! subroutine tdfcnd calculates thermal diffusivity and conductivity ! 2666 ! of the soil for a given point and time. ! 2667 ! ! 2668 ! peters-lidard approach (peters-lidard et al., 1998) ! 2669 ! june 2001 changes: frozen soil condition. ! 2670 ! ! 2671 ! subprogram called: none ! 2672 ! ! 2673 ! use as in peters-lidard, 1998 (modif. from johansen, 1975). ! 2674 ! pablo grunmann, 08/17/98 ! 2675 ! refs.: ! 2676 ! farouki, o.t.,1986: thermal properties of soils. series on rock ! 2677 ! and soil mechanics, vol. 11, trans tech, 136 pp. ! 2678 ! johansen, o., 1975: thermal conductivity of soils. ph.d. thesis, ! 2679 ! university of trondheim, ! 2680 ! peters-lidard, c. d., et al., 1998: the effect of soil thermal ! 2681 ! conductivity parameterization on surface energy fluxes ! 2682 ! and temperatures. journal of the atmospheric sciences, ! 2683 ! vol. 55, pp. 1209-1224. ! 2684 ! ! 2685 ! ==================== defination of variables ==================== ! 2686 ! ! 2687 ! inputs: size ! 2688 ! smc - real, top layer total soil moisture 1 ! 2689 ! qz - real, quartz content (soil type dependent) 1 ! 2690 ! smcmax - real, porosity 1 ! 2691 ! sh2o - real, top layer unfrozen soil moisture 1 ! 2692 ! ! 2693 ! outputs: ! 2694 ! df - real, soil thermal diffusivity and conductivity 1 ! 2695 ! ! 2696 ! locals: ! 2697 ! thkw - water thermal conductivity 1 ! 2698 ! thkqtz - thermal conductivity for quartz 1 ! 2699 ! thko - thermal conductivity for other soil components 1 ! 2700 ! thkqtz - thermal conductivity for the solids combined 1 ! 2701 ! thkice - ice thermal conductivity 1 ! 2702 ! ! 2703 ! ! 2704 ! ==================== end of description ===================== ! 2705 ! 2706 ! --- input: 2707 real (kind=kind_phys), intent(in) :: smc, qz, smcmax, sh2o 2708 2709 ! --- output: 2710 real (kind=kind_phys), intent(out) :: df 2711 2712 ! --- locals: 2713 real (kind=kind_phys) :: gammd, thkdry, ake, thkice, thko, & 2714 & thkqtz, thksat, thks, thkw, satratio, xu, xunfroz 2715 ! 2716 !===> ... begin here 2717 ! 2718 ! --- ... if the soil has any moisture content compute a partial sum/product 2719 ! otherwise use a constant value which works well with most soils 2720 2721 ! --- ... saturation ratio: 2722 2723 satratio = smc / smcmax 2724 2725 ! --- ... parameters w/(m.k) 2726 thkice = 2.2 2727 thkw = 0.57 2728 thko = 2.0 2729 ! if (qz <= 0.2) thko = 3.0 2730 thkqtz = 7.7 2731 2732 ! --- ... solids' conductivity 2733 2734 thks = (thkqtz**qz) * (thko**(1.0-qz)) 2735 2736 ! --- ... unfrozen fraction (from 1., i.e., 100%liquid, to 0. (100% frozen)) 2737 2738 xunfroz = (sh2o + 1.e-9) / (smc + 1.e-9) 2739 2740 ! --- ... unfrozen volume for saturation (porosity*xunfroz) 2741 2742 xu=xunfroz*smcmax 2743 2744 ! --- ... saturated thermal conductivity 2745 2746 thksat = thks**(1.-smcmax) * thkice**(smcmax-xu) * thkw**(xu) 2747 2748 ! --- ... dry density in kg/m3 2749 2750 gammd = (1.0 - smcmax) * 2700.0 2751 2752 ! --- ... dry thermal conductivity in w.m-1.k-1 2753 2754 thkdry = (0.135*gammd + 64.7) / (2700.0 - 0.947*gammd) 2755 2756 if ( sh2o+0.0005 < smc ) then ! frozen 2757 2758 ake = satratio 2759 2760 else ! unfrozen 2761 2762 ! --- ... range of validity for the kersten number (ake) 2763 if ( satratio > 0.1 ) then 2764 2765 ! --- ... kersten number (using "fine" formula, valid for soils containing 2766 ! at least 5% of particles with diameter less than 2.e-6 meters.) 2767 ! (for "coarse" formula, see peters-lidard et al., 1998). 2768 2769 ake = log10( satratio ) + 1.0 2770 2771 else 2772 2773 ! --- ... use k = kdry 2774 ake = 0.0 2775 2776 endif ! end if_satratio_block 2777 2778 endif ! end if_sh2o+0.0005_block 2779 2780 ! --- ... thermal conductivity 2781 2782 df = ake * (thksat - thkdry) + thkdry 2783 ! 2784 return 2785 !................................... 2786 end subroutine tdfcnd 2787 !----------------------------------- 2788 2789 2790 !*********************************************! 2791 ! section-2 2nd level subprograms ! 2792 !*********************************************! 2793 2794 2795 !----------------------------------- 2796 subroutine evapo & 2797 !................................... 2798 ! --- inputs: 2799 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, & 2800 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, & 2801 & shdfac, cfactr, rtdis, fxexp, & 2802 ! --- outputs: 2803 & eta1, edir1, ec1, et1, ett1 & 2804 & ) 2805 2806 ! ===================================================================== ! 2807 ! description: ! 2808 ! ! 2809 ! subroutine evapo calculates soil moisture flux. the soil moisture ! 2810 ! content (smc - a per unit volume measurement) is a dependent variable! 2811 ! that is updated with prognostic eqns. the canopy moisture content ! 2812 ! (cmc) is also updated. frozen ground version: new states added: ! 2813 ! sh2o, and frozen ground correction factor, frzfact and parameter ! 2814 ! slope. ! 2815 ! ! 2816 ! ! 2817 ! subprogram called: devap, transp ! 2818 ! ! 2819 ! ==================== defination of variables ==================== ! 2820 ! ! 2821 ! inputs from calling program: size ! 2822 ! nsoil - integer, number of soil layers 1 ! 2823 ! nroot - integer, number of root layers 1 ! 2824 ! cmc - real, canopy moisture content 1 ! 2825 ! cmcmax - real, maximum canopy water parameters 1 ! 2826 ! etp1 - real, potential evaporation 1 ! 2827 ! dt - real, time step 1 ! 2828 ! zsoil - real, soil layer depth below ground nsoil ! 2829 ! sh2o - real, unfrozen soil moisture nsoil ! 2830 ! smcmax - real, porosity 1 ! 2831 ! smcwlt - real, wilting point 1 ! 2832 ! smcref - real, soil mois threshold 1 ! 2833 ! smcdry - real, dry soil mois threshold 1 ! 2834 ! pc - real, plant coeff 1 ! 2835 ! cfactr - real, canopy water parameters 1 ! 2836 ! rtdis - real, root distribution nsoil ! 2837 ! fxexp - real, bare soil evaporation exponent 1 ! 2838 ! ! 2839 ! outputs to calling program: ! 2840 ! eta1 - real, latent heat flux 1 ! 2841 ! edir1 - real, direct soil evaporation 1 ! 2842 ! ec1 - real, canopy water evaporation 1 ! 2843 ! et1 - real, plant transpiration nsoil ! 2844 ! ett1 - real, total plant transpiration 1 ! 2845 ! ! 2846 ! ==================== end of description ===================== ! 2847 ! 2848 ! --- inputs: 2849 integer, intent(in) :: nsoil, nroot 2850 2851 real (kind=kind_phys), intent(in) :: cmc, cmcmax, etp1, dt, pc, & 2852 & smcmax, smcwlt, smcref, smcdry, shdfac, cfactr, fxexp, & 2853 & zsoil(nsoil), sh2o(nsoil), rtdis(nsoil) 2854 2855 ! --- outputs: 2856 real (kind=kind_phys), intent(out) :: eta1, edir1, ec1, ett1, & 2857 & et1(nsoil) 2858 2859 ! --- locals: 2860 real (kind=kind_phys) :: cmc2ms 2861 2862 integer :: i, k 2863 2864 ! 2865 !===> ... begin here 2866 ! 2867 ! --- ... executable code begins here if the potential evapotranspiration 2868 ! is greater than zero. 2869 2870 edir1 = 0.0 2871 ec1 = 0.0 2872 2873 do k = 1, nsoil 2874 et1(k) = 0.0 2875 enddo 2876 ett1 = 0.0 2877 2878 if (etp1 > 0.0) then 2879 2880 ! --- ... retrieve direct evaporation from soil surface. call this function 2881 ! only if veg cover not complete. 2882 ! frozen ground version: sh2o states replace smc states. 2883 2884 if (shdfac < 1.0) then 2885 2886 call devap & 2887 ! --- inputs: 2888 & ( etp1, sh2o(1), shdfac, smcmax, smcdry, fxexp, & 2889 ! --- outputs: 2890 & edir1 & 2891 & ) 2892 2893 endif 2894 2895 ! --- ... initialize plant total transpiration, retrieve plant transpiration, 2896 ! and accumulate it for all soil layers. 2897 2898 if (shdfac > 0.0) then 2899 2900 call transp & 2901 ! --- inputs: 2902 & ( nsoil, nroot, etp1, sh2o, smcwlt, smcref, & 2903 & cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis, & 2904 ! --- outputs: 2905 & et1 & 2906 & ) 2907 2908 do k = 1, nsoil 2909 ett1 = ett1 + et1(k) 2910 enddo 2911 2912 ! --- ... calculate canopy evaporation. 2913 ! if statements to avoid tangent linear problems near cmc=0.0. 2914 2915 if (cmc > 0.0) then 2916 ec1 = shdfac * ( (cmc/cmcmax)**cfactr ) * etp1 2917 else 2918 ec1 = 0.0 2919 endif 2920 2921 ! --- ... ec should be limited by the total amount of available water 2922 ! on the canopy. -f.chen, 18-oct-1994 2923 2924 cmc2ms = cmc / dt 2925 ec1 = min ( cmc2ms, ec1 ) 2926 endif 2927 2928 endif ! end if_etp1_block 2929 2930 ! --- ... total up evap and transp types to obtain actual evapotransp 2931 2932 eta1 = edir1 + ett1 + ec1 2933 2934 ! 2935 return 2936 !................................... 2937 end subroutine evapo 2938 !----------------------------------- 2939 2940 2941 !----------------------------------- 2942 subroutine shflx & 2943 !................................... 2944 ! --- inputs: 2945 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, & 2946 & psisat, bexp, df1, ice, quartz, csoil, & 2947 ! --- input/outputs: 2948 & stc, t1, tbot, sh2o, & 2949 ! --- outputs: 2950 & ssoil & 2951 & ) 2952 2953 ! ===================================================================== ! 2954 ! description: ! 2955 ! ! 2956 ! subroutine shflx updates the temperature state of the soil column ! 2957 ! based on the thermal diffusion equation and update the frozen soil ! 2958 ! moisture content based on the temperature. ! 2959 ! ! 2960 ! subprogram called: hstep, hrtice, hrt ! 2961 ! ! 2962 ! ! 2963 ! ==================== defination of variables ==================== ! 2964 ! ! 2965 ! inputs: size ! 2966 ! nsoil - integer, number of soil layers 1 ! 2967 ! smc - real, total soil moisture nsoil ! 2968 ! smcmax - real, porosity (sat val of soil mois) 1 ! 2969 ! dt - real, time step 1 ! 2970 ! yy - real, soil temperature at the top of column 1 ! 2971 ! zz1 - real, 1 ! 2972 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 2973 ! zbot - real, specify depth of lower bd soil 1 ! 2974 ! psisat - real, saturated soil potential 1 ! 2975 ! bexp - real, soil type "b" parameter 1 ! 2976 ! df1 - real, thermal diffusivity and conductivity 1 ! 2977 ! ice - integer, sea-ice flag (=1: sea-ice, =0: land) 1 ! 2978 ! quartz - real, soil quartz content 1 ! 2979 ! csoil - real, soil heat capacity 1 ! 2980 ! ! 2981 ! input/outputs: ! 2982 ! stc - real, soil temp nsoil ! 2983 ! t1 - real, ground/canopy/snowpack eff skin temp 1 ! 2984 ! tbot - real, bottom soil temp 1 ! 2985 ! sh2o - real, unfrozen soil moisture nsoil ! 2986 ! ! 2987 ! outputs: ! 2988 ! ssoil - real, upward soil heat flux 1 ! 2989 ! ! 2990 ! ==================== end of description ===================== ! 2991 ! 2992 ! --- parameter constants: 2993 real (kind=kind_phys), parameter :: ctfil1 = 0.5 2994 real (kind=kind_phys), parameter :: ctfil2 = 1.0 - ctfil1 2995 2996 ! --- inputs: 2997 integer, intent(in) :: nsoil, ice 2998 2999 real (kind=kind_phys), intent(in) :: smc(nsoil), smcmax, dt, yy, & 3000 & zz1, zsoil(nsoil), zbot, psisat, bexp, df1, quartz, csoil 3001 3002 ! --- input/outputs: 3003 real (kind=kind_phys), intent(inout) :: stc(nsoil), t1, tbot, & 3004 & sh2o(nsoil) 3005 3006 ! --- outputs: 3007 real (kind=kind_phys), intent(out) :: ssoil 3008 3009 ! --- locals: 3010 real (kind=kind_phys) :: ai(nsold), bi(nsold), ci(nsold), oldt1, & 3011 & rhsts(nsold), stcf(nsold), stsoil(nsoil) 3012 3013 integer :: i 3014 3015 ! 3016 !===> ... begin here 3017 ! 3018 oldt1 = t1 3019 do i = 1, nsoil 3020 stsoil(i) = stc(i) 3021 enddo 3022 3023 ! --- ... hrt routine calcs the right hand side of the soil temp dif eqn 3024 3025 if (ice /= 0) then 3026 3027 ! --- ... sea-ice case, glacial-ice case 3028 3029 call hrtice & 3030 ! --- inputs: 3031 & ( nsoil, stc, zsoil, yy, zz1, df1, ice, & 3032 ! --- input/outputs: 3033 & tbot, & 3034 ! --- outputs: 3035 & rhsts, ai, bi, ci & 3036 & ) 3037 3038 call hstep & 3039 ! --- inputs: 3040 & ( nsoil, stc, dt, & 3041 ! --- input/outputs: 3042 & rhsts, ai, bi, ci, & 3043 ! --- outputs: 3044 & stcf & 3045 & ) 3046 3047 else 3048 3049 ! --- ... land-mass case 3050 3051 call hrt & 3052 ! --- inputs: 3053 & ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, & 3054 & zbot, psisat, dt, bexp, df1, quartz, csoil, & 3055 ! --- input/outputs: 3056 & sh2o, & 3057 ! --- outputs: 3058 & rhsts, ai, bi, ci & 3059 & ) 3060 3061 call hstep & 3062 ! --- inputs: 3063 & ( nsoil, stc, dt, & 3064 ! --- input/outputs: 3065 & rhsts, ai, bi, ci, & 3066 ! --- outputs: 3067 & stcf & 3068 & ) 3069 3070 endif 3071 3072 do i = 1, nsoil 3073 stc(i) = stcf(i) 3074 enddo 3075 3076 ! --- ... in the no snowpack case (via routine nopac branch,) update the grnd 3077 ! (skin) temperature here in response to the updated soil temperature 3078 ! profile above. (note: inspection of routine snopac shows that t1 3079 ! below is a dummy variable only, as skin temperature is updated 3080 ! differently in routine snopac) 3081 3082 t1 = (yy + (zz1 - 1.0)*stc(1)) / zz1 3083 t1 = ctfil1*t1 + ctfil2*oldt1 3084 3085 do i = 1, nsoil 3086 stc(i) = ctfil1*stc(i) + ctfil2*stsoil(i) 3087 enddo 3088 3089 ! --- ... calculate surface soil heat flux 3090 3091 ssoil = df1*(stc(1) - t1) / (0.5*zsoil(1)) 3092 3093 ! 3094 return 3095 !................................... 3096 end subroutine shflx 3097 !----------------------------------- 3098 3099 3100 3101 !----------------------------------- 3102 subroutine smflx & 3103 !................................... 3104 ! --- inputs: 3105 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, & 3106 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, & 3107 & edir1, ec1, et1, & 3108 ! --- input/outputs: 3109 & cmc, sh2o, & 3110 ! --- outputs: 3111 & smc, runoff1, runoff2, runoff3, drip & 3112 & ) 3113 3114 ! ===================================================================== ! 3115 ! description: ! 3116 ! ! 3117 ! subroutine smflx calculates soil moisture flux. the soil moisture ! 3118 ! content (smc - a per unit volume measurement) is a dependent variable! 3119 ! that is updated with prognostic eqns. the canopy moisture content ! 3120 ! (cmc) is also updated. frozen ground version: new states added: sh2o! 3121 ! and frozen ground correction factor, frzx and parameter slope. ! 3122 ! ! 3123 ! ! 3124 ! subprogram called: srt, sstep ! 3125 ! ! 3126 ! ==================== defination of variables ==================== ! 3127 ! ! 3128 ! inputs: size ! 3129 ! nsoil - integer, number of soil layers 1 ! 3130 ! dt - real, time step 1 ! 3131 ! kdt - real, 1 ! 3132 ! smcmax - real, porosity 1 ! 3133 ! smcwlt - real, wilting point 1 ! 3134 ! cmcmax - real, maximum canopy water parameters 1 ! 3135 ! prcp1 - real, effective precip 1 ! 3136 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 3137 ! slope - real, linear reservoir coefficient 1 ! 3138 ! frzx - real, frozen ground parameter 1 ! 3139 ! bexp - real, soil type "b" parameter 1 ! 3140 ! dksat - real, saturated soil hydraulic conductivity 1 ! 3141 ! dwsat - real, saturated soil diffusivity 1 ! 3142 ! shdfac - real, aeral coverage of green veg 1 ! 3143 ! edir1 - real, direct soil evaporation 1 ! 3144 ! ec1 - real, canopy water evaporation 1 ! 3145 ! et1 - real, plant transpiration nsoil ! 3146 ! ! 3147 ! input/outputs: ! 3148 ! cmc - real, canopy moisture content 1 ! 3149 ! sh2o - real, unfrozen soil moisture nsoil ! 3150 ! ! 3151 ! outputs: ! 3152 ! smc - real, total soil moisture nsoil ! 3153 ! runoff1 - real, surface runoff not infiltrating sfc 1 ! 3154 ! runoff2 - real, sub surface runoff (baseflow) 1 ! 3155 ! runoff3 - real, excess of porosity 1 ! 3156 ! drip - real, through-fall of precip and/or dew 1 ! 3157 ! ! 3158 ! ==================== end of description ===================== ! 3159 ! 3160 ! --- inputs: 3161 integer, intent(in) :: nsoil 3162 3163 real (kind=kind_phys), intent(in) :: dt, kdt, smcmax, smcwlt, & 3164 & cmcmax, prcp1, slope, frzx, bexp, dksat, dwsat, shdfac, & 3165 & edir1, ec1, et1(nsoil), zsoil(nsoil) 3166 3167 ! --- input/outputs: 3168 real (kind=kind_phys), intent(inout) :: cmc, sh2o(nsoil) 3169 3170 ! --- outputs: 3171 real (kind=kind_phys), intent(out) :: smc(nsoil), runoff1, & 3172 & runoff2, runoff3, drip 3173 3174 ! --- locals: 3175 real (kind=kind_phys) :: dummy, excess, pcpdrp, rhsct, trhsct, & 3176 & rhstt(nsold), sice(nsold), sh2oa(nsold), sh2ofg(nsold), & 3177 & ai(nsold), bi(nsold), ci(nsold) 3178 3179 integer :: i, k 3180 ! 3181 !===> ... begin here 3182 ! 3183 ! --- ... executable code begins here. 3184 3185 dummy = 0.0 3186 3187 ! --- ... compute the right hand side of the canopy eqn term ( rhsct ) 3188 3189 rhsct = shdfac*prcp1 - ec1 3190 3191 ! --- ... convert rhsct (a rate) to trhsct (an amount) and add it to 3192 ! existing cmc. if resulting amt exceeds max capacity, it becomes 3193 ! drip and will fall to the grnd. 3194 3195 drip = 0.0 3196 trhsct = dt * rhsct 3197 excess = cmc + trhsct 3198 3199 if (excess > cmcmax) drip = excess - cmcmax 3200 3201 ! --- ... pcpdrp is the combined prcp1 and drip (from cmc) that goes into 3202 ! the soil 3203 3204 pcpdrp = (1.0 - shdfac)*prcp1 + drip/dt 3205 3206 ! --- ... store ice content at each soil layer before calling srt & sstep 3207 3208 do i = 1, nsoil 3209 sice(i) = smc(i) - sh2o(i) 3210 enddo 3211 3212 ! --- ... call subroutines srt and sstep to solve the soil moisture 3213 ! tendency equations. 3214 3215 ! --- if the infiltrating precip rate is nontrivial, 3216 ! (we consider nontrivial to be a precip total over the time step 3217 ! exceeding one one-thousandth of the water holding capacity of 3218 ! the first soil layer) 3219 ! then call the srt/sstep subroutine pair twice in the manner of 3220 ! time scheme "f" (implicit state, averaged coefficient) 3221 ! of section 2 of kalnay and kanamitsu (1988, mwr, vol 116, 3222 ! pages 1945-1958)to minimize 2-delta-t oscillations in the 3223 ! soil moisture value of the top soil layer that can arise because 3224 ! of the extreme nonlinear dependence of the soil hydraulic 3225 ! diffusivity coefficient and the hydraulic conductivity on the 3226 ! soil moisture state 3227 ! otherwise call the srt/sstep subroutine pair once in the manner of 3228 ! time scheme "d" (implicit state, explicit coefficient) 3229 ! of section 2 of kalnay and kanamitsu 3230 ! pcpdrp is units of kg/m**2/s or mm/s, zsoil is negative depth in m 3231 3232 ! if ( pcpdrp .gt. 0.0 ) then 3233 if ( (pcpdrp*dt) > (0.001*1000.0*(-zsoil(1))*smcmax) ) then 3234 3235 ! --- ... frozen ground version: 3236 ! smc states replaced by sh2o states in srt subr. sh2o & sice states 3237 ! included in sstep subr. frozen ground correction factor, frzx 3238 ! added. all water balance calculations using unfrozen water 3239 3240 call srt & 3241 ! --- inputs: 3242 & ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat, & 3243 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, & 3244 ! --- outputs: 3245 & rhstt, runoff1, runoff2, ai, bi, ci & 3246 & ) 3247 3248 call sstep & 3249 ! --- inputs: 3250 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, & 3251 ! --- input/outputs: 3252 & dummy, rhstt, ai, bi, ci, & 3253 ! --- outputs: 3254 & sh2ofg, runoff3, smc & 3255 & ) 3256 3257 do k = 1, nsoil 3258 sh2oa(k) = (sh2o(k) + sh2ofg(k)) * 0.5 3259 enddo 3260 3261 call srt & 3262 ! --- inputs: 3263 & ( nsoil, edir1, et1, sh2o, sh2oa, pcpdrp, zsoil, dwsat, & 3264 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, & 3265 ! --- outputs: 3266 & rhstt, runoff1, runoff2, ai, bi, ci & 3267 & ) 3268 3269 call sstep & 3270 ! --- inputs: 3271 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, & 3272 ! --- input/outputs: 3273 & cmc, rhstt, ai, bi, ci, & 3274 ! --- outputs: 3275 & sh2o, runoff3, smc & 3276 & ) 3277 3278 else 3279 3280 call srt & 3281 ! --- inputs: 3282 & ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat, & 3283 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, & 3284 ! --- outputs: 3285 & rhstt, runoff1, runoff2, ai, bi, ci & 3286 & ) 3287 3288 call sstep & 3289 ! --- inputs: 3290 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, & 3291 ! --- input/outputs: 3292 & cmc, rhstt, ai, bi, ci, & 3293 ! --- outputs: 3294 & sh2o, runoff3, smc & 3295 & ) 3296 3297 endif 3298 3299 ! runof = runoff 3300 ! 3301 return 3302 !................................... 3303 end subroutine smflx 3304 !----------------------------------- 3305 3306 3307 !----------------------------------- 3308 subroutine snowpack & 3309 !................................... 3310 ! --- inputs: 3311 & ( esd, dtsec, tsnow, tsoil, & 3312 ! --- input/outputs: 3313 & snowh, sndens & 3314 & ) 3315 3316 ! ===================================================================== ! 3317 ! description: ! 3318 ! ! 3319 ! subroutine snowpack calculates compaction of snowpack under ! 3320 ! conditions of increasing snow density, as obtained from an ! 3321 ! approximate solution of e. anderson's differential equation (3.29),! 3322 ! noaa technical report nws 19, by victor koren, 03/25/95. ! 3323 ! subroutine will return new values of snowh and sndens ! 3324 ! ! 3325 ! ! 3326 ! subprogram called: none ! 3327 ! ! 3328 ! ! 3329 ! ==================== defination of variables ==================== ! 3330 ! ! 3331 ! inputs: size ! 3332 ! esd - real, water equivalent of snow (m) 1 ! 3333 ! dtsec - real, time step (sec) 1 ! 3334 ! tsnow - real, snow surface temperature (k) 1 ! 3335 ! tsoil - real, soil surface temperature (k) 1 ! 3336 ! ! 3337 ! input/outputs: ! 3338 ! snowh - real, snow depth (m) 1 ! 3339 ! sndens - real, snow density 1 ! 3340 ! (g/cm3=dimensionless fraction of h2o density) ! 3341 ! ! 3342 ! ==================== end of description ===================== ! 3343 ! 3344 ! --- parameter constants: 3345 real (kind=kind_phys), parameter :: c1 = 0.01 3346 real (kind=kind_phys), parameter :: c2 = 21.0 3347 real (kind=kind_phys), parameter :: kn = 4000.0 3348 3349 ! --- inputs: 3350 real (kind=kind_phys), intent(in) :: esd, dtsec, tsnow, tsoil 3351 3352 ! --- input/outputs: 3353 real (kind=kind_phys), intent(inout) :: snowh, sndens 3354 3355 ! --- locals: 3356 real (kind=kind_phys) :: bfac, dsx, dthr, dw, snowhc, pexp, & 3357 & tavgc, tsnowc, tsoilc, esdc, esdcx 3358 3359 integer :: ipol, j 3360 ! 3361 !===> ... begin here 3362 ! 3363 ! --- ... conversion into simulation units 3364 3365 snowhc = snowh * 100.0 3366 esdc = esd * 100.0 3367 dthr = dtsec / 3600.0 3368 tsnowc = tsnow - tfreez 3369 tsoilc = tsoil - tfreez 3370 3371 ! --- ... calculating of average temperature of snow pack 3372 3373 tavgc = 0.5 * (tsnowc + tsoilc) 3374 3375 ! --- ... calculating of snow depth and density as a result of compaction 3376 ! sndens=ds0*(exp(bfac*esd)-1.)/(bfac*esd) 3377 ! bfac=dthr*c1*exp(0.08*tavgc-c2*ds0) 3378 ! note: bfac*esd in sndens eqn above has to be carefully treated 3379 ! numerically below: 3380 ! c1 is the fractional increase in density (1/(cm*hr)) 3381 ! c2 is a constant (cm3/g) kojima estimated as 21 cms/g 3382 3383 if (esdc > 1.e-2) then 3384 esdcx = esdc 3385 else 3386 esdcx = 1.e-2 3387 endif 3388 3389 bfac = dthr*c1 * exp(0.08*tavgc - c2*sndens) 3390 3391 ! dsx = sndens * ((dexp(bfac*esdc)-1.0) / (bfac*esdc)) 3392 3393 ! --- ... the function of the form (e**x-1)/x imbedded in above expression 3394 ! for dsx was causing numerical difficulties when the denominator "x" 3395 ! (i.e. bfac*esdc) became zero or approached zero (despite the fact 3396 ! that the analytical function (e**x-1)/x has a well defined limit 3397 ! as "x" approaches zero), hence below we replace the (e**x-1)/x 3398 ! expression with an equivalent, numerically well-behaved 3399 ! polynomial expansion. 3400 3401 ! --- ... number of terms of polynomial expansion, and hence its accuracy, 3402 ! is governed by iteration limit "ipol". 3403 ! ipol greater than 9 only makes a difference on double 3404 ! precision (relative errors given in percent %). 3405 ! ipol=9, for rel.error <~ 1.6 e-6 % (8 significant digits) 3406 ! ipol=8, for rel.error <~ 1.8 e-5 % (7 significant digits) 3407 ! ipol=7, for rel.error <~ 1.8 e-4 % ... 3408 3409 ipol = 4 3410 pexp = 0.0 3411 3412 do j = ipol, 1, -1 3413 ! pexp = (1.0 + pexp)*bfac*esdc /real(j+1) 3414 pexp = (1.0 + pexp)*bfac*esdcx/real(j+1) 3415 enddo 3416 pexp = pexp + 1. 3417 3418 dsx = sndens * pexp 3419 3420 ! --- ... above line ends polynomial substitution 3421 ! end of koren formulation 3422 3423 !! --- ... base formulation (cogley et al., 1990) 3424 ! convert density from g/cm3 to kg/m3 3425 3426 !! dsm = sndens * 1000.0 3427 3428 !! dsx = dsm + dtsec*0.5*dsm*gs2*esd / & 3429 !! & (1.e7*exp(-0.02*dsm + kn/(tavgc+273.16)-14.643)) 3430 3431 !! --- ... convert density from kg/m3 to g/cm3 3432 3433 !! dsx = dsx / 1000.0 3434 3435 !! --- ... end of cogley et al. formulation 3436 3437 ! --- ... set upper/lower limit on snow density 3438 3439 dsx = max( min( dsx, 0.40 ), 0.05 ) 3440 sndens = dsx 3441 3442 ! --- ... update of snow depth and density depending on liquid water 3443 ! during snowmelt. assumed that 13% of liquid water can be 3444 ! stored in snow per day during snowmelt till snow density 0.40. 3445 3446 if (tsnowc >= 0.0) then 3447 dw = 0.13 * dthr / 24.0 3448 sndens = sndens*(1.0 - dw) + dw 3449 if (sndens > 0.40) sndens = 0.40 3450 endif 3451 3452 ! --- ... calculate snow depth (cm) from snow water equivalent and snow 3453 ! density. change snow depth units to meters 3454 3455 snowhc = esdc / sndens 3456 snowh = snowhc * 0.01 3457 3458 ! 3459 return 3460 !................................... 3461 end subroutine snowpack 3462 !----------------------------------- 3463 3464 3465 !*********************************************! 3466 ! section-3 3rd or lower level subprograms ! 3467 !*********************************************! 3468 3469 3470 !----------------------------------- 3471 subroutine devap & 3472 !................................... 3473 ! --- inputs: 3474 & ( etp1, smc, shdfac, smcmax, smcdry, fxexp, & 3475 ! --- outputs: 3476 & edir1 & 3477 & ) 3478 3479 ! ===================================================================== ! 3480 ! description: ! 3481 ! ! 3482 ! subroutine devap calculates direct soil evaporation ! 3483 ! ! 3484 ! ! 3485 ! subprogram called: none ! 3486 ! ! 3487 ! ! 3488 ! ==================== defination of variables ==================== ! 3489 ! ! 3490 ! inputs: size ! 3491 ! etp1 - real, potential evaporation 1 ! 3492 ! smc - real, unfrozen soil moisture 1 ! 3493 ! shdfac - real, aeral coverage of green vegetation 1 ! 3494 ! smcmax - real, porosity (sat val of soil mois) 1 ! 3495 ! smcdry - real, dry soil mois threshold 1 ! 3496 ! fxexp - real, bare soil evaporation exponent 1 ! 3497 ! ! 3498 ! outputs: ! 3499 ! edir1 - real, direct soil evaporation 1 ! 3500 ! ! 3501 ! ==================== end of description ===================== ! 3502 ! 3503 ! --- inputs: 3504 real (kind=kind_phys), intent(in) :: etp1, smc, shdfac, smcmax, & 3505 & smcdry, fxexp 3506 3507 ! --- outputs: 3508 real (kind=kind_phys), intent(out) :: edir1 3509 3510 ! --- locals: 3511 real (kind=kind_phys) :: fx, sratio 3512 ! 3513 !===> ... begin here 3514 ! 3515 ! --- ... direct evap a function of relative soil moisture availability, 3516 ! linear when fxexp=1. 3517 ! fx > 1 represents demand control 3518 ! fx < 1 represents flux control 3519 3520 sratio = (smc - smcdry) / (smcmax - smcdry) 3521 3522 if (sratio > 0.0) then 3523 fx = sratio**fxexp 3524 fx = max ( min ( fx, 1.0 ), 0.0 ) 3525 else 3526 fx = 0.0 3527 endif 3528 3529 ! --- ... allow for the direct-evap-reducing effect of shade 3530 3531 edir1 = fx * ( 1.0 - shdfac ) * etp1 3532 ! 3533 return 3534 !................................... 3535 end subroutine devap 3536 !----------------------------------- 3537 3538 3539 !----------------------------------- 3540 subroutine frh2o & 3541 !................................... 3542 ! --- inputs: 3543 & ( tkelv, smc, sh2o, smcmax, bexp, psis, & 3544 ! --- outputs: 3545 & liqwat & 3546 & ) 3547 3548 ! ===================================================================== ! 3549 ! description: ! 3550 ! ! 3551 ! subroutine frh2o calculates amount of supercooled liquid soil water ! 3552 ! content if temperature is below 273.15k (t0). requires newton-type ! 3553 ! iteration to solve the nonlinear implicit equation given in eqn 17 ! 3554 ! of koren et al (1999, jgr, vol 104(d16), 19569-19585). ! 3555 ! ! 3556 ! new version (june 2001): much faster and more accurate newton ! 3557 ! iteration achieved by first taking log of eqn cited above -- less ! 3558 ! than 4 (typically 1 or 2) iterations achieves convergence. also, ! 3559 ! explicit 1-step solution option for special case of parameter ck=0, ! 3560 ! which reduces the original implicit equation to a simpler explicit ! 3561 ! form, known as the "flerchinger eqn". improved handling of solution ! 3562 ! in the limit of freezing point temperature t0. ! 3563 ! ! 3564 ! subprogram called: none ! 3565 ! ! 3566 ! ! 3567 ! ==================== defination of variables ==================== ! 3568 ! ! 3569 ! inputs: size ! 3570 ! tkelv - real, temperature (k) 1 ! 3571 ! smc - real, total soil moisture content (volumetric) 1 ! 3572 ! sh2o - real, liquid soil moisture content (volumetric) 1 ! 3573 ! smcmax - real, saturation soil moisture content 1 ! 3574 ! bexp - real, soil type "b" parameter 1 ! 3575 ! psis - real, saturated soil matric potential 1 ! 3576 ! ! 3577 ! outputs: ! 3578 ! liqwat - real, supercooled liquid water content 1 ! 3579 ! ! 3580 ! ==================== end of description ===================== ! 3581 ! 3582 ! --- constant parameters: 3583 real (kind=kind_phys), parameter :: ck = 8.0 3584 ! real (kind=kind_phys), parameter :: ck = 0.0 3585 real (kind=kind_phys), parameter :: blim = 5.5 3586 real (kind=kind_phys), parameter :: error = 0.005 3587 3588 ! --- inputs: 3589 real (kind=kind_phys), intent(in) :: tkelv, smc, sh2o, smcmax, & 3590 & bexp, psis 3591 3592 ! --- outputs: 3593 real (kind=kind_phys), intent(out) :: liqwat 3594 3595 ! --- locals: 3596 real (kind=kind_phys) :: bx, denom, df, dswl, fk, swl, swlk 3597 3598 integer :: nlog, kcount 3599 ! 3600 !===> ... begin here 3601 ! 3602 ! --- ... limits on parameter b: b < 5.5 (use parameter blim) 3603 ! simulations showed if b > 5.5 unfrozen water content is 3604 ! non-realistically high at very low temperatures. 3605 3606 bx = bexp 3607 if (bexp > blim) bx = blim 3608 3609 ! --- ... initializing iterations counter and iterative solution flag. 3610 3611 nlog = 0 3612 kcount= 0 3613 3614 ! --- ... if temperature not significantly below freezing (t0), sh2o = smc 3615 3616 if (tkelv > (tfreez-1.e-3)) then 3617 3618 liqwat = smc 3619 3620 else 3621 3622 if (ck /= 0.0) then 3623 3624 ! --- ... option 1: iterated solution for nonzero ck 3625 ! in koren et al, jgr, 1999, eqn 17 3626 3627 ! --- ... initial guess for swl (frozen content) 3628 3629 swl = smc - sh2o 3630 3631 ! --- ... keep within bounds. 3632 3633 swl = max( min( swl, smc-0.02 ), 0.0 ) 3634 3635 ! --- ... start of iterations 3636 3637 do while ( (nlog < 10) .and. (kcount == 0) ) 3638 nlog = nlog + 1 3639 3640 df = alog( (psis*gs2/lsubf) * ( (1.0 + ck*swl)**2.0 ) & 3641 & * (smcmax/(smc-swl))**bx ) - alog(-(tkelv-tfreez)/tkelv) 3642 3643 denom = 2.0*ck/(1.0 + ck*swl) + bx/(smc - swl) 3644 swlk = swl - df/denom 3645 3646 ! --- ... bounds useful for mathematical solution. 3647 3648 swlk = max( min( swlk, smc-0.02 ), 0.0 ) 3649 3650 ! --- ... mathematical solution bounds applied. 3651 3652 dswl = abs(swlk - swl) 3653 swl = swlk 3654 3655 ! --- ... if more than 10 iterations, use explicit method (ck=0 approx.) 3656 ! when dswl less or eq. error, no more iterations required. 3657 3658 if ( dswl <= error ) then 3659 kcount = kcount + 1 3660 endif 3661 enddo ! end do_while_loop 3662 3663 ! --- ... bounds applied within do-block are valid for physical solution. 3664 3665 liqwat = smc - swl 3666 3667 endif ! end if_ck_block 3668 3669 ! --- ... option 2: explicit solution for flerchinger eq. i.e. ck=0 3670 ! in koren et al., jgr, 1999, eqn 17 3671 ! apply physical bounds to flerchinger solution 3672 3673 if (kcount == 0) then 3674 fk = ( ( (lsubf/(gs2*(-psis))) & 3675 & * ((tkelv-tfreez)/tkelv) )**(-1/bx) ) * smcmax 3676 3677 fk = max( fk, 0.02 ) 3678 3679 liqwat = min( fk, smc ) 3680 endif 3681 3682 endif ! end if_tkelv_block 3683 ! 3684 return 3685 !................................... 3686 end subroutine frh2o 3687 !----------------------------------- 3688 3689 3690 !----------------------------------- 3691 subroutine hrt & 3692 !................................... 3693 ! --- inputs: 3694 & ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, & 3695 & zbot, psisat, dt, bexp, df1, quartz, csoil, & 3696 ! --- input/outputs: 3697 & sh2o, & 3698 ! --- outputs: 3699 & rhsts, ai, bi, ci & 3700 & ) 3701 3702 ! ===================================================================== ! 3703 ! description: ! 3704 ! ! 3705 ! subroutine hrt calculates the right hand side of the time tendency ! 3706 ! term of the soil thermal diffusion equation. also to compute ! 3707 ! (prepare) the matrix coefficients for the tri-diagonal matrix of ! 3708 ! the implicit time scheme. ! 3709 ! ! 3710 ! subprogram called: tbnd, snksrc, tmpavg ! 3711 ! ! 3712 ! ! 3713 ! ==================== defination of variables ==================== ! 3714 ! ! 3715 ! inputs: size ! 3716 ! nsoil - integer, number of soil layers 1 ! 3717 ! stc - real, soil temperature nsoil ! 3718 ! smc - real, total soil moisture nsoil ! 3719 ! smcmax - real, porosity 1 ! 3720 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 3721 ! yy - real, 1 ! 3722 ! zz1 - real, soil temperture at the top soil column 1 ! 3723 ! tbot - real, bottom soil temp 1 ! 3724 ! zbot - real, specify depth of lower bd soil 1 ! 3725 ! psisat - real, saturated soil potential 1 ! 3726 ! dt - real, time step 1 ! 3727 ! bexp - real, soil type "b" parameter 1 ! 3728 ! df1 - real, thermal diffusivity 1 ! 3729 ! quartz - real, soil quartz content 1 ! 3730 ! csoil - real, soil heat capacity 1 ! 3731 ! ! 3732 ! input/outputs: ! 3733 ! sh2o - real, unfrozen soil moisture nsoil ! 3734 ! ! 3735 ! outputs: ! 3736 ! rhsts - real, time tendency of soil thermal diffusion nsoil ! 3737 ! ai - real, matrix coefficients nsold ! 3738 ! bi - real, matrix coefficients nsold ! 3739 ! ci - real, matrix coefficients nsold ! 3740 ! ! 3741 ! ==================== end of description ===================== ! 3742 ! 3743 ! --- inputs: 3744 integer, intent(in) :: nsoil 3745 3746 real (kind=kind_phys), intent(in) :: stc(nsoil), smc(nsoil), & 3747 & smcmax, zsoil(nsoil), yy, zz1, tbot, zbot, psisat, dt, & 3748 & bexp, df1, quartz, csoil 3749 3750 ! --- input/outputs: 3751 real (kind=kind_phys), intent(inout) :: sh2o(nsoil) 3752 3753 ! --- outputs: 3754 real (kind=kind_phys), intent(out) :: rhsts(nsoil), ai(nsold), & 3755 & bi(nsold), ci(nsold) 3756 3757 ! --- locals: 3758 real (kind=kind_phys) :: ddz, ddz2, denom, df1n, df1k, dtsdz, & 3759 & dtsdz2, hcpct, qtot, ssoil, sice, tavg, tbk, tbk1, & 3760 & tsnsr, tsurf 3761 3762 integer :: i, k 3763 3764 logical :: itavg 3765 3766 ! 3767 !===> ... begin here 3768 ! 3769 ! --- ... initialize logical for soil layer temperature averaging. 3770 3771 itavg = .true. 3772 ! itavg = .false. 3773 3774 ! === begin section for top soil layer 3775 3776 ! --- ... calc the heat capacity of the top soil layer 3777 3778 hcpct = sh2o(1)*cph2o2 + (1.0 - smcmax)*csoil & 3779 & + (smcmax - smc(1))*cp2 + (smc(1) - sh2o(1))*cpice1 3780 3781 ! --- ... calc the matrix coefficients ai, bi, and ci for the top layer 3782 3783 ddz = 1.0 / ( -0.5*zsoil(2) ) 3784 ai(1) = 0.0 3785 ci(1) = (df1*ddz) / ( zsoil(1)*hcpct ) 3786 bi(1) = -ci(1) + df1 / ( 0.5*zsoil(1)*zsoil(1)*hcpct*zz1 ) 3787 3788 ! --- ... calculate the vertical soil temp gradient btwn the 1st and 2nd soil 3789 ! layers. then calculate the subsurface heat flux. use the temp 3790 ! gradient and subsfc heat flux to calc "right-hand side tendency 3791 ! terms", or "rhsts", for top soil layer. 3792 3793 dtsdz = (stc(1) - stc(2)) / (-0.5*zsoil(2)) 3794 ssoil = df1 * (stc(1) - yy) / (0.5*zsoil(1)*zz1) 3795 rhsts(1) = (df1*dtsdz - ssoil) / (zsoil(1)*hcpct) 3796 3797 ! --- ... next capture the vertical difference of the heat flux at top and 3798 ! bottom of first soil layer for use in heat flux constraint applied to 3799 ! potential soil freezing/thawing in routine snksrc. 3800 3801 qtot = ssoil - df1*dtsdz 3802 3803 ! --- ... if temperature averaging invoked (itavg=true; else skip): 3804 ! set temp "tsurf" at top of soil column (for use in freezing soil 3805 ! physics later in subroutine snksrc). if snowpack content is 3806 ! zero, then tsurf expression below gives tsurf = skin temp. if 3807 ! snowpack is nonzero (hence argument zz1=1), then tsurf expression 3808 ! below yields soil column top temperature under snowpack. then 3809 ! calculate temperature at bottom interface of 1st soil layer for use 3810 ! later in subroutine snksrc 3811 3812 if (itavg) then 3813 3814 tsurf = (yy + (zz1-1)*stc(1)) / zz1 3815 3816 call tbnd & 3817 ! --- inputs: 3818 & ( stc(1), stc(2), zsoil, zbot, 1, nsoil, & 3819 ! --- outputs: 3820 & tbk & 3821 & ) 3822 3823 endif 3824 3825 ! --- ... calculate frozen water content in 1st soil layer. 3826 3827 sice = smc(1) - sh2o(1) 3828 3829 ! --- ... if frozen water present or any of layer-1 mid-point or bounding 3830 ! interface temperatures below freezing, then call snksrc to 3831 ! compute heat source/sink (and change in frozen water content) 3832 ! due to possible soil water phase change 3833 3834 if ( (sice > 0.0) .or. (tsurf < tfreez) .or. & 3835 & (stc(1) < tfreez) .or. (tbk < tfreez) ) then 3836 3837 if (itavg) then 3838 3839 call tmpavg & 3840 ! --- inputs: 3841 & ( tsurf, stc(1), tbk, zsoil, nsoil, 1, & 3842 ! --- outputs: 3843 & tavg & 3844 & ) 3845 3846 else 3847 3848 tavg = stc(1) 3849 3850 endif ! end if_itavg_block 3851 3852 call snksrc & 3853 ! --- inputs: 3854 & ( nsoil, 1, tavg, smc(1), smcmax, psisat, bexp, dt, & 3855 & qtot, zsoil, & 3856 ! --- input/outputs: 3857 & sh2o(1), & 3858 ! --- outputs: 3859 & tsnsr & 3860 & ) 3861 3862 3863 rhsts(1) = rhsts(1) - tsnsr / ( zsoil(1)*hcpct ) 3864 3865 endif ! end if_sice_block 3866 3867 ! === this ends section for top soil layer. 3868 3869 ! --- ... initialize ddz2 3870 3871 ddz2 = 0.0 3872 3873 ! --- ... loop thru the remaining soil layers, repeating the above process 3874 ! (except subsfc or "ground" heat flux not repeated in lower layers) 3875 3876 df1k = df1 3877 3878 do k = 2, nsoil 3879 3880 ! --- ... calculate heat capacity for this soil layer. 3881 3882 hcpct = sh2o(k)*cph2o2 + (1.0 - smcmax)*csoil & 3883 & + (smcmax - smc(k))*cp2 + (smc(k) - sh2o(k))*cpice1 3884 3885 if (k /= nsoil) then 3886 3887 ! --- ... this section for layer 2 or greater, but not last layer. 3888 ! calculate thermal diffusivity for this layer. 3889 3890 call tdfcnd & 3891 ! --- inputs: 3892 & ( smc(k), quartz, smcmax, sh2o(k), & 3893 ! --- outputs: 3894 & df1n & 3895 & ) 3896 3897 ! --- ... calc the vertical soil temp gradient thru this layer 3898 3899 denom = 0.5 * (zsoil(k-1) - zsoil(k+1)) 3900 dtsdz2 = (stc(k) - stc(k+1)) / denom 3901 3902 ! --- ... calc the matrix coef, ci, after calc'ng its partial product 3903 3904 ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1)) 3905 ci(k) = -df1n*ddz2 / ((zsoil(k-1) - zsoil(k)) * hcpct) 3906 3907 ! --- ... if temperature averaging invoked (itavg=true; else skip): 3908 ! calculate temp at bottom of layer. 3909 3910 if (itavg) then 3911 3912 call tbnd & 3913 ! --- inputs: 3914 & ( stc(k), stc(k+1), zsoil, zbot, k, nsoil, & 3915 ! --- outputs: 3916 & tbk1 & 3917 & ) 3918 3919 endif 3920 3921 else 3922 3923 ! --- ... special case of bottom soil layer: calculate thermal diffusivity 3924 ! for bottom layer. 3925 3926 call tdfcnd & 3927 ! --- inputs: 3928 & ( smc(k), quartz, smcmax, sh2o(k), & 3929 ! --- outputs: 3930 & df1n & 3931 & ) 3932 3933 ! --- ... calc the vertical soil temp gradient thru bottom layer. 3934 3935 denom = 0.5 * (zsoil(k-1) + zsoil(k)) - zbot 3936 dtsdz2 = (stc(k) - tbot) / denom 3937 3938 ! --- ... set matrix coef, ci to zero if bottom layer. 3939 3940 ci(k) = 0.0 3941 3942 ! --- ... if temperature averaging invoked (itavg=true; else skip): 3943 ! calculate temp at bottom of last layer. 3944 3945 if (itavg) then 3946 3947 call tbnd & 3948 ! --- inputs: 3949 & ( stc(k), tbot, zsoil, zbot, k, nsoil, & 3950 ! --- outputs: 3951 & tbk1 & 3952 & ) 3953 3954 endif 3955 3956 endif ! end if_k_block 3957 3958 ! --- ... calculate rhsts for this layer after calc'ng a partial product. 3959 3960 denom = (zsoil(k) - zsoil(k-1)) * hcpct 3961 rhsts(k) = ( df1n*dtsdz2 - df1k*dtsdz ) / denom 3962 3963 qtot = -1.0 * denom * rhsts(k) 3964 sice = smc(k) - sh2o(k) 3965 3966 if ( (sice > 0.0) .or. (tbk < tfreez) .or. & 3967 & (stc(k) < tfreez) .or. (tbk1 < tfreez) ) then 3968 3969 if (itavg) then 3970 3971 call tmpavg & 3972 ! --- inputs: 3973 & ( tbk, stc(k), tbk1, zsoil, nsoil, k, & 3974 ! --- outputs: 3975 & tavg & 3976 & ) 3977 3978 else 3979 tavg = stc(k) 3980 endif 3981 3982 call snksrc & 3983 ! --- inputs: 3984 & ( nsoil, k, tavg, smc(k), smcmax, psisat, bexp, dt, & 3985 & qtot, zsoil, & 3986 ! --- input/outputs: 3987 & sh2o(k), & 3988 ! --- outputs: 3989 & tsnsr & 3990 & ) 3991 3992 rhsts(k) = rhsts(k) - tsnsr/denom 3993 endif 3994 3995 ! --- ... calc matrix coefs, ai, and bi for this layer. 3996 3997 ai(k) = - df1 * ddz / ((zsoil(k-1) - zsoil(k)) * hcpct) 3998 bi(k) = -(ai(k) + ci(k)) 3999 4000 ! --- ... reset values of df1, dtsdz, ddz, and tbk for loop to next soil layer. 4001 4002 tbk = tbk1 4003 df1k = df1n 4004 dtsdz = dtsdz2 4005 ddz = ddz2 4006 4007 enddo ! end do_k_loop 4008 4009 ! 4010 return 4011 !................................... 4012 end subroutine hrt 4013 !----------------------------------- 4014 4015 4016 !----------------------------------- 4017 subroutine hrtice & 4018 !................................... 4019 ! --- inputs: 4020 & ( nsoil, stc, zsoil, yy, zz1, df1, ice, & 4021 ! --- input/outputs: 4022 & tbot, & 4023 ! --- outputs: 4024 & rhsts, ai, bi, ci & 4025 & ) 4026 4027 ! ===================================================================== ! 4028 ! description: ! 4029 ! ! 4030 ! subroutine hrtice calculates the right hand side of the time tendency! 4031 ! term of the soil thermal diffusion equation for sea-ice (ice = 1) or ! 4032 ! glacial-ice (ice). compute (prepare) the matrix coefficients for the ! 4033 ! tri-diagonal matrix of the implicit time scheme. ! 4034 ! (note: this subroutine only called for sea-ice or glacial ice, but ! 4035 ! not for non-glacial land (ice = 0). ! 4036 ! ! 4037 ! subprogram called: none ! 4038 ! ! 4039 ! ! 4040 ! ==================== defination of variables ==================== ! 4041 ! ! 4042 ! inputs: size ! 4043 ! nsoil - integer, number of soil layers 1 ! 4044 ! stc - real, soil temperature nsoil ! 4045 ! zsoil - real, soil depth (negative sign, as below grd) nsoil ! 4046 ! yy - real, soil temperature at the top of column 1 ! 4047 ! zz1 - real, 1 ! 4048 ! df1 - real, thermal diffusivity and conductivity 1 ! 4049 ! ice - integer, sea-ice flag (=1: sea-ice, =0: land) 1 ! 4050 ! ! 4051 ! input/outputs: ! 4052 ! tbot - real, bottom soil temperature 1 ! 4053 ! ! 4054 ! outputs: ! 4055 ! rhsts - real, time tendency of soil thermal diffusion nsoil ! 4056 ! ai - real, matrix coefficients nsold ! 4057 ! bi - real, matrix coefficients nsold ! 4058 ! ci - real, matrix coefficients nsold ! 4059 ! ! 4060 ! ==================== end of description ===================== ! 4061 ! 4062 ! --- inputs: 4063 integer, intent(in) :: nsoil, ice 4064 4065 real (kind=kind_phys), intent(in) :: stc(nsoil), zsoil(nsoil), & 4066 & yy, zz1, df1 4067 4068 ! --- input/outputs: 4069 real (kind=kind_phys), intent(inout) :: tbot 4070 4071 ! --- outputs: 4072 real (kind=kind_phys), intent(out) :: rhsts(nsoil), ai(nsold), & 4073 & bi(nsold), ci(nsold) 4074 4075 ! --- locals: 4076 real (kind=kind_phys) :: ddz, ddz2, denom, dtsdz, dtsdz2, & 4077 & hcpct, ssoil, zbot 4078 4079 integer :: k 4080 4081 ! 4082 !===> ... begin here 4083 ! 4084 ! --- ... set a nominal universal value of the sea-ice specific heat capacity, 4085 ! hcpct = 1880.0*917.0 = 1.72396e+6 (source: fei chen, 1995) 4086 ! set bottom of sea-ice pack temperature: tbot = 271.16 4087 ! set a nominal universal value of glacial-ice specific heat capacity, 4088 ! hcpct = 2100.0*900.0 = 1.89000e+6 (source: bob grumbine, 2005) 4089 ! tbot passed in as argument, value from global data set 4090 4091 if (ice == 1) then 4092 ! --- ... sea-ice 4093 hcpct = 1.72396e+6 4094 tbot = 271.16 4095 else 4096 ! --- ... glacial-ice 4097 hcpct = 1.89000e+6 4098 endif 4099 4100 ! --- ... the input argument df1 is a universally constant value of sea-ice 4101 ! and glacial-ice thermal diffusivity, set in sflx as df1 = 2.2. 4102 4103 ! --- ... set ice pack depth. use tbot as ice pack lower boundary temperature 4104 ! (that of unfrozen sea water at bottom of sea ice pack). assume ice 4105 ! pack is of n=nsoil layers spanning a uniform constant ice pack 4106 ! thickness as defined by zsoil(nsoil) in routine sflx. 4107 ! if glacial-ice, set zbot = -25 meters 4108 4109 if (ice == 1) then 4110 ! --- ... sea-ice 4111 zbot = zsoil(nsoil) 4112 else 4113 ! --- ... glacial-ice 4114 zbot = -25.0 4115 endif 4116 4117 ! --- ... calc the matrix coefficients ai, bi, and ci for the top layer 4118 4119 ddz = 1.0 / (-0.5*zsoil(2)) 4120 ai(1) = 0.0 4121 ci(1) = (df1*ddz) / (zsoil(1)*hcpct) 4122 bi(1) = -ci(1) + df1 / (0.5*zsoil(1)*zsoil(1)*hcpct*zz1) 4123 4124 ! --- ... calc the vertical soil temp gradient btwn the top and 2nd soil 4125 ! layers. recalc/adjust the soil heat flux. use the gradient and 4126 ! flux to calc rhsts for the top soil layer. 4127 4128 dtsdz = (stc(1) - stc(2)) / (-0.5*zsoil(2)) 4129 ssoil = df1 * (stc(1) - yy) / (0.5*zsoil(1)*zz1) 4130 rhsts(1) = (df1*dtsdz - ssoil) / (zsoil(1)*hcpct) 4131 4132 ! --- ... initialize ddz2 4133 4134 ddz2 = 0.0 4135 4136 ! --- ... loop thru the remaining soil layers, repeating the above process 4137 4138 do k = 2, nsoil 4139 4140 if (k /= nsoil) then 4141 4142 ! --- ... calc the vertical soil temp gradient thru this layer. 4143 4144 denom = 0.5 * (zsoil(k-1) - zsoil(k+1)) 4145 dtsdz2 = (stc(k) - stc(k+1)) / denom 4146 4147 ! --- ... calc the matrix coef, ci, after calc'ng its partial product. 4148 4149 ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1)) 4150 ci(k) = -df1*ddz2 / ((zsoil(k-1) - zsoil(k))*hcpct) 4151 4152 else 4153 4154 ! --- ... calc the vertical soil temp gradient thru the lowest layer. 4155 4156 dtsdz2 = (stc(k) - tbot) & 4157 & / (0.5*(zsoil(k-1) + zsoil(k)) - zbot) 4158 4159 ! --- ... set matrix coef, ci to zero. 4160 4161 ci(k) = 0.0 4162 4163 endif ! end if_k_block 4164 4165 ! --- ... calc rhsts for this layer after calc'ng a partial product. 4166 4167 denom = (zsoil(k) - zsoil(k-1)) * hcpct 4168 rhsts(k) = (df1*dtsdz2 - df1*dtsdz) / denom 4169 4170 ! --- ... calc matrix coefs, ai, and bi for this layer. 4171 4172 ai(k) = - df1*ddz / ((zsoil(k-1) - zsoil(k)) * hcpct) 4173 bi(k) = -(ai(k) + ci(k)) 4174 4175 ! --- ... reset values of dtsdz and ddz for loop to next soil lyr. 4176 4177 dtsdz = dtsdz2 4178 ddz = ddz2 4179 4180 enddo ! end do_k_loop 4181 ! 4182 return 4183 !................................... 4184 end subroutine hrtice 4185 !----------------------------------- 4186 4187 4188 !----------------------------------- 4189 subroutine hstep & 4190 !................................... 4191 ! --- inputs: 4192 & ( nsoil, stcin, dt, & 4193 ! --- input/outputs: 4194 & rhsts, ai, bi, ci, & 4195 ! --- outputs: 4196 & stcout & 4197 & ) 4198 4199 ! ===================================================================== ! 4200 ! description: ! 4201 ! ! 4202 ! subroutine hstep calculates/updates the soil temperature field. ! 4203 ! ! 4204 ! subprogram called: rosr12 ! 4205 ! ! 4206 ! ! 4207 ! ==================== defination of variables ==================== ! 4208 ! ! 4209 ! inputs: size ! 4210 ! nsoil - integer, number of soil layers 1 ! 4211 ! stcin - real, soil temperature nsoil ! 4212 ! dt - real, time step 1 ! 4213 ! ! 4214 ! input/outputs: ! 4215 ! rhsts - real, time tendency of soil thermal diffusion nsoil ! 4216 ! ai - real, matrix coefficients nsold ! 4217 ! bi - real, matrix coefficients nsold ! 4218 ! ci - real, matrix coefficients nsold ! 4219 ! ! 4220 ! outputs: ! 4221 ! stcout - real, updated soil temperature nsoil ! 4222 ! ! 4223 ! ==================== end of description ===================== ! 4224 ! 4225 ! --- inputs: 4226 integer, intent(in) :: nsoil 4227 4228 real (kind=kind_phys), intent(in) :: stcin(nsoil), dt 4229 4230 ! --- input/outputs: 4231 real (kind=kind_phys), intent(inout) :: rhsts(nsoil), & 4232 & ai(nsold), bi(nsold), ci(nsold) 4233 4234 ! --- outputs: 4235 real (kind=kind_phys), intent(out) :: stcout(nsoil) 4236 4237 ! --- locals: 4238 integer :: k 4239 4240 real (kind=kind_phys) :: ciin(nsold), rhstsin(nsoil) 4241 4242 ! 4243 !===> ... begin here 4244 ! 4245 ! --- ... create finite difference values for use in rosr12 routine 4246 4247 do k = 1, nsoil 4248 rhsts(k) = rhsts(k) * dt 4249 ai(k) = ai(k) * dt 4250 bi(k) = 1.0 + bi(k)*dt 4251 ci(k) = ci(k) * dt 4252 enddo 4253 4254 ! --- ... copy values for input variables before call to rosr12 4255 4256 do k = 1, nsoil 4257 rhstsin(k) = rhsts(k) 4258 enddo 4259 4260 do k = 1, nsold 4261 ciin(k) = ci(k) 4262 enddo 4263 4264 ! --- ... solve the tri-diagonal matrix equation 4265 4266 call rosr12 & 4267 ! --- inputs: 4268 & ( nsoil, ai, bi, rhstsin, & 4269 ! --- input/outputs: 4270 & ciin, & 4271 ! --- outputs: 4272 & ci, rhsts & 4273 & ) 4274 4275 ! --- ... calc/update the soil temps using matrix solution 4276 4277 do k = 1, nsoil 4278 stcout(k) = stcin(k) + ci(k) 4279 enddo 4280 ! 4281 return 4282 !................................... 4283 end subroutine hstep 4284 !----------------------------------- 4285 4286 4287 !----------------------------------- 4288 subroutine rosr12 & 4289 !................................... 4290 ! --- inputs: 4291 & ( nsoil, a, b, d, & 4292 ! --- input/outputs: 4293 & c, & 4294 ! --- outputs: 4295 & p, delta & 4296 & ) 4297 4298 ! ===================================================================== ! 4299 ! description: ! 4300 ! ! 4301 ! subroutine rosr12 inverts (solve) the tri-diagonal matrix problem ! 4302 ! shown below: ! 4303 ! ! 4304 ! ### ### ### ### ### ###! 4305 ! #b(1), c(1), 0 , 0 , 0 , . . . , 0 # # # # #! 4306 ! #a(2), b(2), c(2), 0 , 0 , . . . , 0 # # # # #! 4307 ! # 0 , a(3), b(3), c(3), 0 , . . . , 0 # # # # d(3) #! 4308 ! # 0 , 0 , a(4), b(4), c(4), . . . , 0 # # p(4) # # d(4) #! 4309 ! # 0 , 0 , 0 , a(5), b(5), . . . , 0 # # p(5) # # d(5) #! 4310 ! # . . # # . # = # . #! 4311 ! # . . # # . # # . #! 4312 ! # . . # # . # # . #! 4313 ! # 0 , . . . , 0 , a(m-2), b(m-2), c(m-2), 0 # #p(m-2)# #d(m-2)#! 4314 ! # 0 , . . . , 0 , 0 , a(m-1), b(m-1), c(m-1)# #p(m-1)# #d(m-1)#! 4315 ! # 0 , . . . , 0 , 0 , 0 , a(m) , b(m) # # p(m) # # d(m) #! 4316 ! ### ### ### ### ### ###! 4317 ! ! 4318 ! subprogram called: none ! 4319 ! ! 4320 ! ! 4321 ! ==================== defination of variables ==================== ! 4322 ! ! 4323 ! inputs: size ! 4324 ! nsoil - integer, number of soil layers 1 ! 4325 ! a - real, matrix coefficients nsoil ! 4326 ! b - real, matrix coefficients nsoil ! 4327 ! d - real, soil water time tendency nsoil ! 4328 ! ! 4329 ! input/outputs: ! 4330 ! c - real, matrix coefficients nsoil ! 4331 ! ! 4332 ! outputs: ! 4333 ! p - real, nsoil ! 4334 ! delta - real, nsoil ! 4335 ! ! 4336 ! ==================== end of description ===================== ! 4337 ! 4338 ! --- inputs: 4339 integer, intent(in) :: nsoil 4340 4341 real (kind=kind_phys), dimension(nsoil), intent(in) :: a, b, d 4342 4343 ! --- input/outputs: 4344 real (kind=kind_phys), dimension(nsoil), intent(inout) :: c 4345 4346 ! --- outputs: 4347 real (kind=kind_phys), dimension(nsoil), intent(out) :: p, delta 4348 4349 ! --- locals: 4350 integer :: k, kk 4351 4352 ! 4353 !===> ... begin here 4354 ! 4355 ! --- ... initialize eqn coef c for the lowest soil layer 4356 4357 c(nsoil) = 0.0 4358 4359 ! --- ... solve the coefs for the 1st soil layer 4360 4361 p(1) = -c(1) / b(1) 4362 delta(1) = d(1) / b(1) 4363 4364 ! --- ... solve the coefs for soil layers 2 thru nsoil 4365 4366 do k = 2, nsoil 4367 p(k) = -c(k) * ( 1.0 / (b(k) + a (k)*p(k-1)) ) 4368 delta(k) = (d(k) - a(k)*delta(k-1)) & 4369 & * ( 1.0 / (b(k) + a(k)*p(k-1)) ) 4370 enddo 4371 4372 ! --- ... set p to delta for lowest soil layer 4373 4374 p(nsoil) = delta(nsoil) 4375 4376 ! --- ... adjust p for soil layers 2 thru nsoil 4377 4378 do k = 2, nsoil 4379 kk = nsoil - k + 1 4380 p(kk) = p(kk)*p(kk+1) + delta(kk) 4381 enddo 4382 ! 4383 return 4384 !................................... 4385 end subroutine rosr12 4386 !----------------------------------- 4387 4388 4389 !----------------------------------- 4390 subroutine snksrc & 4391 !................................... 4392 ! --- inputs: 4393 & ( nsoil, k, tavg, smc, smcmax, psisat, bexp, dt, & 4394 & qtot, zsoil, & 4395 ! --- input/outputs: 4396 & sh2o, & 4397 ! --- outputs: 4398 & tsrc & 4399 & ) 4400 4401 ! ===================================================================== ! 4402 ! description: ! 4403 ! ! 4404 ! subroutine snksrc calculates sink/source term of the termal ! 4405 ! diffusion equation. ! 4406 ! ! 4407 ! subprograms called: frh2o ! 4408 ! ! 4409 ! ! 4410 ! ==================== defination of variables ==================== ! 4411 ! ! 4412 ! inputs: size ! 4413 ! nsoil - integer, number of soil layers 1 ! 4414 ! k - integer, index of soil layers 1 ! 4415 ! tavg - real, soil layer average temperature 1 ! 4416 ! smc - real, total soil moisture 1 ! 4417 ! smcmax - real, porosity 1 ! 4418 ! psisat - real, saturated soil potential 1 ! 4419 ! bexp - real, soil type "b" parameter 1 ! 4420 ! dt - real, time step 1 ! 4421 ! qtot - real, tot vertical diff of heat flux 1 ! 4422 ! zsoil - real, soil layer depth below ground (negative) nsoil ! 4423 ! ! 4424 ! input/outputs: ! 4425 ! sh2o - real, available liqued water 1 ! 4426 ! ! 4427 ! outputs: ! 4428 ! tsrc - real, heat source/sink 1 ! 4429 ! ! 4430 ! ==================== end of description ===================== ! 4431 ! 4432 ! --- parameter constants: 4433 real (kind=kind_phys), parameter :: dh2o = 1.0000e3 4434 4435 ! --- inputs: 4436 integer, intent(in) :: nsoil, k 4437 4438 real (kind=kind_phys), intent(in) :: tavg, smc, smcmax, psisat, & 4439 & bexp, dt, qtot, zsoil(nsoil) 4440 4441 ! --- input/outputs: 4442 real (kind=kind_phys), intent(inout) :: sh2o 4443 4444 ! --- outputs: 4445 real (kind=kind_phys), intent(out) :: tsrc 4446 4447 ! --- locals: 4448 real (kind=kind_phys) :: dz, free, xh2o 4449 4450 ! --- external functions: 4451 ! real (kind=kind_phys) :: frh2o 4452 4453 ! 4454 !===> ... begin here 4455 ! 4456 if (k == 1) then 4457 dz = -zsoil(1) 4458 else 4459 dz = zsoil(k-1) - zsoil(k) 4460 endif 4461 4462 ! --- ... via function frh2o, compute potential or 'equilibrium' unfrozen 4463 ! supercooled free water for given soil type and soil layer temperature. 4464 ! function frh20 invokes eqn (17) from v. koren et al (1999, jgr, vol. 4465 ! 104, pg 19573). (aside: latter eqn in journal in centigrade units. 4466 ! routine frh2o use form of eqn in kelvin units.) 4467 4468 ! free = frh2o( tavg,smc,sh2o,smcmax,bexp,psisat ) 4469 4470 call frh2o & 4471 ! --- inputs: 4472 & ( tavg, smc, sh2o, smcmax, bexp, psisat, & 4473 ! --- outputs: 4474 & free & 4475 & ) 4476 4477 4478 ! --- ... in next block of code, invoke eqn 18 of v. koren et al (1999, jgr, 4479 ! vol. 104, pg 19573.) that is, first estimate the new amountof liquid 4480 ! water, 'xh2o', implied by the sum of (1) the liquid water at the begin 4481 ! of current time step, and (2) the freeze of thaw change in liquid 4482 ! water implied by the heat flux 'qtot' passed in from routine hrt. 4483 ! second, determine if xh2o needs to be bounded by 'free' (equil amt) or 4484 ! if 'free' needs to be bounded by xh2o. 4485 4486 xh2o = sh2o + qtot*dt / (dh2o*lsubf*dz) 4487 4488 ! --- ... first, if freezing and remaining liquid less than lower bound, then 4489 ! reduce extent of freezing, thereby letting some or all of heat flux 4490 ! qtot cool the soil temp later in routine hrt. 4491 4492 if ( xh2o < sh2o .and. xh2o < free) then 4493 if ( free > sh2o ) then 4494 xh2o = sh2o 4495 else 4496 xh2o = free 4497 endif 4498 endif 4499 4500 ! --- ... second, if thawing and the increase in liquid water greater than 4501 ! upper bound, then reduce extent of thaw, thereby letting some or 4502 ! all of heat flux qtot warm the soil temp later in routine hrt. 4503 4504 if ( xh2o > sh2o .and. xh2o > free ) then 4505 if ( free < sh2o ) then 4506 xh2o = sh2o 4507 else 4508 xh2o = free 4509 endif 4510 endif 4511 4512 xh2o = max( min( xh2o, smc ), 0.0 ) 4513 4514 ! --- ... calculate phase-change heat source/sink term for use in routine hrt 4515 ! and update liquid water to reflcet final freeze/thaw increment. 4516 4517 tsrc = -dh2o * lsubf * dz * (xh2o - sh2o) / dt 4518 sh2o = xh2o 4519 ! 4520 return 4521 !................................... 4522 end subroutine snksrc 4523 !----------------------------------- 4524 4525 4526 !----------------------------------- 4527 subroutine srt & 4528 !................................... 4529 ! --- inputs: 4530 & ( nsoil, edir, et, sh2o, sh2oa, pcpdrp, zsoil, dwsat, & 4531 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, & 4532 ! --- outputs: 4533 & rhstt, runoff1, runoff2, ai, bi, ci & 4534 & ) 4535 4536 ! ===================================================================== ! 4537 ! description: ! 4538 ! subroutine srt calculates the right hand side of the time tendency ! 4539 ! term of the soil water diffusion equation. also to compute ! 4540 ! ( prepare ) the matrix coefficients for the tri-diagonal matrix ! 4541 ! of the implicit time scheme. ! 4542 ! ! 4543 ! subprogram called: wdfcnd ! 4544 ! ! 4545 ! ! 4546 ! ==================== defination of variables ==================== ! 4547 ! ! 4548 ! inputs: size ! 4549 ! nsoil - integer, number of soil layers 1 ! 4550 ! edir - real, direct soil evaporation 1 ! 4551 ! et - real, plant transpiration nsoil ! 4552 ! sh2o - real, unfrozen soil moisture nsoil ! 4553 ! sh2oa - real, nsoil ! 4554 ! pcpdrp - real, combined prcp and drip 1 ! 4555 ! zsoil - real, soil layer depth below ground nsoil ! 4556 ! dwsat - real, saturated soil diffusivity 1 ! 4557 ! dksat - real, saturated soil hydraulic conductivity 1 ! 4558 ! smcmax - real, porosity 1 ! 4559 ! bexp - real, soil type "b" parameter 1 ! 4560 ! dt - real, time step 1 ! 4561 ! smcwlt - real, wilting point 1 ! 4562 ! slope - real, linear reservoir coefficient 1 ! 4563 ! kdt - real, 1 ! 4564 ! frzx - real, frozen ground parameter 1 ! 4565 ! sice - real, ice content at each soil layer nsoil ! 4566 ! ! 4567 ! outputs: ! 4568 ! rhstt - real, soil water time tendency nsoil ! 4569 ! runoff1 - real, surface runoff not infiltrating sfc 1 ! 4570 ! runoff2 - real, sub surface runoff (baseflow) 1 ! 4571 ! ai - real, matrix coefficients nsold ! 4572 ! bi - real, matrix coefficients nsold ! 4573 ! ci - real, matrix coefficients nsold ! 4574 ! ! 4575 ! ==================== end of description ===================== ! 4576 ! 4577 ! --- inputs: 4578 integer, intent(in) :: nsoil 4579 4580 real (kind=kind_phys), dimension(nsoil), intent(in) :: et, & 4581 & sh2o, sh2oa, zsoil, sice 4582 4583 real (kind=kind_phys), intent(in) :: edir, pcpdrp, dwsat, dksat, & 4584 & smcmax, smcwlt, bexp, dt, slope, kdt, frzx 4585 4586 ! --- outputs: 4587 real (kind=kind_phys), intent(out) :: runoff1, runoff2, & 4588 & rhstt(nsoil), ai(nsold), bi(nsold), ci(nsold) 4589 4590 4591 ! --- locals: 4592 real (kind=kind_phys) :: acrt, dd, ddt, ddz, ddz2, denom, denom2, & 4593 & dice, dsmdz, dsmdz2, dt1, fcr, infmax, mxsmc, mxsmc2, px, & 4594 & numer, pddum, sicemax, slopx, smcav, sstt, sum, val, wcnd, & 4595 & wcnd2, wdf, wdf2, dmax(nsold) 4596 4597 integer :: cvfrz, ialp1, iohinf, j, jj, k, ks 4598 ! 4599 !===> ... begin here 4600 ! 4601 ! --- ... frozen ground version: 4602 ! reference frozen ground parameter, cvfrz, is a shape parameter 4603 ! of areal distribution function of soil ice content which equals 4604 ! 1/cv. cv is a coefficient of spatial variation of soil ice content. 4605 ! based on field data cv depends on areal mean of frozen depth, and 4606 ! it close to constant = 0.6 if areal mean frozen depth is above 20 cm. 4607 ! that is why parameter cvfrz = 3 (int{1/0.6*0.6}). current logic 4608 ! doesn't allow cvfrz be bigger than 3 4609 4610 parameter (cvfrz = 3) 4611 4612 c ---------------------------------------------------------------------- 4613 ! --- ... determine rainfall infiltration rate and runoff. include 4614 ! the infiltration formule from schaake and koren model. 4615 ! modified by q duan 4616 4617 iohinf = 1 4618 4619 ! --- ... let sicemax be the greatest, if any, frozen water content within 4620 ! soil layers. 4621 4622 sicemax = 0.0 4623 do ks = 1, nsoil 4624 if (sice(ks) > sicemax) sicemax = sice(ks) 4625 enddo 4626 4627 ! --- ... determine rainfall infiltration rate and runoff 4628 4629 pddum = pcpdrp 4630 runoff1 = 0.0 4631 4632 if (pcpdrp /= 0.0) then 4633 4634 ! --- ... modified by q. duan, 5/16/94 4635 4636 dt1 = dt/86400. 4637 smcav = smcmax - smcwlt 4638 dmax(1) = -zsoil(1) * smcav 4639 4640 ! --- ... frozen ground version: 4641 4642 dice = -zsoil(1) * sice(1) 4643 4644 dmax(1) = dmax(1)*(1.0 - (sh2oa(1)+sice(1)-smcwlt)/smcav) 4645 dd = dmax(1) 4646 4647 do ks = 2, nsoil 4648 4649 ! --- ... frozen ground version: 4650 4651 dice = dice + ( zsoil(ks-1) - zsoil(ks) ) * sice(ks) 4652 4653 dmax(ks) = (zsoil(ks-1)-zsoil(ks))*smcav 4654 dmax(ks) = dmax(ks)*(1.0 - (sh2oa(ks)+sice(ks)-smcwlt)/smcav) 4655 dd = dd + dmax(ks) 4656 enddo 4657 4658 ! --- ... val = (1.-exp(-kdt*sqrt(dt1))) 4659 ! in below, remove the sqrt in above 4660 4661 val = 1.0 - exp(-kdt*dt1) 4662 ddt = dd * val 4663 4664 px = pcpdrp * dt 4665 if (px < 0.0) px = 0.0 4666 4667 infmax = (px*(ddt/(px+ddt)))/dt 4668 4669 ! --- ... frozen ground version: 4670 ! reduction of infiltration based on frozen ground parameters 4671 4672 fcr = 1. 4673 4674 if (dice > 1.e-2) then 4675 acrt = cvfrz * frzx / dice 4676 sum = 1. 4677 4678 ialp1 = cvfrz - 1 4679 do j = 1, ialp1 4680 k = 1 4681 4682 do jj = j+1,ialp1 4683 k = k * jj 4684 enddo 4685 4686 sum = sum + (acrt**( cvfrz-j)) / float (k) 4687 enddo 4688 4689 fcr = 1.0 - exp(-acrt) * sum 4690 endif 4691 4692 infmax = infmax * fcr 4693 4694 ! --- ... correction of infiltration limitation: 4695 ! if infmax .le. hydrolic conductivity assign infmax the value 4696 ! of hydrolic conductivity 4697 4698 ! mxsmc = max ( sh2oa(1), sh2oa(2) ) 4699 mxsmc = sh2oa(1) 4700 4701 call wdfcnd & 4702 ! --- inputs: 4703 & ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax, & 4704 ! --- outputs: 4705 & wdf, wcnd & 4706 & ) 4707 4708 infmax = max( infmax, wcnd ) 4709 infmax = min( infmax, px ) 4710 4711 if (pcpdrp > infmax) then 4712 runoff1 = pcpdrp - infmax 4713 pddum = infmax 4714 endif 4715 4716 endif ! end if_pcpdrp_block 4717 4718 ! --- ... to avoid spurious drainage behavior, 'upstream differencing' 4719 ! in line below replaced with new approach in 2nd line: 4720 ! 'mxsmc = max(sh2oa(1), sh2oa(2))' 4721 4722 mxsmc = sh2oa(1) 4723 4724 call wdfcnd & 4725 ! --- inputs: 4726 & ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax, & 4727 ! --- outputs: 4728 & wdf, wcnd & 4729 & ) 4730 4731 ! --- ... calc the matrix coefficients ai, bi, and ci for the top layer 4732 4733 ddz = 1.0 / ( -.5*zsoil(2) ) 4734 ai(1) = 0.0 4735 bi(1) = wdf * ddz / ( -zsoil(1) ) 4736 ci(1) = -bi(1) 4737 4738 ! --- ... calc rhstt for the top layer after calc'ng the vertical soil 4739 ! moisture gradient btwn the top and next to top layers. 4740 4741 dsmdz = ( sh2o(1) - sh2o(2) ) / ( -.5*zsoil(2) ) 4742 rhstt(1) = (wdf*dsmdz + wcnd - pddum + edir + et(1)) / zsoil(1) 4743 sstt = wdf * dsmdz + wcnd + edir + et(1) 4744 4745 ! --- ... initialize ddz2 4746 4747 ddz2 = 0.0 4748 4749 ! --- ... loop thru the remaining soil layers, repeating the abv process 4750 4751 do k = 2, nsoil 4752 denom2 = (zsoil(k-1) - zsoil(k)) 4753 4754 if (k /= nsoil) then 4755 slopx = 1.0 4756 4757 ! --- ... again, to avoid spurious drainage behavior, 'upstream differencing' 4758 ! in line below replaced with new approach in 2nd line: 4759 ! 'mxsmc2 = max (sh2oa(k), sh2oa(k+1))' 4760 4761 mxsmc2 = sh2oa(k) 4762 4763 call wdfcnd & 4764 ! --- inputs: 4765 & ( mxsmc2, smcmax, bexp, dksat, dwsat, sicemax, & 4766 ! --- outputs: 4767 & wdf2, wcnd2 & 4768 & ) 4769 4770 ! --- ... calc some partial products for later use in calc'ng rhstt 4771 4772 denom = (zsoil(k-1) - zsoil(k+1)) 4773 dsmdz2 = (sh2o(k) - sh2o(k+1)) / (denom * 0.5) 4774 4775 ! --- ... calc the matrix coef, ci, after calc'ng its partial product 4776 4777 ddz2 = 2.0 / denom 4778 ci(k) = -wdf2 * ddz2 / denom2 4779 4780 else ! if_k_block 4781 4782 ! --- ... slope of bottom layer is introduced 4783 4784 slopx = slope 4785 4786 ! --- ... retrieve the soil water diffusivity and hydraulic conductivity 4787 ! for this layer 4788 4789 call wdfcnd & 4790 ! --- inputs: 4791 & ( sh2oa(nsoil), smcmax, bexp, dksat, dwsat, sicemax, & 4792 ! --- outputs: 4793 & wdf2, wcnd2 & 4794 & ) 4795 4796 ! --- ... calc a partial product for later use in calc'ng rhstt 4797 dsmdz2 = 0.0 4798 4799 ! --- ... set matrix coef ci to zero 4800 4801 ci(k) = 0.0 4802 4803 endif ! end if_k_block 4804 4805 ! --- ... calc rhstt for this layer after calc'ng its numerator 4806 4807 numer = wdf2*dsmdz2 + slopx*wcnd2 - wdf*dsmdz - wcnd + et(k) 4808 rhstt(k) = numer / (-denom2) 4809 4810 ! --- ... calc matrix coefs, ai, and bi for this layer 4811 4812 ai(k) = -wdf * ddz / denom2 4813 bi(k) = -( ai(k) + ci(k) ) 4814 4815 ! --- ... reset values of wdf, wcnd, dsmdz, and ddz for loop to next lyr 4816 ! runoff2: sub-surface or baseflow runoff 4817 4818 if (k == nsoil) then 4819 runoff2 = slopx * wcnd2 4820 endif 4821 4822 if (k /= nsoil) then 4823 wdf = wdf2 4824 wcnd = wcnd2 4825 dsmdz= dsmdz2 4826 ddz = ddz2 4827 endif 4828 enddo ! end do_k_loop 4829 ! 4830 return 4831 !................................... 4832 end subroutine srt 4833 !----------------------------------- 4834 4835 4836 !----------------------------------- 4837 subroutine sstep & 4838 !................................... 4839 ! --- inputs: 4840 & ( nsoil, sh2oin, rhsct, dt, smcmax, cmcmax, zsoil, sice, & 4841 ! --- input/outputs: 4842 & cmc, rhstt, ai, bi, ci, & 4843 ! --- outputs: 4844 & sh2oout, runoff3, smc & 4845 & ) 4846 4847 ! ===================================================================== ! 4848 ! description: ! 4849 ! subroutine sstep calculates/updates soil moisture content values ! 4850 ! and canopy moisture content values. ! 4851 ! ! 4852 ! subprogram called: rosr12 ! 4853 ! ! 4854 ! ! 4855 ! ==================== defination of variables ==================== ! 4856 ! ! 4857 ! inputs: size ! 4858 ! nsoil - integer, number of soil layers 1 ! 4859 ! sh2oin - real, unfrozen soil moisture nsoil ! 4860 ! rhsct - real, 1 ! 4861 ! dt - real, time step 1 ! 4862 ! smcmax - real, porosity 1 ! 4863 ! cmcmax - real, maximum canopy water parameters 1 ! 4864 ! zsoil - real, soil layer depth below ground nsoil ! 4865 ! sice - real, ice content at each soil layer nsoil ! 4866 ! ! 4867 ! input/outputs: ! 4868 ! cmc - real, canopy moisture content 1 ! 4869 ! rhstt - real, soil water time tendency nsoil ! 4870 ! ai - real, matrix coefficients nsold ! 4871 ! bi - real, matrix coefficients nsold ! 4872 ! ci - real, matrix coefficients nsold ! 4873 ! ! 4874 ! outputs: ! 4875 ! sh2oout - real, updated soil moisture content nsoil ! 4876 ! runoff3 - real, excess of porosity 1 ! 4877 ! smc - real, total soil moisture nsoil ! 4878 ! ! 4879 ! ==================== end of description ===================== ! 4880 ! 4881 ! --- input: 4882 integer, intent(in) :: nsoil 4883 4884 real (kind=kind_phys), dimension(nsoil), intent(in) :: sh2oin, & 4885 & zsoil, sice 4886 4887 real (kind=kind_phys), intent(in) :: rhsct, dt, smcmax, cmcmax 4888 4889 ! --- inout/outputs: 4890 real (kind=kind_phys), intent(inout) :: cmc, rhstt(nsoil), & 4891 & ai(nsold), bi(nsold), ci(nsold) 4892 4893 ! --- outputs: 4894 real (kind=kind_phys), intent(out) :: sh2oout(nsoil), runoff3, & 4895 & smc(nsoil) 4896 4897 ! --- locals: 4898 real (kind=kind_phys) :: ciin(nsold), rhsttin(nsoil), ddz, stot, & 4899 & wplus 4900 4901 integer :: i, k, kk11 4902 ! 4903 !===> ... begin here 4904 ! 4905 ! --- ... create 'amount' values of variables to be input to the 4906 ! tri-diagonal matrix routine. 4907 4908 do k = 1, nsoil 4909 rhstt(k) = rhstt(k) * dt 4910 ai(k) = ai(k) * dt 4911 bi(k) = 1. + bi(k) * dt 4912 ci(k) = ci(k) * dt 4913 enddo 4914 4915 ! --- ... copy values for input variables before call to rosr12 4916 4917 do k = 1, nsoil 4918 rhsttin(k) = rhstt(k) 4919 enddo 4920 4921 do k = 1, nsold 4922 ciin(k) = ci(k) 4923 enddo 4924 4925 ! --- ... call rosr12 to solve the tri-diagonal matrix 4926 4927 call rosr12 & 4928 ! --- inputs: 4929 & ( nsoil, ai, bi, rhsttin, & 4930 ! --- input/outputs: 4931 & ciin, & 4932 ! --- outputs: 4933 & ci, rhstt & 4934 & ) 4935 4936 ! --- ... sum the previous smc value and the matrix solution to get 4937 ! a new value. min allowable value of smc will be 0.02. 4938 ! runoff3: runoff within soil layers 4939 4940 wplus = 0.0 4941 runoff3 = 0.0 4942 ddz = -zsoil(1) 4943 4944 do k = 1, nsoil 4945 if (k /= 1) ddz = zsoil(k - 1) - zsoil(k) 4946 4947 sh2oout(k) = sh2oin(k) + ci(k) + wplus/ddz 4948 4949 stot = sh2oout(k) + sice(k) 4950 if (stot > smcmax) then 4951 if (k == 1) then 4952 ddz = -zsoil(1) 4953 else 4954 kk11 = k - 1 4955 ddz = -zsoil(k) + zsoil(kk11) 4956 endif 4957 4958 wplus = (stot - smcmax) * ddz 4959 else 4960 wplus = 0.0 4961 endif 4962 4963 smc(k) = max( min( stot, smcmax ), 0.02 ) 4964 sh2oout(k) = max( smc(k)-sice(k), 0.0 ) 4965 enddo 4966 4967 runoff3 = wplus 4968 4969 ! --- ... update canopy water content/interception (cmc). convert rhsct to 4970 ! an 'amount' value and add to previous cmc value to get new cmc. 4971 4972 cmc = cmc + dt*rhsct 4973 if (cmc < 1.e-20) cmc = 0.0 4974 cmc = min( cmc, cmcmax ) 4975 ! 4976 return 4977 !................................... 4978 end subroutine sstep 4979 !----------------------------------- 4980 4981 4982 !----------------------------------- 4983 subroutine tbnd & 4984 !................................... 4985 ! --- inputs: 4986 & ( tu, tb, zsoil, zbot, k, nsoil, & 4987 ! --- outputs: 4988 & tbnd1 & 4989 & ) 4990 4991 ! ===================================================================== ! 4992 ! description: ! 4993 ! ! 4994 ! subroutine tbnd calculates temperature on the boundary of the ! 4995 ! layer by interpolation of the middle layer temperatures ! 4996 ! ! 4997 ! subprogram called: none ! 4998 ! ! 4999 ! ==================== defination of variables ==================== ! 5000 ! ! 5001 ! inputs: size ! 5002 ! tu - real, soil temperature 1 ! 5003 ! tb - real, bottom soil temp 1 ! 5004 ! zsoil - real, soil layer depth nsoil ! 5005 ! zbot - real, specify depth of lower bd soil 1 ! 5006 ! k - integer, soil layer index 1 ! 5007 ! nsoil - integer, number of soil layers 1 ! 5008 ! ! 5009 ! outputs: ! 5010 ! tbnd1 - real, temperature at bottom interface of the lyr 1 ! 5011 ! ! 5012 ! ==================== end of description ===================== ! 5013 ! 5014 ! --- input: 5015 integer, intent(in) :: k, nsoil 5016 5017 real (kind=kind_phys), intent(in) :: tu, tb, zbot, zsoil(nsoil) 5018 5019 ! --- output: 5020 real (kind=kind_phys), intent(out) :: tbnd1 5021 5022 ! --- locals: 5023 real (kind=kind_phys) :: zb, zup 5024 5025 ! --- ... use surface temperature on the top of the first layer 5026 5027 if (k == 1) then 5028 zup = 0.0 5029 else 5030 zup = zsoil(k-1) 5031 endif 5032 5033 ! --- ... use depth of the constant bottom temperature when interpolate 5034 ! temperature into the last layer boundary 5035 5036 if (k == nsoil) then 5037 zb = 2.0*zbot - zsoil(k) 5038 else 5039 zb = zsoil(k+1) 5040 endif 5041 5042 ! --- ... linear interpolation between the average layer temperatures 5043 5044 tbnd1 = tu + (tb-tu)*(zup-zsoil(k))/(zup-zb) 5045 ! 5046 return 5047 !................................... 5048 end subroutine tbnd 5049 !----------------------------------- 5050 5051 5052 !----------------------------------- 5053 subroutine tmpavg & 5054 !................................... 5055 ! --- inputs: 5056 & ( tup, tm, tdn, zsoil, nsoil, k, & 5057 ! --- outputs: 5058 & tavg & 5059 & ) 5060 5061 ! ===================================================================== ! 5062 ! description: ! 5063 ! subroutine tmpavg calculates soil layer average temperature (tavg) ! 5064 ! in freezing/thawing layer using up, down, and middle layer ! 5065 ! temperatures (tup, tdn, tm), where tup is at top boundary of ! 5066 ! layer, tdn is at bottom boundary of layer. tm is layer prognostic ! 5067 ! state temperature. ! 5068 ! ! 5069 ! ! 5070 ! subprogram called: none ! 5071 ! ! 5072 ! 5073 ! ==================== defination of variables ==================== ! 5074 ! ! 5075 ! inputs: size ! 5076 ! tup - real, temperature ar top boundary of layer 1 ! 5077 ! tm - real, layer prognostic state temperature 1 ! 5078 ! tdn - real, temperature ar bottom boundary of layer 1 ! 5079 ! zsoil - real, soil layer depth nsoil ! 5080 ! nsoil - integer, number of soil layers 1 ! 5081 ! k - integer, layer index 1 ! 5082 ! outputs: ! 5083 ! tavg - real, soil layer average temperature 1 ! 5084 ! ! 5085 ! ==================== end of description ===================== ! 5086 ! 5087 ! --- input: 5088 integer, intent(in) :: nsoil, k 5089 5090 real (kind=kind_phys), intent(in) :: tup, tm, tdn, zsoil(nsoil) 5091 5092 ! --- output: 5093 real (kind=kind_phys), intent(out) :: tavg 5094 5095 ! --- locals: 5096 real (kind=kind_phys) :: dz, dzh, x0, xdn, xup 5097 ! 5098 !===> ... begin here 5099 ! 5100 if (k == 1) then 5101 dz = -zsoil(1) 5102 else 5103 dz = zsoil(k-1) - zsoil(k) 5104 endif 5105 5106 dzh = dz * 0.5 5107 5108 if (tup < tfreez) then 5109 5110 if (tm < tfreez) then 5111 if (tdn < tfreez) then ! tup, tm, tdn < t0 5112 tavg = (tup + 2.0*tm + tdn) / 4.0 5113 else ! tup & tm < t0, tdn >= t0 5114 x0 = (tfreez - tm) * dzh / (tdn - tm) 5115 tavg = 0.5*(tup*dzh + tm*(dzh+x0)+tfreez*(2.*dzh-x0)) / dz 5116 endif 5117 else 5118 if (tdn < tfreez) then ! tup < t0, tm >= t0, tdn < t0 5119 xup = (tfreez-tup) * dzh / (tm-tup) 5120 xdn = dzh - (tfreez-tm) * dzh / (tdn-tm) 5121 tavg = 0.5*(tup*xup + tfreez*(2.*dz-xup-xdn)+tdn*xdn) / dz 5122 else ! tup < t0, tm >= t0, tdn >= t0 5123 xup = (tfreez-tup) * dzh / (tm-tup) 5124 tavg = 0.5*(tup*xup + tfreez*(2.*dz-xup)) / dz 5125 endif 5126 endif 5127 5128 else ! if_tup_block 5129 5130 if (tm < tfreez) then 5131 if (tdn < tfreez) then ! tup >= t0, tm < t0, tdn < t0 5132 xup = dzh - (tfreez-tup) * dzh / (tm-tup) 5133 tavg = 0.5*(tfreez*(dz-xup) + tm*(dzh+xup)+tdn*dzh) / dz 5134 else ! tup >= t0, tm < t0, tdn >= t0 5135 xup = dzh - (tfreez-tup) * dzh / (tm-tup) 5136 xdn = (tfreez-tm) * dzh / (tdn-tm) 5137 tavg = 0.5 * (tfreez*(2.*dz-xup-xdn) + tm*(xup+xdn)) / dz 5138 endif 5139 else 5140 if (tdn < tfreez) then ! tup >= t0, tm >= t0, tdn < t0 5141 xdn = dzh - (tfreez-tm) * dzh / (tdn-tm) 5142 tavg = (tfreez*(dz-xdn) + 0.5*(tfreez+tdn)*xdn) / dz 5143 else ! tup >= t0, tm >= t0, tdn >= t0 5144 tavg = (tup + 2.0*tm + tdn) / 4.0 5145 endif 5146 endif 5147 5148 endif ! end if_tup_block 5149 ! 5150 return 5151 !................................... 5152 end subroutine tmpavg 5153 !----------------------------------- 5154 5155 5156 !----------------------------------- 5157 subroutine transp & 5158 !................................... 5159 ! --- inputs: 5160 & ( nsoil, nroot, etp1, smc, smcwlt, smcref, & 5161 & cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis, & 5162 ! --- outputs: 5163 & et1 & 5164 & ) 5165 5166 ! ===================================================================== ! 5167 ! description: ! 5168 ! subroutine transp calculates transpiration for the veg class. ! 5169 ! ! 5170 ! subprogram called: none ! 5171 ! ! 5172 ! ! 5173 ! ==================== defination of variables ==================== ! 5174 ! ! 5175 ! inputs: size ! 5176 ! nsoil - integer, number of soil layers 1 ! 5177 ! nroot - integer, number of root layers 1 ! 5178 ! etp1 - real, potential evaporation 1 ! 5179 ! smc - real, unfrozen soil moisture nsoil ! 5180 ! smcwlt - real, wilting point 1 ! 5181 ! smcref - real, soil mois threshold 1 ! 5182 ! cmc - real, canopy moisture content 1 ! 5183 ! cmcmax - real, maximum canopy water parameters 1 ! 5184 ! zsoil - real, soil layer depth below ground nsoil ! 5185 ! shdfac - real, aeral coverage of green vegetation 1 ! 5186 ! pc - real, plant coeff 1 ! 5187 ! cfactr - real, canopy water parameters 1 ! 5188 ! rtdis - real, root distribution nsoil ! 5189 ! ! 5190 ! outputs: ! 5191 ! et1 - real, plant transpiration nsoil ! 5192 ! ! 5193 ! ==================== end of description ===================== ! 5194 ! 5195 ! --- input: 5196 integer, intent(in) :: nsoil, nroot 5197 5198 real (kind=kind_phys), intent(in) :: etp1, smcwlt, smcref, & 5199 & cmc, cmcmax, shdfac, pc, cfactr 5200 5201 real (kind=kind_phys), dimension(nsoil), intent(in) :: smc, & 5202 & zsoil, rtdis 5203 5204 ! --- output: 5205 real (kind=kind_phys), dimension(nsoil), intent(out) :: et1 5206 5207 ! --- locals: 5208 real (kind=kind_phys) :: denom, etp1a, rtx, sgx, gx(7) 5209 5210 integer :: i, k 5211 ! 5212 !===> ... begin here 5213 ! 5214 ! --- ... initialize plant transp to zero for all soil layers. 5215 5216 do k = 1, nsoil 5217 et1(k) = 0.0 5218 enddo 5219 5220 ! --- ... calculate an 'adjusted' potential transpiration 5221 ! if statement below to avoid tangent linear problems near zero 5222 ! note: gx and other terms below redistribute transpiration by layer, 5223 ! et(k), as a function of soil moisture availability, while preserving 5224 ! total etp1a. 5225 5226 if (cmc /= 0.0) then 5227 etp1a = shdfac * pc * etp1 * (1.0 - (cmc /cmcmax) ** cfactr) 5228 else 5229 etp1a = shdfac * pc * etp1 5230 endif 5231 5232 sgx = 0.0 5233 do i = 1, nroot 5234 gx(i) = ( smc(i) - smcwlt ) / ( smcref - smcwlt ) 5235 gx(i) = max ( min ( gx(i), 1.0 ), 0.0 ) 5236 sgx = sgx + gx(i) 5237 enddo 5238 sgx = sgx / nroot 5239 5240 denom = 0.0 5241 do i = 1, nroot 5242 rtx = rtdis(i) + gx(i) - sgx 5243 gx(i) = gx(i) * max ( rtx, 0.0 ) 5244 denom = denom + gx(i) 5245 enddo 5246 if (denom <= 0.0) denom = 1.0 5247 5248 do i = 1, nroot 5249 et1(i) = etp1a * gx(i) / denom 5250 enddo 5251 5252 ! --- ... above code assumes a vertically uniform root distribution 5253 ! code below tests a variable root distribution 5254 5255 ! et(1) = ( zsoil(1) / zsoil(nroot) ) * gx * etp1a 5256 ! et(1) = ( zsoil(1) / zsoil(nroot) ) * etp1a 5257 5258 ! --- ... using root distribution as weighting factor 5259 5260 ! et(1) = rtdis(1) * etp1a 5261 ! et(1) = etp1a * part(1) 5262 5263 ! --- ... loop down thru the soil layers repeating the operation above, 5264 ! but using the thickness of the soil layer (rather than the 5265 ! absolute depth of each layer) in the final calculation. 5266 5267 ! do k = 2, nroot 5268 ! gx = ( smc(k) - smcwlt ) / ( smcref - smcwlt ) 5269 ! gx = max ( min ( gx, 1.0 ), 0.0 ) 5270 ! --- ... test canopy resistance 5271 ! gx = 1.0 5272 ! et(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot))*gx*etp1a 5273 ! et(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot))*etp1a 5274 5275 ! --- ... using root distribution as weighting factor 5276 5277 ! et(k) = rtdis(k) * etp1a 5278 ! et(k) = etp1a*part(k) 5279 ! enddo 5280 5281 ! 5282 return 5283 !................................... 5284 end subroutine transp 5285 !----------------------------------- 5286 5287 5288 !----------------------------------- 5289 subroutine wdfcnd & 5290 !................................... 5291 ! --- inputs: 5292 & ( smc, smcmax, bexp, dksat, dwsat, sicemax, & 5293 ! --- outputs: 5294 & wdf, wcnd & 5295 & ) 5296 5297 ! ===================================================================== ! 5298 ! description: ! 5299 ! subroutine wdfcnd calculates soil water diffusivity and soil ! 5300 ! hydraulic conductivity. ! 5301 ! ! 5302 ! subprogram called: none ! 5303 ! ! 5304 ! ! 5305 ! ==================== defination of variables ==================== ! 5306 ! ! 5307 ! inputs: size ! 5308 ! smc - real, layer total soil moisture 1 ! 5309 ! smcmax - real, porosity 1 ! 5310 ! bexp - real, soil type "b" parameter 1 ! 5311 ! dksat - real, saturated soil hydraulic conductivity 1 ! 5312 ! dwsat - real, saturated soil diffusivity 1 ! 5313 ! sicemax - real, max frozen water content in soil layer 1 ! 5314 ! ! 5315 ! outputs: ! 5316 ! wdf - real, soil water diffusivity 1 ! 5317 ! wcnd - real, soil hydraulic conductivity 1 ! 5318 ! ! 5319 ! ==================== end of description ===================== ! 5320 ! 5321 ! --- input: 5322 real (kind=kind_phys), intent(in) :: smc, smcmax, bexp, dksat, & 5323 & dwsat, sicemax 5324 5325 ! --- output: 5326 real (kind=kind_phys), intent(out) :: wdf, wcnd 5327 5328 ! --- locals: 5329 real (kind=kind_phys) :: expon, factr1, factr2, vkwgt 5330 ! 5331 !===> ... begin here 5332 ! 5333 ! --- ... calc the ratio of the actual to the max psbl soil h2o content 5334 5335 factr1 = 0.2 / smcmax 5336 factr2 = smc / smcmax 5337 5338 ! --- ... prep an expntl coef and calc the soil water diffusivity 5339 5340 expon = bexp + 2.0 5341 wdf = dwsat * factr2 ** expon 5342 5343 ! --- ... frozen soil hydraulic diffusivity. very sensitive to the vertical 5344 ! gradient of unfrozen water. the latter gradient can become very 5345 ! extreme in freezing/thawing situations, and given the relatively 5346 ! few and thick soil layers, this gradient sufferes serious 5347 ! trunction errors yielding erroneously high vertical transports of 5348 ! unfrozen water in both directions from huge hydraulic diffusivity. 5349 ! therefore, we found we had to arbitrarily constrain wdf 5350 ! 5351 ! version d_10cm: ....... factr1 = 0.2/smcmax 5352 ! weighted approach....... pablo grunmann, 28_sep_1999. 5353 5354 if (sicemax > 0.0) then 5355 vkwgt = 1.0 / (1.0 + (500.0*sicemax)**3.0) 5356 wdf = vkwgt*wdf + (1.0- vkwgt)*dwsat*factr1**expon 5357 endif 5358 5359 ! --- ... reset the expntl coef and calc the hydraulic conductivity 5360 5361 expon = (2.0 * bexp) + 3.0 5362 wcnd = dksat * factr2 ** expon 5363 ! 5364 return 5365 !................................... 5366 end subroutine wdfcnd 5367 !----------------------------------- 5368 5369 ! =========================== ! 5370 ! end contain programs ! 5371 ! =========================== ! 5372 5373 !................................... 5374 end subroutine sflx 5375 !----------------------------------- 5376