! ===================================================================== !
! description: !
! !
! usage: !
! !
! call sfc_drv !
! --- inputs: !
! ( im, km, ps, u1, v1, t1, q1, soiltyp, vegtype, sigmaf, !
! sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, !
! prsl1, prslki, zf, islimsk, ddvel, slopetyp, !
! shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, !
! --- in/outs: !
! weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, !
! canopy, trans, tsurf, !
! --- outputs: !
! sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, !
! cmm, chh, evbs, evcw, sbsno, snowc, stm, snohf, !
! smcwlt2, smcref2, zorl, wet1 ) !
! !
! !
! subprogram called: sflx !
! !
! program history log: !
! xxxx -- created !
! 200x -- sarah lu modified !
! oct 2006 -- h. wei modified !
! apr 2009 -- y.-t. hou modified to include surface emissivity !
! effect on lw radiation. replaced the comfussing !
! slrad (net sw + dlw) with sfc net sw snet=dsw-usw !
! sep 2009 -- s. moorthi modification to remove rcl and unit change!
! nov 2011 -- sarah lu corrected wet1 calculation
! !
! ==================== defination of variables ==================== !
! !
! inputs: size !
! im - integer, horiz dimention and num of used pts 1 !
! km - integer, vertical soil layer dimension 1 !
! ps - real, surface pressure (pa) im !
! u1, v1 - real, u/v component of surface layer wind im !
! t1 - real, surface layer mean temperature (k) im !
! q1 - real, surface layer mean specific humidity im !
! soiltyp - integer, soil type (integer index) im !
! vegtype - integer, vegetation type (integer index) im !
! sigmaf - real, areal fractional cover of green vegetation im !
! sfcemis - real, sfc lw emissivity ( fraction ) im !
! dlwflx - real, total sky sfc downward lw flux ( w/m**2 ) im !
! dswflx - real, total sky sfc downward sw flux ( w/m**2 ) im !
! snet - real, total sky sfc netsw flx into ground(w/m**2) im !
! delt - real, time interval (second) 1 !
! tg3 - real, deep soil temperature (k) im !
! cm - real, surface exchange coeff for momentum (m/s) im !
! ch - real, surface exchange coeff heat & moisture(m/s) im !
! prsl1 - real, sfc layer 1 mean pressure (pa) im !
! prslki - real, im !
! zf - real, height of bottom layer (m) im !
! islimsk - integer, sea/land/ice mask (=0/1/2) im !
! ddvel - real, im !
! slopetyp - integer, class of sfc slope (integer index) im !
! shdmin - real, min fractional coverage of green veg im !
! shdmax - real, max fractnl cover of green veg (not used) im !
! snoalb - real, upper bound on max albedo over deep snow im !
! sfalb - real, mean sfc diffused sw albedo (fractional) im !
! flag_iter- logical, im !
! flag_guess-logical, im !
! !
! input/outputs: !
! weasd - real, water equivalent accumulated snow depth (mm) im !
! snwdph - real, snow depth (water equiv) over land im !
! tskin - real, ground surface skin temperature ( k ) im !
! tprcp - real, total precipitation im !
! srflag - real, snow/rain flag for precipitation im !
! smc - real, total soil moisture content (fractional) im,km !
! stc - real, soil temp (k) im,km !
! slc - real, liquid soil moisture im,km !
! canopy - real, canopy moisture content (m) im !
! trans - real, total plant transpiration (m/s) im !
! tsurf - real, surface skin temperature (after iteration) im !
! !
! outputs: !
! sncovr1 - real, snow cover over land (fractional) im !
! qsurf - real, specific humidity at sfc im !
! gflux - real, soil heat flux (w/m**2) im !
! drain - real, subsurface runoff (mm/s) im !
! evap - real, evaperation from latent heat flux im !
! hflx - real, sensible heat flux im !
! ep - real, potential evaporation im !
! runoff - real, surface runoff (m/s) im !
! cmm - real, im !
! chh - real, im !
! evbs - real, direct soil evaporation (m/s) im !
! evcw - real, canopy water evaporation (m/s) im !
! sbsno - real, sublimation/deposit from snopack (m/s) im !
! snowc - real, fractional snow cover im !
! stm - real, total soil column moisture content (m) im !
! snohf - real, snow/freezing-rain latent heat flux (w/m**2)im !
! smcwlt2 - real, dry soil moisture threshold im !
! smcref2 - real, soil moisture threshold im !
! zorl - real, surface roughness im !
! wet1 - real, normalized soil wetness im !
! !
! ==================== end of description ===================== !
!-----------------------------------
subroutine sfc_drv &
!...................................
! --- inputs:
& ( im, km, ps, u1, v1, t1, q1, soiltyp, vegtype, sigmaf, &
& sfcemis, dlwflx, dswsfc, snet, delt, tg3, cm, ch, &
& prsl1, prslki, zf, islimsk, ddvel, slopetyp, &
& shdmin, shdmax, snoalb, sfalb, flag_iter, flag_guess, &
! --- in/outs:
& weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, &
& canopy, trans, tsurf, &
! --- outputs:
& sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, &
& cmm, chh, evbs, evcw, sbsno, snowc, stm, snohf, &
& smcwlt2, smcref2, zorl, wet1 &
& )
!
use machine , only : kind_phys
use funcphys, only : fpvs
use physcons, only : con_g, con_hvap, con_cp, con_jcal, &
& con_eps, con_epsm1, con_fvirt, con_rd
implicit none
! --- constant parameters:
real(kind=kind_phys), parameter :: cpinv = 1.0/con_cp
real(kind=kind_phys), parameter :: hvapi = 1.0/con_hvap
real(kind=kind_phys), parameter :: elocp = con_hvap/con_cp
real(kind=kind_phys), parameter :: rhoh2o = 1000.0
! real(kind=kind_phys), parameter :: convrad = con_jcal*1.e4/60.0
real(kind=kind_phys), parameter :: a2 = 17.2693882
real(kind=kind_phys), parameter :: a3 = 273.16
real(kind=kind_phys), parameter :: a4 = 35.86
real(kind=kind_phys), parameter :: a23m4 = a2*(a3-a4)
real(kind=kind_phys), save :: zsoil_noah(4)
data zsoil_noah / -0.1, -0.4, -1.0, -2.0 /
! --- input:
integer, intent(in) :: im, km
integer, dimension(im), intent(in) :: soiltyp, vegtype, slopetyp
real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, &
& t1, q1, sigmaf, sfcemis, dlwflx, dswsfc, snet, tg3, cm, &
& ch, prsl1, prslki, ddvel, shdmin, shdmax, &
& snoalb, sfalb, zf
integer, dimension(im), intent(in) :: islimsk
real (kind=kind_phys), intent(in) :: delt
logical, dimension(im), intent(in) :: flag_iter, flag_guess
! --- in/out:
real (kind=kind_phys), dimension(im), intent(inout) :: weasd, &
& snwdph, tskin, tprcp, srflag, canopy, trans, tsurf
real (kind=kind_phys), dimension(im,km), intent(inout) :: &
& smc, stc, slc
! --- output:
real (kind=kind_phys), dimension(im), intent(out) :: sncovr1, &
& qsurf, gflux, drain, evap, hflx, ep, runoff, cmm, chh, &
& evbs, evcw, sbsno, snowc, stm, snohf, smcwlt2, smcref2, &
& zorl, wet1
! --- locals:
real (kind=kind_phys), dimension(im) :: rch, rho, &
& q0, qs1, theta1, tv1, wind, weasd_old, snwdph_old, &
& tprcp_old, srflag_old, tskin_old, canopy_old
real (kind=kind_phys), dimension(km) :: et, sldpth, stsoil, &
& smsoil, slsoil
real (kind=kind_phys), dimension(im,km) :: zsoil, smc_old, &
& stc_old, slc_old
real (kind=kind_phys) :: alb, albedo, beta, chx, cmx, cmc, &
& dew, drip, dqsdt2, ec, edir, ett, eta, esnow, etp, &
& flx1, flx2, flx3, ffrozp, lwdn, pc, prcp, ptu, q2, &
& q2sat, solnet, rc, rcs, rct, rcq, rcsoil, rsmin, &
& runoff1, runoff2, runoff3, sfcspd, sfcprs, sfctmp, &
& sfcems, sheat, shdfac, shdmin1d, shdmax1d, smcwlt, &
& smcdry, smcref, smcmax, sneqv, snoalb1d, snowh, &
& snomlt, sncovr, soilw, soilm, ssoil, tsea, th2, tbot, &
& xlai, zlvl, swdn, tem,z0
integer :: couple, ice, nsoil, nroot, slope, stype, vtype
integer :: i, k
logical :: flag(im)
!
!===> ... begin here
!
! --- ... set flag for land points
do i = 1, im
flag(i) = (islimsk(i) == 1)
enddo
! --- ... save land-related prognostic fields for guess run
do i = 1, im
if (flag(i) .and. flag_guess(i)) then
weasd_old(i) = weasd(i)
snwdph_old(i) = snwdph(i)
tskin_old(i) = tskin(i)
canopy_old(i) = canopy(i)
tprcp_old(i) = tprcp(i)
srflag_old(i) = srflag(i)
do k = 1, km
smc_old(i,k) = smc(i,k)
stc_old(i,k) = stc(i,k)
slc_old(i,k) = slc(i,k)
enddo
endif
enddo
! --- ... initialization block
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
ep(i) = 0.0
evap (i) = 0.0
hflx (i) = 0.0
gflux(i) = 0.0
drain(i) = 0.0
canopy(i) = max(canopy(i), 0.0)
evbs (i) = 0.0
evcw (i) = 0.0
trans(i) = 0.0
sbsno(i) = 0.0
snowc(i) = 0.0
snohf(i) = 0.0
endif
enddo
! --- ... initialize variables
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
wind(i) = sqrt( u1(i)*u1(i) + v1(i)*v1(i) ) &
& + max(0.0, min(ddvel(i), 30.0))
wind(i) = max(wind(i), 1.0)
q0(i) = max(q1(i), 1.e-8) !* q1=specific humidity at level 1 (kg/kg)
theta1(i) = t1(i) * prslki(i) !* adiabatic temp at level 1 (k)
tv1(i) = t1(i) * (1.0 + con_fvirt*q0(i))
rho(i) = prsl1(i) / (con_rd * tv1(i))
qs1(i) = fpvs( t1(i) ) !* qs1=sat. humidity at level 1 (kg/kg)
qs1(i) = con_eps*qs1(i) / (prsl1(i) + con_epsm1*qs1(i))
qs1(i) = max(qs1(i), 1.e-8)
q0 (i) = min(qs1(i), q0(i))
endif
enddo
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
do k = 1, km
zsoil(i,k) = zsoil_noah(k)
enddo
endif
enddo
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
! --- ... noah: prepare variables to run noah lsm
! 1. configuration information (c):
! ------------------------------
! couple - couple-uncouple flag (=1: coupled, =0: uncoupled)
! ffrozp - flag for snow-rain detection (1.=snow, 0.=rain)
! ice - sea-ice flag (=1: sea-ice, =0: land)
! dt - timestep (sec) (dt should not exceed 3600 secs) = delt
! zlvl - height (m) above ground of atmospheric forcing variables
! nsoil - number of soil layers (at least 2)
! sldpth - the thickness of each soil layer (m)
couple = 1 ! run noah lsm in 'couple' mode
if (srflag(i) == 1.0) then ! snow phase
ffrozp = 1.0
elseif (srflag(i) == 0.0) then ! rain phase
ffrozp = 0.0
endif
ice = 0
zlvl = zf(i)
nsoil = km
sldpth(1) = - zsoil(i,1)
do k = 2, km
sldpth(k) = zsoil(i,k-1) - zsoil(i,k)
enddo
! 2. forcing data (f):
! -----------------
! lwdn - lw dw radiation flux (w/m2)
! solnet - net sw radiation flux (dn-up) (w/m2)
! sfcprs - pressure at height zlvl above ground (pascals)
! prcp - precip rate (kg m-2 s-1)
! sfctmp - air temperature (k) at height zlvl above ground
! th2 - air potential temperature (k) at height zlvl above ground
! q2 - mixing ratio at height zlvl above ground (kg kg-1)
lwdn = dlwflx(i) !..downward lw flux at sfc in w/m2
swdn = dswsfc(i) !..downward sw flux at sfc in w/m2
solnet = snet(i) !..net sw rad flx (dn-up) at sfc in w/m2
sfcems = sfcemis(i)
sfcprs = prsl1(i)
prcp = rhoh2o * tprcp(i) / delt
sfctmp = t1(i)
th2 = theta1(i)
q2 = q0(i)
! 3. other forcing (input) data (i):
! ------------------------------
! sfcspd - wind speed (m s-1) at height zlvl above ground
! q2sat - sat mixing ratio at height zlvl above ground (kg kg-1)
! dqsdt2 - slope of sat specific humidity curve at t=sfctmp (kg kg-1 k-1)
sfcspd = wind(i)
q2sat = qs1(i)
dqsdt2 = q2sat * a23m4/(sfctmp-a4)**2
! 4. canopy/soil characteristics (s):
! --------------------------------
! vegtyp - vegetation type (integer index) -> vtype
! soiltyp - soil type (integer index) -> stype
! slopetyp- class of sfc slope (integer index) -> slope
! shdfac - areal fractional coverage of green vegetation (0.0-1.0)
! shdmin - minimum areal fractional coverage of green vegetation -> shdmin1d
! ptu - photo thermal unit (plant phenology for annuals/crops)
! alb - backround snow-free surface albedo (fraction)
! snoalb - upper bound on maximum albedo over deep snow -> snoalb1d
! tbot - bottom soil temperature (local yearly-mean sfc air temp)
vtype = vegtype(i)
stype = soiltyp(i)
slope = slopetyp(i)
shdfac= sigmaf(i)
shdmin1d = shdmin(i)
shdmax1d = shdmax(i)
snoalb1d = snoalb(i)
ptu = 0.0
alb = sfalb(i)
tbot = tg3(i)
! 5. history (state) variables (h):
! ------------------------------
! cmc - canopy moisture content (m)
! t1 - ground/canopy/snowpack) effective skin temperature (k) -> tsea
! stc(nsoil) - soil temp (k) -> stsoil
! smc(nsoil) - total soil moisture content (volumetric fraction) -> smsoil
! sh2o(nsoil)- unfrozen soil moisture content (volumetric fraction) -> slsoil
! snowh - actual snow depth (m)
! sneqv - liquid water-equivalent snow depth (m)
! albedo - surface albedo including snow effect (unitless fraction)
! ch - surface exchange coefficient for heat and moisture (m s-1) -> chx
! cm - surface exchange coefficient for momentum (m s-1) -> cmx
cmc = canopy(i)/1000. ! convert from mm to m
tsea = tsurf(i) ! clu_q2m_iter
do k = 1, km
stsoil(k) = stc(i,k)
smsoil(k) = smc(i,k)
slsoil(k) = slc(i,k)
enddo
snowh = snwdph(i) * 0.001 ! convert from mm to m
sneqv = weasd(i) * 0.001 ! convert from mm to m
if (sneqv /= 0.0 .and. snowh == 0.0) then
snowh = 10.0 * sneqv
endif
chx = ch(i) * wind(i) ! compute conductance
cmx = cm(i) * wind(i)
chh(i) = chx * rho(i)
cmm(i) = cmx
! --- ... call noah lsm
call sflx &
! --- inputs:
& ( nsoil, couple, ice, ffrozp, delt, zlvl, sldpth, &
& swdn, solnet, lwdn, sfcems, sfcprs, sfctmp, &
& sfcspd, prcp, q2, q2sat, dqsdt2, th2, &
& vtype, stype, slope, shdmin1d, alb, snoalb1d, &
! --- input/outputs:
& tbot, cmc, tsea, stsoil, smsoil, slsoil, sneqv, chx, cmx, &
! --- outputs:
& 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, &
& z0 )
! --- ... noah: prepare variables for return to parent mode
! 6. output (o):
! -----------
! eta - actual latent heat flux (w m-2: positive, if upward from sfc)
! sheat - sensible heat flux (w m-2: positive, if upward from sfc)
! beta - ratio of actual/potential evap (dimensionless)
! etp - potential evaporation (w m-2)
! ssoil - soil heat flux (w m-2: negative if downward from surface)
! runoff1 - surface runoff (m s-1), not infiltrating the surface
! runoff2 - subsurface runoff (m s-1), drainage out bottom
evap(i) = eta
hflx(i) = sheat
gflux(i) = ssoil
evbs(i) = edir
evcw(i) = ec
trans(i) = ett
sbsno(i) = esnow
snowc(i) = sncovr
stm(i) = soilm
snohf(i) = flx1 + flx2 + flx3
smcwlt2(i) = smcwlt
smcref2(i) = smcref
ep(i) = etp
tsurf(i) = tsea
do k = 1, km
stc(i,k) = stsoil(k)
smc(i,k) = smsoil(k)
slc(i,k) = slsoil(k)
enddo
wet1(i) = smsoil(1) / smcmax !Sarah Lu added 09/09/2010 (for GOCART)
! --- ... unit conversion (from m s-1 to mm s-1)
runoff(i) = runoff1 * 1000.0
drain (i) = runoff2 * 1000.0
! --- ... unit conversion (from m to mm)
canopy(i) = cmc * 1000.0
snwdph(i) = snowh * 1000.0
weasd(i) = sneqv * 1000.0
sncovr1(i) = sncovr
! ---- ... outside sflx, roughness uses cm as unit
zorl(i) = z0*100.
! --- ... do not return the following output fields to parent model
! ec - canopy water evaporation (m s-1)
! edir - direct soil evaporation (m s-1)
! et(nsoil)-plant transpiration from a particular root layer (m s-1)
! ett - total plant transpiration (m s-1)
! esnow - sublimation from (or deposition to if <0) snowpack (m s-1)
! drip - through-fall of precip and/or dew in excess of canopy
! water-holding capacity (m)
! dew - dewfall (or frostfall for t<273.15) (m)
! beta - ratio of actual/potential evap (dimensionless)
! flx1 - precip-snow sfc (w m-2)
! flx2 - freezing rain latent heat flux (w m-2)
! flx3 - phase-change heat flux from snowmelt (w m-2)
! snomlt - snow melt (m) (water equivalent)
! sncovr - fractional snow cover (unitless fraction, 0-1)
! runoff3 - numerical trunctation in excess of porosity (smcmax)
! for a given soil layer at the end of a time step
! rc - canopy resistance (s m-1)
! pc - plant coefficient (unitless fraction, 0-1) where pc*etp
! = actual transp
! xlai - leaf area index (dimensionless)
! rsmin - minimum canopy resistance (s m-1)
! rcs - incoming solar rc factor (dimensionless)
! rct - air temperature rc factor (dimensionless)
! rcq - atmos vapor pressure deficit rc factor (dimensionless)
! rcsoil - soil moisture rc factor (dimensionless)
! soilw - available soil moisture in root zone (unitless fraction
! between smcwlt and smcmax)
! soilm - total soil column moisture content (frozen+unfrozen) (m)
! smcwlt - wilting point (volumetric)
! smcdry - dry soil moisture threshold where direct evap frm top
! layer ends (volumetric)
! smcref - soil moisture threshold where transpiration begins to
! stress (volumetric)
! smcmax - porosity, i.e. saturated value of soil moisture
! (volumetric)
! nroot - number of root layers, a function of veg type, determined
! in subroutine redprm.
endif ! end if_flag_iter_and_flag_block
enddo ! end do_i_loop
! --- ... compute qsurf (specific humidity at sfc)
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
rch(i) = rho(i) * con_cp * ch(i) * wind(i)
qsurf(i) = q1(i) + evap(i) / (elocp * rch(i))
endif
enddo
do i = 1, im
if (flag_iter(i) .and. flag(i)) then
tem = 1.0 / rho(i)
hflx(i) = hflx(i) * tem * cpinv
evap(i) = evap(i) * tem * hvapi
endif
enddo
! --- ... restore land-related prognostic fields for guess run
do i = 1, im
if (flag(i)) then
if (flag_guess(i)) then
weasd(i) = weasd_old(i)
snwdph(i) = snwdph_old(i)
tskin(i) = tskin_old(i)
canopy(i) = canopy_old(i)
tprcp(i) = tprcp_old(i)
srflag(i) = srflag_old(i)
do k = 1, km
smc(i,k) = smc_old(i,k)
stc(i,k) = stc_old(i,k)
slc(i,k) = slc_old(i,k)
enddo
else
tskin(i) = tsurf(i)
endif
endif
enddo
!
return
!...................................
end subroutine sfc_drv
!-----------------------------------