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