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, &
122 & tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm,z0, &
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, &
271 use machine ,
only : kind_phys
273 use physcons,
only : con_cp, con_rd, con_t0c, con_g, con_pi, &
274 & con_cliq, con_csol, con_hvap, con_hfus, &
285 integer,
parameter :: nsold = 4
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
294 real (kind=kind_phys),
parameter :: elcp = 2.4888e+3
296 real (kind=kind_phys),
parameter :: rd1 = 287.04
297 real (kind=kind_phys),
parameter :: cp = con_cp
298 real (kind=kind_phys),
parameter :: cp1 = 1004.5
299 real (kind=kind_phys),
parameter :: cp2 = 1004.0
301 real (kind=kind_phys),
parameter :: cph2o1 = 4.218e+3
302 real (kind=kind_phys),
parameter :: cph2o2 = 4.2e6
303 real (kind=kind_phys),
parameter :: cpice = con_csol
304 real (kind=kind_phys),
parameter :: cpice1 = 2.106e6
306 real (kind=kind_phys),
parameter :: sigma1 = 5.67e-8
309 integer,
intent(in) :: nsoil, couple, icein, vegtyp, soiltyp, &
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 &
317 logical,
intent(in) :: lheatstrg, exticeden
320 real (kind=kind_phys),
intent(inout) :: tbot, cmc, t1, sneqv, &
321 & stc(nsoil), smc(nsoil), sh2o(nsoil), ch, cm
324 integer,
intent(out) :: nroot
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, &
332 character(len=*),
intent(out) :: errmsg
333 integer,
intent(out) :: errflg
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
344 real (kind=kind_phys) :: shdfac0
345 real (kind=kind_phys),
dimension(nsold) :: rtdis, zsoil
347 logical :: frzgra, snowng
349 integer :: ice, k, kz
377 if(ivegsrc == 2)
then
378 if (vegtyp == 13)
then
384 if(ivegsrc == 1)
then
385 if (vegtyp == 15)
then
400 zsoil(kz) = -3.0 * float(kz) / float(nsoil)
409 zsoil(1) = -sldpth(1)
411 zsoil(kz) = -sldpth(kz) + zsoil(kz-1)
422 call redprm(errmsg, errflg)
423 if(ivegsrc == 1)
then
433 rsmin=400.0*(1-shdfac0)+40.0*shdfac0
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
463 bexp = bexp * max(1.+bexpp, 0.)
465 if( bexpp >= 0.)
then
466 bexp = bexp * min(1.+bexpp, 2.)
470 xlai = xlai * (1.+xlaip)
471 xlai = max(xlai, .75)
487 if (sneqv < 0.01)
then
493 elseif (ice == -1)
then
495 if (sneqv < 0.10)
then
519 if (sneqv .eq. 0.0)
then
524 sndens = sneqv / snowh
525 sndens = max( 0.0, min( 1.0, sndens ))
542 if (ffrozp > 0.)
then
545 if (t1 <= tfreez) frzgra = .true.
556 if (snowng .or. frzgra)
then
560 sn_new = ffrozp*prcp * dt * 0.001
561 sneqv = sneqv + sn_new
562 prcp1 = (1.-ffrozp)*prcp
566 sn_new = prcp * dt * 0.001
567 sneqv = sneqv + sn_new
611 if (sneqv == 0.0)
then
671 & ( smc(1), quartz, smcmax, sh2o(1), &
688 if((.not.lheatstrg) .and. &
689 & (ivegsrc == 1 .and. vegtyp == 13))
then
690 df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
692 df1 = df1 * exp( sbeta*shdfac )
701 dsoil = -0.5 * zsoil(1)
703 if (sneqv == 0.0)
then
705 ssoil = df1 * (t1 - stc(1)) / dsoil
710 frcsno = snowh / dtot
711 frcsoi = dsoil / dtot
721 df1a = frcsno*sncond + frcsoi*df1
728 df1 = df1a*sncovr + df1 *(1.0-sncovr)
734 ssoil = df1 * (t1 - stc(1)) / dtot
742 if (sncovr > 0.0)
then
756 t2v = sfctmp * (1.0 + 0.61*q2)
788 if (couple == 0)
then
793 t1v = t1 * (1.0 + 0.61 * q2)
794 th2v = th2 * (1.0 + 0.61 * q2)
834 if (shdfac > 0.)
then
854 if (sneqv .eq. 0.0)
then
896 sheat = -(ch*cp1*sfcprs) / (rd1*t2v) * (th2 - t1)
908 et(k) = et(k) * lsubc
912 esnow = esnow * lsubs
913 etp = etp * ((1.0 - sncovr)*lsubc + sncovr*lsubs)
916 eta = edir + ec + ett + esnow
940 runoff3 = runoff3 / dt
941 runoff2 = runoff2 + runoff3
950 runoff1 = snomlt / dt
956 soilm = -1.0 * smc(1) * zsoil(1)
958 soilm = soilm + smc(k)*(zsoil(k-1) - zsoil(k))
961 soilwm = -1.0 * (smcmax - smcwlt) * zsoil(1)
962 soilww = -1.0 * (smc(1) - smcwlt) * zsoil(1)
964 soilwm = soilwm + (smcmax - smcwlt) * (zsoil(k-1) - zsoil(k))
965 soilww = soilww + (smc(k) - smcwlt) * (zsoil(k-1) - zsoil(k))
968 soilw = soilww / soilwm
1033 albedo = alb + sncovr*(snoalb - alb)
1142 real (kind=kind_phys) :: delta, ff, gx, rr, part(nsold)
1159 ff = 0.55 * 2.0 * swdn / (rgl*xlai)
1160 rcs = (ff + rsmin/rsmax) / (1.0 + ff)
1161 rcs = max( rcs, 0.0001 )
1166 rct = 1.0 - 0.0016 * (topt - sfctmp)**2.0
1167 rct = max( rct, 0.0001 )
1172 rcq = 1.0 / (1.0 + hs*(q2sat-q2))
1173 rcq = max( rcq, 0.01 )
1178 gx = (sh2o(1) - smcwlt) / (smcref - smcwlt)
1179 gx = max( 0.0, min( 1.0, gx ) )
1182 part(1) = (zsoil(1)/zsoil(nroot)) * gx
1189 gx = (sh2o(k) - smcwlt) / (smcref - smcwlt)
1190 gx = max( 0.0, min( 1.0, gx ) )
1193 part(k) = ((zsoil(k) - zsoil(k-1)) / zsoil(nroot)) * gx
1201 rcsoil = rcsoil + part(k)
1203 rcsoil = max( rcsoil, 0.0001 )
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
1215 pc = (rr + delta) / (rr*(1.0 + rc*ch) + delta)
1252 real (kind=kind_phys),
parameter :: unit = 0.11631
1261 real (kind=kind_phys) :: c
1269 c = 0.328 * 10**(2.25*sndens)
1282 end subroutine csnow
1406 real (kind=kind_phys) :: df1, eta1, etp1, prcp1, yy, yynum, &
1407 & zz1, ec1, edir1, et1(nsoil), ett1
1438 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, &
1439 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, &
1440 & shdfac, cfactr, rtdis, fxexp, &
1442 & eta1, edir1, ec1, et1, ett1 &
1447 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, &
1448 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, &
1449 & edir1, ec1, et1, &
1453 & runoff1, runoff2, runoff3, drip &
1470 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, &
1471 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, &
1472 & edir1, ec1, et1, &
1476 & runoff1, runoff2, runoff3, drip &
1484 edir = edir1 * 1000.0
1488 et(k) = et1(k) * 1000.0
1495 if ( etp <= 0.0 )
then
1497 if ( etp < 0.0 )
then
1510 & ( smc(1), quartz, smcmax, sh2o(1), &
1528 if((.not.lheatstrg) .and. &
1529 & (ivegsrc == 1 .and. vegtyp == 13))
then
1530 df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
1532 df1 = df1 * exp( sbeta*shdfac )
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
1544 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, &
1545 & psisat, bexp, df1, ice, quartz, csoil, vegtyp, &
1546 & shdfac, lheatstrg, &
1548 & stc, t1, tbot, sh2o, &
1562 end subroutine nopac
1629 real (kind=kind_phys) :: a, delta, fnet, rad, rho
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)
1647 if (.not. snowng)
then
1648 if (prcp > 0.0) rr = rr + cph2o1*prcp/rch
1651 rr = rr + (cpice*ffrozp+cph2o1*(1.-ffrozp)) &
1655 fnet = fdown - sfcems*sigma1*t24 - ssoil
1661 flx2 = -lsubf * prcp
1667 rad = fnet/rch + th2 - sfctmp
1668 a = elcp * (q2sat - q2)
1669 epsca = (a*rr + rad*delta) / (delta + rr)
1670 etp = epsca * rch / lsubc
1870 character(len=*),
intent(out) :: errmsg
1871 integer,
intent(out) :: errflg
1875 real (kind=kind_phys) :: frzfact, frzk, refdk
1886 if (soiltyp > defined_soil)
then
1887 write(*,*)
'warning: too many soil types,soiltyp=',soiltyp, &
1888 &
'defined_soil=',defined_soil
1890 errmsg =
'ERROR(sflx.f): too many soil types'
1894 if (vegtyp > defined_veg)
then
1895 write(*,*)
'warning: too many veg types'
1897 errmsg =
'ERROR(sflx.f): too many veg types'
1901 if (slopetyp > defined_slope)
then
1902 write(*,*)
'warning: too many slope types'
1904 errmsg =
'ERROR(sflx.f): too many slope types'
1913 cfactr = cfactr_data
1914 cmcmax = cmcmax_data
1921 refkdt = refkdt_data
1928 dksat = satdk(soiltyp)
1929 dwsat = satdw(soiltyp)
1931 kdt = refkdt * dksat / refdk
1933 psisat = satpsi(soiltyp)
1934 quartz = qtz(soiltyp)
1935 smcdry = drysmc(soiltyp)
1936 smcmax = maxsmc(soiltyp)
1937 smcref = refsmc(soiltyp)
1938 smcwlt = wltsmc(soiltyp)
1940 frzfact = (smcmax / smcref) * (0.412 / 0.468)
1944 frzx = frzk * frzfact
1948 nroot = nroot_data(vegtyp)
1949 snup = snupx(vegtyp)
1950 rsmin = rsmtbl(vegtyp)
1952 rgl = rgltbl(vegtyp)
1956 xlai= lai_data(vegtyp)
1958 if (vegtyp == bare) shdfac = 0.0
1960 if (nroot > nsoil)
then
1961 write(*,*)
'warning: too many root layers'
1963 errmsg =
'ERROR(sflx.f): too many root layers'
1971 rtdis(i) = -sldpth(i) / zsoil(nroot)
1976 slope = slope_data(slopetyp)
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
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
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)
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, &
2061 integer :: ilech, itr
2065 real (kind=kind_phys) :: pslmu, pslms, pslhu, pslhs, zz
2066 real (kind=kind_phys) :: pspmu, pspms, psphu, psphs, xx, yy
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))
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
2094 zilfc = -czil * vkrm * sqvisc
2101 du2 = max( sfcspd*sfcspd, epsu2 )
2108 if (btgh*ch*dthv /= 0.0)
then
2109 wstar2 = wwst2 * abs( btgh*ch*dthv )**(2.0/3.0)
2114 ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust )
2118 zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0
2127 rlogu = log( zslu/zu )
2128 rlogt = log( zslt/zt )
2130 rlmo = elfc*ch*dthv / ustar**3
2142 zetalt = max( zslt*rlmo, ztmin )
2143 rlmo = zetalt / zslt
2144 zetalu = zslu * rlmo
2148 if (ilech == 0)
then
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
2156 xlu = sqrt( sqrt( xlu4 ) )
2157 xlt = sqrt( sqrt( xlt4 ) )
2158 xu = sqrt( sqrt( xu4 ) )
2159 xt = sqrt( sqrt( xt4 ) )
2169 simm = pspmu( xlu ) - psmz + rlogu
2171 simh = psphu( xlt ) - pshz + rlogt
2173 zetalu = min( zetalu, ztmax )
2174 zetalt = min( zetalt, ztmax )
2175 psmz = pspms( zetau )
2183 simm = pspms( zetalu ) - psmz + rlogu
2184 pshz = psphs( zetat )
2185 simh = psphs( zetalt ) - pshz + rlogt
2192 if (rlmo < 0.0)
then
2193 psmz = pslmu( zetau )
2201 simm = pslmu( zetalu ) - psmz + rlogu
2202 pshz = pslhu( zetat )
2203 simh = pslhu( zetalt ) - pshz + rlogt
2205 zetalu = min( zetalu, ztmax )
2206 zetalt = min( zetalt, ztmax )
2208 psmz = pslms( zetau )
2216 simm = pslms( zetalu ) - psmz + rlogu
2217 pshz = pslhs( zetat )
2218 simh = pslhs( zetalt ) - pshz + rlogt
2225 ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust )
2229 zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0
2232 rlogt = log( zslt/zt )
2234 ustark = ustar * vkrm
2235 cm = max( ustark/simm, cxch )
2236 ch = max( ustark/simh, cxch )
2240 if (btgh*ch*dthv /= 0.0)
then
2241 wstar2 = wwst2 * abs(btgh*ch*dthv) ** (2.0/3.0)
2246 rlmn = elfc*ch*dthv / ustar**3
2247 rlma = rlmo*wold + rlmn*wnew
2311 real (kind=kind_phys) :: rsnow, z0n
2319 if (sneqv < snup)
then
2320 rsnow = sneqv / snup
2321 sncovr = 1.0 - (exp(-salp*rsnow) - rsnow*exp(-salp))
2451 real,
parameter :: esdmin = 1.e-6
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, &
2502 prcp1 = prcp1 * 0.001
2532 etanrg = etp * ((1.0-sncovr)*lsubc + sncovr*lsubs)
2541 esnow1 = esnow * 0.001
2542 esnow2 = esnow1 * dt
2543 etanrg = esnow * lsubs
2547 if (sncovr < 1.0)
then
2551 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, &
2552 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, &
2553 & shdfac, cfactr, rtdis, fxexp, &
2555 & etns1, edir1, ec1, et1, ett1 &
2558 edir1 = edir1 * (1.0 - sncovr)
2559 ec1 = ec1 * (1.0 - sncovr)
2562 et1(k) = et1(k) * (1.0 - sncovr)
2565 ett1 = ett1 * (1.0 - sncovr)
2566 etns1 = etns1 * (1.0 - sncovr)
2568 edir = edir1 * 1000.0
2572 et(k) = et1(k) * 1000.0
2576 etns = etns1 * 1000.0
2580 esnow = etp * sncovr
2582 esnow1 = esnow * 0.001
2583 esnow2 = esnow1 * dt
2584 etanrg = esnow*lsubs + etns*lsubc
2598 flx1 = (cpice* ffrozp + cph2o1*(1.-ffrozp)) &
2599 & * prcp * (t1 - sfctmp)
2601 if (prcp > 0.0) flx1 = cph2o1 * prcp * (t1 - sfctmp)
2611 dsoil = -0.5 * zsoil(1)
2612 dtot = snowh + dsoil
2613 denom = 1.0 + df1 / (dtot * rr * rch)
2617 t12a = ( (fdown - flx1 - flx2 - sfcems*sigma1*t24) / rch &
2618 & + th2 - sfctmp - etanrg/rch ) / rr
2620 t12b = df1 * stc(1) / (dtot * rr * rch)
2621 t12 = (sfctmp + t12a + t12b) / denom
2632 if (t12 <= tfreez)
then
2635 ssoil = df1 * (t1 - stc(1)) / dtot
2637 sneqv = max(0.0, sneqv-esnow2)
2662 t1 = tfreez * max(0.01,sncovr**snoexp) + &
2663 & t12 * (1.0 - max(0.01,sncovr**snoexp))
2666 ssoil = df1 * (t1 - stc(1)) / dtot
2672 if (sneqv-esnow2 <= esdmin)
then
2684 sneqv = sneqv - esnow2
2685 seh = rch * (t1 - th2)
2690 flx3 = fdown - flx1 - flx2 - sfcems*sigma1*t14 &
2691 & - ssoil - seh - etanrg
2692 if (flx3 <= 0.0) flx3 = 0.0
2694 ex = flx3 * 0.001 / lsubf
2707 if (sneqv-snomlt >= esdmin)
then
2709 sneqv = sneqv - snomlt
2716 flx3 = ex * 1000.0 * lsubf
2733 if (ice == 0) prcp1 = prcp1 + ex
2750 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, &
2751 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, &
2752 & edir1, ec1, et1, &
2756 & runoff1, runoff2, runoff3, drip &
2769 yy = stc(1) - 0.5*ssoil*zsoil(1)*zz1 / df1
2780 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, &
2781 & psisat, bexp, df1, ice, quartz, csoil, vegtyp, &
2782 & shdfac, lheatstrg, &
2784 & stc, t11, tbot, sh2o, &
2794 if (sneqv > 0.0)
then
2798 & ( sneqv, dt, t1, yy, &
2819 elseif (ice == 1)
then
2821 if (sneqv >= 0.01)
then
2825 & ( sneqv, dt, t1, yy, &
2843 if (sneqv >= 0.10)
then
2847 & ( sneqv, dt, t1, yy, &
2912 real(kind=kind_phys) :: dsnew, snowhc, hnewc, newsnc, tempc
2919 snowhc = snowh * 100.0
2920 newsnc = sn_new * 100.0
2921 tempc = sfctmp - tfreez
2928 if(.not. exticeden)
then
2929 if (tempc <= -15.0)
then
2932 dsnew = 0.05 + 0.0017*(tempc + 15.0)**1.5
2935 dsnew = rhonewsn*0.001
2940 hnewc = newsnc / dsnew
2941 sndens = (snowhc*sndens + hnewc*dsnew) / (snowhc + hnewc)
2942 snowhc = snowhc + hnewc
2943 snowh = snowhc * 0.01
2986 real(kind=kind_phys) :: z0s
2994 z0 = (1.0 - sncovr)*z0 + sncovr*z0s
3009 & ( smc, qz, smcmax, sh2o, &
3059 real (kind=kind_phys),
intent(in) :: smc, qz, smcmax, sh2o
3062 real (kind=kind_phys),
intent(out) :: df
3065 real (kind=kind_phys) :: gammd, thkdry, ake, thkice, thko, &
3066 & thkqtz, thksat, thks, thkw, satratio, xu, xunfroz
3075 satratio = smc / smcmax
3086 thks = (thkqtz**qz) * (thko**(1.0-qz))
3090 xunfroz = (sh2o + 1.e-9) / (smc + 1.e-9)
3098 thksat = thks**(1.-smcmax) * thkice**(smcmax-xu) * thkw**(xu)
3102 gammd = (1.0 - smcmax) * 2700.0
3106 thkdry = (0.135*gammd + 64.7) / (2700.0 - 0.947*gammd)
3108 if ( sh2o+0.0005 < smc )
then
3115 if ( satratio > 0.1 )
then
3121 ake = log10( satratio ) + 1.0
3134 df = ake * (thksat - thkdry) + thkdry
3156 & ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil, &
3157 & sh2o, smcmax, smcwlt, smcref, smcdry, pc, &
3158 & shdfac, cfactr, rtdis, fxexp, &
3160 & eta1, edir1, ec1, et1, ett1 &
3206 integer,
intent(in) :: nsoil, nroot
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)
3213 real (kind=kind_phys),
intent(out) :: eta1, edir1, ec1, ett1, &
3217 real (kind=kind_phys) :: cmc2ms
3235 if (etp1 > 0.0)
then
3241 if (shdfac < 1.0)
then
3245 & ( etp1, sh2o(1), shdfac, smcmax, smcdry, fxexp, &
3255 if (shdfac > 0.0)
then
3259 & ( nsoil, nroot, etp1, sh2o, smcwlt, smcref, &
3260 & cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis, &
3266 ett1 = ett1 + et1(k)
3273 ec1 = shdfac * ( (cmc/cmcmax)**cfactr ) * etp1
3282 ec1 = min( cmc2ms, ec1 )
3289 eta1 = edir1 + ett1 + ec1
3294 end subroutine evapo
3305 & ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot, &
3306 & psisat, bexp, df1, ice, quartz, csoil, vegtyp, &
3307 & shdfac, lheatstrg, &
3309 & stc, t1, tbot, sh2o, &
3358 real (kind=kind_phys),
parameter :: ctfil1 = 0.5
3359 real (kind=kind_phys),
parameter :: ctfil2 = 1.0 - ctfil1
3362 integer,
intent(in) :: nsoil, ice, vegtyp
3364 real (kind=kind_phys),
intent(in) :: smc(nsoil), smcmax, dt, yy, &
3365 & zz1, zsoil(nsoil), zbot, psisat, bexp, df1, quartz, csoil, &
3368 logical,
intent(in) :: lheatstrg
3371 real (kind=kind_phys),
intent(inout) :: stc(nsoil), t1, tbot, &
3375 real (kind=kind_phys),
intent(out) :: ssoil
3378 real (kind=kind_phys) :: ai(nsold), bi(nsold), ci(nsold), oldt1, &
3379 & rhsts(nsold), stcf(nsold), stsoil(nsoil)
3399 & ( nsoil, stc, zsoil, yy, zz1, df1, ice, &
3403 & rhsts, ai, bi, ci &
3408 & ( nsoil, stc, dt, &
3410 & rhsts, ai, bi, ci, &
3421 & ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, &
3422 & zbot, psisat, dt, bexp, df1, quartz, csoil,vegtyp, &
3423 & shdfac, lheatstrg, &
3427 & rhsts, ai, bi, ci &
3432 & ( nsoil, stc, dt, &
3434 & rhsts, ai, bi, ci, &
3451 t1 = (yy + (zz1 - 1.0)*stc(1)) / zz1
3452 t1 = ctfil1*t1 + ctfil2*oldt1
3455 stc(i) = ctfil1*stc(i) + ctfil2*stsoil(i)
3460 ssoil = df1*(stc(1) - t1) / (0.5*zsoil(1))
3465 end subroutine shflx
3479 & ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1, &
3480 & zsoil, slope, frzx, bexp, dksat, dwsat, shdfac, &
3481 & edir1, ec1, et1, &
3485 & runoff1, runoff2, runoff3, drip &
3535 integer,
intent(in) :: nsoil
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)
3542 real (kind=kind_phys),
intent(inout) :: cmc, sh2o(nsoil), &
3546 real (kind=kind_phys),
intent(out) :: runoff1, runoff2, &
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)
3564 rhsct = shdfac*prcp1 - ec1
3572 excess = cmc + trhsct
3574 if (excess > cmcmax) drip = excess - cmcmax
3579 pcpdrp = (1.0 - shdfac)*prcp1 + drip/dt
3584 sice(i) = smc(i) - sh2o(i)
3608 if ( (pcpdrp*dt) > (0.001*1000.0*(-zsoil(1))*smcmax) )
then
3617 & ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat, &
3618 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, &
3620 & rhstt, runoff1, runoff2, ai, bi, ci &
3625 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, &
3627 & dummy, rhstt, ai, bi, ci, &
3629 & sh2ofg, runoff3, smc &
3633 sh2oa(k) = (sh2o(k) + sh2ofg(k)) * 0.5
3638 & ( nsoil, edir1, et1, sh2o, sh2oa, pcpdrp, zsoil, dwsat, &
3639 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, &
3641 & rhstt, runoff1, runoff2, ai, bi, ci &
3646 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, &
3648 & cmc, rhstt, ai, bi, ci, &
3650 & sh2o, runoff3, smc &
3657 & ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat, &
3658 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, &
3660 & rhstt, runoff1, runoff2, ai, bi, ci &
3665 & ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice, &
3667 & cmc, rhstt, ai, bi, ci, &
3669 & sh2o, runoff3, smc &
3678 end subroutine smflx
3691 & ( esd, dtsec, tsnow, tsoil, &
3725 real (kind=kind_phys),
parameter :: c1 = 0.01
3726 real (kind=kind_phys),
parameter :: c2 = 21.0
3729 real (kind=kind_phys),
intent(in) :: esd, dtsec, tsnow, tsoil
3732 real (kind=kind_phys),
intent(inout) :: snowh, sndens
3735 real (kind=kind_phys) :: bfac, dsx, dthr, dw, snowhc, pexp, &
3736 & tavgc, tsnowc, tsoilc, esdc, esdcx
3744 snowhc = snowh * 100.0
3746 dthr = dtsec / 3600.0
3747 tsnowc = tsnow - tfreez
3748 tsoilc = tsoil - tfreez
3752 tavgc = 0.5 * (tsnowc + tsoilc)
3762 if (esdc > 1.e-2)
then
3768 bfac = dthr*c1 * exp(0.08*tavgc - c2*sndens)
3793 pexp = (1.0 + pexp)*bfac*esdcx/real(j+1)
3818 dsx = max( min( dsx, 0.40 ), 0.05 )
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
3834 snowhc = esdc / sndens
3835 snowh = snowhc * 0.01
3854 & ( etp1, smc, shdfac, smcmax, smcdry, fxexp, &
3884 real (kind=kind_phys),
intent(in) :: etp1, smc, shdfac, smcmax, &
3888 real (kind=kind_phys),
intent(out) :: edir1
3891 real (kind=kind_phys) :: fx, sratio
3900 sratio = (smc - smcdry) / (smcmax - smcdry)
3902 if (sratio > 0.0)
then
3904 fx = max( min( fx, 1.0 ), 0.0 )
3911 edir1 = fx * ( 1.0 - shdfac ) * etp1
3915 end subroutine devap
3934 & ( tkelv, smc, sh2o, smcmax, bexp, psis, &
3974 real (kind=kind_phys),
parameter :: ck = 8.0
3976 real (kind=kind_phys),
parameter :: blim = 5.5
3977 real (kind=kind_phys),
parameter ::
error = 0.005
3980 real (kind=kind_phys),
intent(in) :: tkelv, smc, sh2o, smcmax, &
3984 real (kind=kind_phys),
intent(out) :: liqwat
3987 real (kind=kind_phys) :: bx, denom, df, dswl, fk, swl, swlk
3989 integer :: nlog, kcount
3998 if (bexp > blim) bx = blim
4007 if (tkelv > (tfreez-1.e-3))
then
4024 swl = max( min( swl, smc-0.02 ), 0.0 )
4028 do while ( (nlog < 10) .and. (kcount == 0) )
4031 df = log( (psis*gs2/lsubf) * ( (1.0 + ck*swl)**2.0 ) &
4032 & * (smcmax/(smc-swl))**bx ) - log(-(tkelv-tfreez)/tkelv)
4034 denom = 2.0*ck/(1.0 + ck*swl) + bx/(smc - swl)
4035 swlk = swl - df/denom
4039 swlk = max( min( swlk, smc-0.02 ), 0.0 )
4043 dswl = abs(swlk - swl)
4049 if ( dswl <=
error )
then
4064 if (kcount == 0)
then
4065 fk = ( ( (lsubf/(gs2*(-psis))) &
4066 & * ((tkelv-tfreez)/tkelv) )**(-1/bx) ) * smcmax
4068 fk = max( fk, 0.02 )
4070 liqwat = min( fk, smc )
4077 end subroutine frh2o
4089 & ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot, &
4090 & zbot, psisat, dt, bexp, df1, quartz, csoil, vegtyp, &
4091 & shdfac, lheatstrg, &
4095 & rhsts, ai, bi, ci &
4144 integer,
intent(in) :: nsoil, vegtyp
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
4150 logical,
intent(in) :: lheatstrg
4153 real (kind=kind_phys),
intent(inout) :: sh2o(nsoil)
4156 real (kind=kind_phys),
intent(out) :: rhsts(nsoil), ai(nsold), &
4157 & bi(nsold), ci(nsold)
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
4173 if (.not.lheatstrg .and. ivegsrc == 1)
then
4178 if( vegtyp == 13 )
then
4180 csoil_loc=3.0e6*(1.-shdfac)+csoil*shdfac
4193 hcpct = sh2o(1)*cph2o2 + (1.0 - smcmax)*csoil_loc &
4194 & + (smcmax - smc(1))*cp2 + (smc(1) - sh2o(1))*cpice1
4198 ddz = 1.0 / ( -0.5*zsoil(2) )
4200 ci(1) = (df1*ddz) / ( zsoil(1)*hcpct )
4201 bi(1) = -ci(1) + df1 / ( 0.5*zsoil(1)*zsoil(1)*hcpct*zz1 )
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)
4216 qtot = ssoil - df1*dtsdz
4229 tsurf = (yy + (zz1-1)*stc(1)) / zz1
4233 & ( stc(1), stc(2), zsoil, zbot, 1, nsoil, &
4242 sice = smc(1) - sh2o(1)
4249 if ( (sice > 0.0) .or. (tsurf < tfreez) .or. &
4250 & (stc(1) < tfreez) .or. (tbk < tfreez) )
then
4256 & ( tsurf, stc(1), tbk, zsoil, nsoil, 1, &
4269 & ( nsoil, 1, tavg, smc(1), smcmax, psisat, bexp, dt, &
4278 rhsts(1) = rhsts(1) - tsnsr / ( zsoil(1)*hcpct )
4297 hcpct = sh2o(k)*cph2o2 + (1.0 - smcmax)*csoil_loc &
4298 & + (smcmax - smc(k))*cp2 + (smc(k) - sh2o(k))*cpice1
4300 if (k /= nsoil)
then
4307 & ( smc(k), quartz, smcmax, sh2o(k), &
4319 if((.not.lheatstrg) .and.
4320 & (ivegsrc == 1 .and. vegtyp == 13))
then
4321 df1n = 3.24*(1.-shdfac) + shdfac*df1n
4326 denom = 0.5 * (zsoil(k-1) - zsoil(k+1))
4327 dtsdz2 = (stc(k) - stc(k+1)) / denom
4331 ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1))
4332 ci(k) = -df1n*ddz2 / ((zsoil(k-1) - zsoil(k)) * hcpct)
4341 & ( stc(k), stc(k+1), zsoil, zbot, k, nsoil, &
4355 & ( smc(k), quartz, smcmax, sh2o(k), &
4367 if((.not.lheatstrg) .and.
4368 & (ivegsrc == 1 .and. vegtyp == 13))
then
4369 df1n = 3.24*(1.-shdfac) + shdfac*df1n
4374 denom = 0.5 * (zsoil(k-1) + zsoil(k)) - zbot
4375 dtsdz2 = (stc(k) - tbot) / denom
4388 & ( stc(k), tbot, zsoil, zbot, k, nsoil, &
4399 denom = (zsoil(k) - zsoil(k-1)) * hcpct
4400 rhsts(k) = ( df1n*dtsdz2 - df1k*dtsdz ) / denom
4402 qtot = -1.0 * denom * rhsts(k)
4403 sice = smc(k) - sh2o(k)
4405 if ( (sice > 0.0) .or. (tbk < tfreez) .or. &
4406 & (stc(k) < tfreez) .or. (tbk1 < tfreez) )
then
4412 & ( tbk, stc(k), tbk1, zsoil, nsoil, k, &
4423 & ( nsoil, k, tavg, smc(k), smcmax, psisat, bexp, dt, &
4431 rhsts(k) = rhsts(k) - tsnsr/denom
4436 ai(k) = - df1 * ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
4437 bi(k) = -(ai(k) + ci(k))
4462 & ( nsoil, stc, zsoil, yy, zz1, df1, ice, &
4466 & rhsts, ai, bi, ci &
4505 integer,
intent(in) :: nsoil, ice
4507 real (kind=kind_phys),
intent(in) :: stc(nsoil), zsoil(nsoil), &
4511 real (kind=kind_phys),
intent(inout) :: tbot
4514 real (kind=kind_phys),
intent(out) :: rhsts(nsoil), ai(nsold), &
4515 & bi(nsold), ci(nsold)
4518 real (kind=kind_phys) :: ddz, ddz2, denom, dtsdz, dtsdz2, &
4519 & hcpct, ssoil, zbot
4561 ddz = 1.0 / (-0.5*zsoil(2))
4563 ci(1) = (df1*ddz) / (zsoil(1)*hcpct)
4564 bi(1) = -ci(1) + df1 / (0.5*zsoil(1)*zsoil(1)*hcpct*zz1)
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)
4582 if (k /= nsoil)
then
4586 denom = 0.5 * (zsoil(k-1) - zsoil(k+1))
4587 dtsdz2 = (stc(k) - stc(k+1)) / denom
4591 ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1))
4592 ci(k) = -df1*ddz2 / ((zsoil(k-1) - zsoil(k))*hcpct)
4598 dtsdz2 = (stc(k) - tbot) &
4599 & / (0.5*(zsoil(k-1) + zsoil(k)) - zbot)
4609 denom = (zsoil(k) - zsoil(k-1)) * hcpct
4610 rhsts(k) = (df1*dtsdz2 - df1*dtsdz) / denom
4614 ai(k) = - df1*ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
4615 bi(k) = -(ai(k) + ci(k))
4635 & ( nsoil, stcin, dt, &
4637 & rhsts, ai, bi, ci, &
4669 integer,
intent(in) :: nsoil
4671 real (kind=kind_phys),
intent(in) :: stcin(nsoil), dt
4674 real (kind=kind_phys),
intent(inout) :: rhsts(nsoil), &
4675 & ai(nsold), bi(nsold), ci(nsold)
4678 real (kind=kind_phys),
intent(out) :: stcout(nsoil)
4683 real (kind=kind_phys) :: ciin(nsold), rhstsin(nsoil)
4691 rhsts(k) = rhsts(k) * dt
4693 bi(k) = 1.0 + bi(k)*dt
4700 rhstsin(k) = rhsts(k)
4711 & ( nsoil, ai, bi, rhstsin, &
4721 stcout(k) = stcin(k) + ci(k)
4726 end subroutine hstep
4735 & ( nsoil, a, b, d, &
4783 integer,
intent(in) :: nsoil
4785 real (kind=kind_phys),
dimension(nsoil),
intent(in) :: a, b, d
4788 real (kind=kind_phys),
dimension(nsoil),
intent(inout) :: c
4791 real (kind=kind_phys),
dimension(nsoil),
intent(out) :: p, delta
4806 delta(1) = d(1) / b(1)
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)) )
4818 p(nsoil) = delta(nsoil)
4824 p(kk) = p(kk)*p(kk+1) + delta(kk)
4838 & ( nsoil, k, tavg, smc, smcmax, psisat, bexp, dt, &
4878 real (kind=kind_phys),
parameter :: dh2o = 1.0000e3
4881 integer,
intent(in) :: nsoil, k
4883 real (kind=kind_phys),
intent(in) :: tavg, smc, smcmax, psisat, &
4884 & bexp, dt, qtot, zsoil(nsoil)
4887 real (kind=kind_phys),
intent(inout) :: sh2o
4890 real (kind=kind_phys),
intent(out) :: tsrc
4893 real (kind=kind_phys) :: dz, free, xh2o
4903 dz = zsoil(k-1) - zsoil(k)
4916 & ( tavg, smc, sh2o, smcmax, bexp, psisat, &
4930 xh2o = sh2o + qtot*dt / (dh2o*lsubf*dz)
4936 if ( xh2o < sh2o .and. xh2o < free)
then
4937 if ( free > sh2o )
then
4948 if ( xh2o > sh2o .and. xh2o > free )
then
4949 if ( free < sh2o )
then
4956 xh2o = max( min( xh2o, smc ), 0.0 )
4961 tsrc = -dh2o * lsubf * dz * (xh2o - sh2o) / dt
4978 & ( nsoil, edir, et, sh2o, sh2oa, pcpdrp, zsoil, dwsat, &
4979 & dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice, &
4981 & rhstt, runoff1, runoff2, ai, bi, ci &
5026 integer,
intent(in) :: nsoil
5028 real (kind=kind_phys),
dimension(nsoil),
intent(in) :: et, &
5029 & sh2o, sh2oa, zsoil, sice
5031 real (kind=kind_phys),
intent(in) :: edir, pcpdrp, dwsat, dksat, &
5032 & smcmax, smcwlt, bexp, dt, slope, kdt, frzx
5035 real (kind=kind_phys),
intent(out) :: runoff1, runoff2, &
5036 & rhstt(nsoil), ai(nsold), bi(nsold), ci(nsold)
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)
5045 integer :: cvfrz, ialp1, iohinf, j, jj, k, ks
5058 parameter(cvfrz = 3)
5060c ----------------------------------------------------------------------
5072 if (sice(ks) > sicemax) sicemax = sice(ks)
5080 if (pcpdrp /= 0.0)
then
5085 smcav = smcmax - smcwlt
5086 dmax(1) = -zsoil(1) * smcav
5090 dice = -zsoil(1) * sice(1)
5092 dmax(1) = dmax(1)*(1.0 - (sh2oa(1)+sice(1)-smcwlt)/smcav)
5099 dice = dice + ( zsoil(ks-1) - zsoil(ks) ) * sice(ks)
5101 dmax(ks) = (zsoil(ks-1)-zsoil(ks))*smcav
5102 dmax(ks) = dmax(ks)*(1.0 - (sh2oa(ks)+sice(ks)-smcwlt)/smcav)
5109 val = 1.0 - exp(-kdt*dt1)
5113 if (px < 0.0) px = 0.0
5115 infmax = (px*(ddt/(px+ddt)))/dt
5122 if (dice > 1.e-2)
then
5123 acrt = cvfrz * frzx / dice
5134 sum = sum + (acrt**( cvfrz-j)) / float(k)
5137 fcr = 1.0 - exp(-acrt) * sum
5140 infmax = infmax * fcr
5151 & ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax, &
5156 infmax = max( infmax, wcnd )
5157 infmax = min( infmax, px )
5159 if (pcpdrp > infmax)
then
5160 runoff1 = pcpdrp - infmax
5174 & ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax, &
5181 ddz = 1.0 / ( -.5*zsoil(2) )
5183 bi(1) = wdf * ddz / ( -zsoil(1) )
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)
5200 denom2 = (zsoil(k-1) - zsoil(k))
5202 if (k /= nsoil)
then
5213 & ( mxsmc2, smcmax, bexp, dksat, dwsat, sicemax, &
5220 denom = (zsoil(k-1) - zsoil(k+1))
5221 dsmdz2 = (sh2o(k) - sh2o(k+1)) / (denom * 0.5)
5226 ci(k) = -wdf2 * ddz2 / denom2
5239 & ( sh2oa(nsoil), smcmax, bexp, dksat, dwsat, sicemax, &
5255 numer = wdf2*dsmdz2 + slopx*wcnd2 - wdf*dsmdz - wcnd + et(k)
5256 rhstt(k) = numer / (-denom2)
5260 ai(k) = -wdf * ddz / denom2
5261 bi(k) = -( ai(k) + ci(k) )
5266 if (k == nsoil)
then
5267 runoff2 = slopx * wcnd2
5270 if (k /= nsoil)
then
5290 & ( nsoil, sh2oin, rhsct, dt, smcmax, cmcmax, zsoil, sice, &
5292 & cmc, rhstt, ai, bi, ci, &
5294 & sh2oout, runoff3, smc &
5332 integer,
intent(in) :: nsoil
5334 real (kind=kind_phys),
dimension(nsoil),
intent(in) :: sh2oin, &
5337 real (kind=kind_phys),
intent(in) :: rhsct, dt, smcmax, cmcmax
5340 real (kind=kind_phys),
intent(inout) :: cmc, rhstt(nsoil), &
5341 & ai(nsold), bi(nsold), ci(nsold)
5344 real (kind=kind_phys),
intent(out) :: sh2oout(nsoil), runoff3, &
5348 real (kind=kind_phys) :: ciin(nsold), rhsttin(nsoil), ddz, stot, &
5351 integer :: i, k, kk11
5359 rhstt(k) = rhstt(k) * dt
5361 bi(k) = 1. + bi(k) * dt
5368 rhsttin(k) = rhstt(k)
5379 & ( nsoil, ai, bi, rhsttin, &
5395 if (k /= 1) ddz = zsoil(k - 1) - zsoil(k)
5397 sh2oout(k) = sh2oin(k) + ci(k) + wplus/ddz
5399 stot = sh2oout(k) + sice(k)
5400 if (stot > smcmax)
then
5405 ddz = -zsoil(k) + zsoil(kk11)
5408 wplus = (stot - smcmax) * ddz
5413 smc(k) = max( min( stot, smcmax ), 0.02 )
5414 sh2oout(k) = max( smc(k)-sice(k), 0.0 )
5422 cmc = cmc + dt*rhsct
5423 if (cmc < 1.e-20) cmc = 0.0
5424 cmc = min( cmc, cmcmax )
5428 end subroutine sstep
5438 & ( tu, tb, zsoil, zbot, k, nsoil, &
5467 integer,
intent(in) :: k, nsoil
5469 real (kind=kind_phys),
intent(in) :: tu, tb, zbot, zsoil(nsoil)
5472 real (kind=kind_phys),
intent(out) :: tbnd1
5475 real (kind=kind_phys) :: zb, zup
5488 if (k == nsoil)
then
5489 zb = 2.0*zbot - zsoil(k)
5496 tbnd1 = tu + (tb-tu)*(zup-zsoil(k))/(zup-zb)
5513 & ( tup, tm, tdn, zsoil, nsoil, k, &
5545 integer,
intent(in) :: nsoil, k
5547 real (kind=kind_phys),
intent(in) :: tup, tm, tdn, zsoil(nsoil)
5550 real (kind=kind_phys),
intent(out) :: tavg
5553 real (kind=kind_phys) :: dz, dzh, x0, xdn, xup
5560 dz = zsoil(k-1) - zsoil(k)
5565 if (tup < tfreez)
then
5567 if (tm < tfreez)
then
5568 if (tdn < tfreez)
then
5569 tavg = (tup + 2.0*tm + tdn) / 4.0
5571 x0 = (tfreez - tm) * dzh / (tdn - tm)
5572 tavg = 0.5*(tup*dzh + tm*(dzh+x0)+tfreez*(2.*dzh-x0)) / dz
5575 if (tdn < tfreez)
then
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
5580 xup = (tfreez-tup) * dzh / (tm-tup)
5581 tavg = 0.5*(tup*xup + tfreez*(2.*dz-xup)) / dz
5587 if (tm < tfreez)
then
5588 if (tdn < tfreez)
then
5589 xup = dzh - (tfreez-tup) * dzh / (tm-tup)
5590 tavg = 0.5*(tfreez*(dz-xup) + tm*(dzh+xup)+tdn*dzh) / dz
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
5597 if (tdn < tfreez)
then
5598 xdn = dzh - (tfreez-tm) * dzh / (tdn-tm)
5599 tavg = (tfreez*(dz-xdn) + 0.5*(tfreez+tdn)*xdn) / dz
5601 tavg = (tup + 2.0*tm + tdn) / 4.0
5618 & ( nsoil, nroot, etp1, smc, smcwlt, smcref, &
5619 & cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis, &
5654 integer,
intent(in) :: nsoil, nroot
5656 real (kind=kind_phys),
intent(in) :: etp1, smcwlt, smcref, &
5657 & cmc, cmcmax, shdfac, pc, cfactr
5659 real (kind=kind_phys),
dimension(nsoil),
intent(in) :: smc, &
5663 real (kind=kind_phys),
dimension(nsoil),
intent(out) :: et1
5666 real (kind=kind_phys) :: denom, etp1a, rtx, sgx, gx(7)
5684 if (cmc /= 0.0)
then
5685 etp1a = shdfac * pc * etp1 * (1.0 - (cmc /cmcmax) ** cfactr)
5687 etp1a = shdfac * pc * etp1
5692 gx(i) = ( smc(i) - smcwlt ) / ( smcref - smcwlt )
5693 gx(i) = max( min( gx(i), 1.0 ), 0.0 )
5700 rtx = rtdis(i) + gx(i) - sgx
5701 gx(i) = gx(i) * max( rtx, 0.0 )
5702 denom = denom + gx(i)
5704 if (denom <= 0.0) denom = 1.0
5707 et1(i) = etp1a * gx(i) / denom
5752 & ( smc, smcmax, bexp, dksat, dwsat, sicemax, &
5782 real (kind=kind_phys),
intent(in) :: smc, smcmax, bexp, dksat, &
5786 real (kind=kind_phys),
intent(out) :: wdf, wcnd
5789 real (kind=kind_phys) :: expon, factr1, factr2, vkwgt
5795 factr1 = min(1.0, max(0.0, 0.2/smcmax))
5796 factr2 = min(1.0, max(0.0, smc/smcmax))
5801 wdf = dwsat * factr2 ** expon
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
5821 expon = (2.0 * bexp) + 3.0
5822 wcnd = dksat * factr2 ** expon
subroutine csnow
This subroutine calculates snow termal conductivity.
subroutine alcalc
This subroutine calculates albedo including snow effect (0 -> 1).
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.
subroutine redprm(errmsg, errflg)
This subroutine internally sets default values or optionally read-in via namelist i/o,...
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.
subroutine penman
This subroutine calculates potential evaporation for the current point. various partial sums/products...
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...
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...
subroutine rosr12(nsoil, a, b, d, c, p, delta)
This subroutine inverts (solve) the tri-diagonal matrix problem.
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...
subroutine transp(nsoil, nroot, etp1, smc, smcwlt, smcref, cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis, et1)
This subroutine calculates transpiration for the veg class.
subroutine frh2o(tkelv, smc, sh2o, smcmax, bexp, psis, liqwat)
This subroutine calculates amount of supercooled liquid soil water content if temperature is below 27...
subroutine tdfcnd(smc, qz, smcmax, sh2o, df)
This subroutine calculates thermal diffusivity and conductivity of the soil for a given point and tim...
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...
subroutine nopac
This subroutine calculates soil moisture and heat flux values and update soil moisture content and so...
subroutine snowz0
This subroutine calculates total roughness length over snow.
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...
subroutine canres
This subroutine calculates canopy resistance which depends on incoming solar radiation,...
subroutine snowpack(esd, dtsec, tsnow, tsoil, snowh, sndens)
This subroutine calculates compaction of a snowpack under conditions of increasing snow density,...
subroutine snow_new
This subroutine calculates snow depth and densitity to account for the new snowfall....
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...
subroutine devap(etp1, smc, shdfac, smcmax, smcdry, fxexp, edir1)
This subrtouine calculates direct soil evaporation.
subroutine wdfcnd(smc, smcmax, bexp, dksat, dwsat, sicemax, wdf, wcnd)
This subroutine calculates soil water diffusivity and soil hydraulic conductivity.
subroutine hstep(nsoil, stcin, dt, rhsts, ai, bi, ci, stcout)
This subroutine calculates/updates the soil temperature field.
subroutine snopac
This subroutine calculates soil moisture and heat flux values and update soil moisture content and so...
subroutine tmpavg(tup, tm, tdn, zsoil, nsoil, k, tavg)
This subroutine calculates soil layer average temperature (tavg) in freezing/thawing layer using up,...
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 ...
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...
subroutine snfrac
This subroutine calculates snow fraction (0->1).
subroutine sfcdif
This subroutine calculates surface layer exchange coefficients via iterative process(see Chen et al....
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.