#include <misc.h>
#include <preproc.h>
module Biogeophysics1Mod 1,1
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: Biogeophysics1Mod
!
! !DESCRIPTION:
! Performs calculation of leaf temperature and surface fluxes.
! Biogeophysics2.F90 then determines soil/snow and ground
! temperatures and updates the surface fluxes for the new ground
! temperature.
!
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
!
! !PUBLIC TYPES:
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
public :: Biogeophysics1 ! Calculate leaf temperature and surface fluxes
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: Biogeophysics1
!
! !INTERFACE:
subroutine Biogeophysics1(lbg, ubg, lbc, ubc, lbp, ubp, & 1,7
num_nolakec, filter_nolakec, num_nolakep, filter_nolakep)
!
! !DESCRIPTION:
! This is the main subroutine to execute the calculation of leaf temperature
! and surface fluxes. Biogeophysics2.F90 then determines soil/snow and ground
! temperatures and updates the surface fluxes for the new ground
! temperature.
!
! Calling sequence is:
! Biogeophysics1: surface biogeophysics driver
! -> QSat: saturated vapor pressure, specific humidity, and
! derivatives at ground surface and derivatives at
! leaf surface using updated leaf temperature
! Leaf temperature
! Foliage energy conservation is given by the foliage energy budget
! equation:
! Rnet - Hf - LEf = 0
! The equation is solved by Newton-Raphson iteration, in which this
! iteration includes the calculation of the photosynthesis and
! stomatal resistance, and the integration of turbulent flux profiles.
! The sensible and latent heat transfer between foliage and atmosphere
! and ground is linked by the equations:
! Ha = Hf + Hg and Ea = Ef + Eg
!
! !USES:
use clmtype
use clm_atmlnd , only : clm_a2l
use clm_varcon , only : denh2o, denice, roverg, hvap, hsub, &
istice, istice_mec, istwet, istsoil, isturb, istdlak, &
zlnd, zsno, tfrz, &
icol_roof, icol_sunwall, icol_shadewall, &
icol_road_imperv, icol_road_perv, tfrz, spval, istdlak
use clm_varpar , only : nlevgrnd, nlevurb, nlevsno, max_pft_per_gcell, nlevsoi
use QSatMod , only : QSat
use shr_const_mod , only : SHR_CONST_PI
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbg, ubg ! gridcell-index bounds
integer, intent(in) :: lbc, ubc ! column-index bounds
integer, intent(in) :: lbp, ubp ! pft-index bounds
integer, intent(in) :: num_nolakec ! number of column non-lake points in column filter
integer, intent(in) :: filter_nolakec(ubc-lbc+1) ! column filter for non-lake points
integer, intent(in) :: num_nolakep ! number of column non-lake points in pft filter
integer, intent(in) :: filter_nolakep(ubp-lbp+1) ! pft filter for non-lake points
!
! !CALLED FROM:
! subroutine clm_driver1
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision
! Migrated to clm2.0 by Keith Oleson and Mariana Vertenstein
! Migrated to clm2.1 new data structures by Peter Thornton and M. Vertenstein
! 27 February 2008: Keith Oleson; weighted soil/snow emissivity
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
integer , pointer :: ivt(:) !pft vegetation type
integer , pointer :: ityplun(:) !landunit type
integer , pointer :: clandunit(:) !column's landunit index
integer , pointer :: cgridcell(:) !column's gridcell index
real(r8), pointer :: pwtgcell(:) !weight relative to gridcell for each pft
integer , pointer :: ctype(:) !column type
real(r8), pointer :: forc_pbot(:) !atmospheric pressure (Pa)
real(r8), pointer :: forc_q(:) !atmospheric specific humidity (kg/kg)
real(r8), pointer :: forc_t(:) !atmospheric temperature (Kelvin)
real(r8), pointer :: forc_hgt_t(:) !observational height of temperature [m]
real(r8), pointer :: forc_hgt_u(:) !observational height of wind [m]
real(r8), pointer :: forc_hgt_q(:) !observational height of specific humidity [m]
integer , pointer :: npfts(:) !number of pfts on gridcell
integer , pointer :: pfti(:) !initial pft on gridcell
integer , pointer :: plandunit(:) !pft's landunit index
real(r8), pointer :: forc_hgt_u_pft(:) !observational height of wind at pft level [m]
real(r8), pointer :: forc_hgt_t_pft(:) !observational height of temperature at pft level [m]
real(r8), pointer :: forc_hgt_q_pft(:) !observational height of specific humidity at pft level [m]
integer , pointer :: frac_veg_nosno(:) !fraction of vegetation not covered by snow (0 OR 1) [-]
integer , pointer :: pgridcell(:) !pft's gridcell index
integer , pointer :: pcolumn(:) !pft's column index
real(r8), pointer :: z_0_town(:) !momentum roughness length of urban landunit (m)
real(r8), pointer :: z_d_town(:) !displacement height of urban landunit (m)
real(r8), pointer :: forc_th(:) !atmospheric potential temperature (Kelvin)
real(r8), pointer :: forc_u(:) !atmospheric wind speed in east direction (m/s)
real(r8), pointer :: forc_v(:) !atmospheric wind speed in north direction (m/s)
real(r8), pointer :: smpmin(:) !restriction for min of soil potential (mm)
integer , pointer :: snl(:) !number of snow layers
real(r8), pointer :: frac_sno(:) !fraction of ground covered by snow (0 to 1)
real(r8), pointer :: h2osno(:) !snow water (mm H2O)
real(r8), pointer :: elai(:) !one-sided leaf area index with burying by snow
real(r8), pointer :: esai(:) !one-sided stem area index with burying by snow
real(r8), pointer :: z0mr(:) !ratio of momentum roughness length to canopy top height (-)
real(r8), pointer :: displar(:) !ratio of displacement height to canopy top height (-)
real(r8), pointer :: htop(:) !canopy top (m)
real(r8), pointer :: dz(:,:) !layer depth (m)
real(r8), pointer :: t_soisno(:,:) !soil temperature (Kelvin)
real(r8), pointer :: h2osoi_liq(:,:) !liquid water (kg/m2)
real(r8), pointer :: h2osoi_ice(:,:) !ice lens (kg/m2)
real(r8), pointer :: watsat(:,:) !volumetric soil water at saturation (porosity)
real(r8), pointer :: sucsat(:,:) !minimum soil suction (mm)
real(r8), pointer :: bsw(:,:) !Clapp and Hornberger "b"
real(r8), pointer :: watfc(:,:) !volumetric soil water at field capacity
real(r8), pointer :: watopt(:,:) !volumetric soil moisture corresponding to no restriction on ET from urban pervious surface
real(r8), pointer :: watdry(:,:) !volumetric soil moisture corresponding to no restriction on ET from urban pervious surface
real(r8), pointer :: rootfr_road_perv(:,:) !fraction of roots in each soil layer for urban pervious road
real(r8), pointer :: rootr_road_perv(:,:) !effective fraction of roots in each soil layer for urban pervious road
!
! local pointers to implicit out arguments
!
real(r8), pointer :: t_grnd(:) !ground temperature (Kelvin)
real(r8), pointer :: qg(:) !ground specific humidity [kg/kg]
real(r8), pointer :: dqgdT(:) !d(qg)/dT
real(r8), pointer :: emg(:) !ground emissivity
real(r8), pointer :: htvp(:) !latent heat of vapor of water (or sublimation) [j/kg]
real(r8), pointer :: beta(:) !coefficient of convective velocity [-]
real(r8), pointer :: zii(:) !convective boundary height [m]
real(r8), pointer :: thm(:) !intermediate variable (forc_t+0.0098*forc_hgt_t_pft)
real(r8), pointer :: thv(:) !virtual potential temperature (kelvin)
real(r8), pointer :: z0mg(:) !roughness length over ground, momentum [m]
real(r8), pointer :: z0hg(:) !roughness length over ground, sensible heat [m]
real(r8), pointer :: z0qg(:) !roughness length over ground, latent heat [m]
real(r8), pointer :: emv(:) !vegetation emissivity
real(r8), pointer :: z0m(:) !momentum roughness length (m)
real(r8), pointer :: displa(:) !displacement height (m)
real(r8), pointer :: z0mv(:) !roughness length over vegetation, momentum [m]
real(r8), pointer :: z0hv(:) !roughness length over vegetation, sensible heat [m]
real(r8), pointer :: z0qv(:) !roughness length over vegetation, latent heat [m]
real(r8), pointer :: eflx_sh_tot(:) !total sensible heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_sh_tot_u(:) !urban total sensible heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_sh_tot_r(:) !rural total sensible heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_lh_tot(:) !total latent heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_lh_tot_u(:) !urban total latent heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_lh_tot_r(:) !rural total latent heat flux (W/m**2) [+ to atm]
real(r8), pointer :: eflx_sh_veg(:) !sensible heat flux from leaves (W/m**2) [+ to atm]
real(r8), pointer :: qflx_evap_tot(:) !qflx_evap_soi + qflx_evap_veg + qflx_tran_veg
real(r8), pointer :: qflx_evap_veg(:) !vegetation evaporation (mm H2O/s) (+ = to atm)
real(r8), pointer :: qflx_tran_veg(:) !vegetation transpiration (mm H2O/s) (+ = to atm)
real(r8), pointer :: cgrnd(:) !deriv. of soil energy flux wrt to soil temp [w/m2/k]
real(r8), pointer :: cgrnds(:) !deriv. of soil sensible heat flux wrt soil temp [w/m2/k]
real(r8), pointer :: cgrndl(:) !deriv. of soil latent heat flux wrt soil temp [w/m**2/k]
real(r8) ,pointer :: tssbef(:,:) !soil/snow temperature before update
real(r8) ,pointer :: soilalpha(:) !factor that reduces ground saturated specific humidity (-)
real(r8) ,pointer :: soilbeta(:) !factor that reduces ground evaporation
real(r8) ,pointer :: soilalpha_u(:) !Urban factor that reduces ground saturated specific humidity (-)
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
integer :: g,l,c,p !indices
integer :: j !soil/snow level index
integer :: fp !lake filter pft index
integer :: fc !lake filter column index
real(r8) :: qred !soil surface relative humidity
real(r8) :: avmuir !ir inverse optical depth per unit leaf area
real(r8) :: eg !water vapor pressure at temperature T [pa]
real(r8) :: qsatg !saturated humidity [kg/kg]
real(r8) :: degdT !d(eg)/dT
real(r8) :: qsatgdT !d(qsatg)/dT
real(r8) :: fac !soil wetness of surface layer
real(r8) :: psit !negative potential of soil
real(r8) :: hr !relative humidity
real(r8) :: hr_road_perv !relative humidity for urban pervious road
real(r8) :: wx !partial volume of ice and water of surface layer
real(r8) :: fac_fc !soil wetness of surface layer relative to field capacity
real(r8) :: eff_porosity ! effective porosity in layer
real(r8) :: vol_ice ! partial volume of ice lens in layer
real(r8) :: vol_liq ! partial volume of liquid water in layer
integer :: pi !index
!------------------------------------------------------------------------------
! Assign local pointers to derived type members (gridcell-level)
forc_hgt_t => clm_a2l%forc_hgt_t
forc_u => clm_a2l%forc_u
forc_v => clm_a2l%forc_v
forc_hgt_u => clm_a2l%forc_hgt_u
forc_hgt_q => clm_a2l%forc_hgt_q
npfts => clm3%g%npfts
pfti => clm3%g%pfti
! Assign local pointers to derived type members (landunit-level)
ityplun => clm3%g%l%itype
z_0_town => clm3%g%l%z_0_town
z_d_town => clm3%g%l%z_d_town
! Assign local pointers to derived type members (column-level)
forc_pbot => clm3%g%l%c%cps%forc_pbot
forc_q => clm3%g%l%c%cws%forc_q
forc_t => clm3%g%l%c%ces%forc_t
forc_th => clm3%g%l%c%ces%forc_th
cgridcell => clm3%g%l%c%gridcell
clandunit => clm3%g%l%c%landunit
ctype => clm3%g%l%c%itype
beta => clm3%g%l%c%cps%beta
dqgdT => clm3%g%l%c%cws%dqgdT
emg => clm3%g%l%c%cps%emg
frac_sno => clm3%g%l%c%cps%frac_sno
h2osno => clm3%g%l%c%cws%h2osno
htvp => clm3%g%l%c%cps%htvp
qg => clm3%g%l%c%cws%qg
smpmin => clm3%g%l%c%cps%smpmin
snl => clm3%g%l%c%cps%snl
t_grnd => clm3%g%l%c%ces%t_grnd
thv => clm3%g%l%c%ces%thv
z0hg => clm3%g%l%c%cps%z0hg
z0mg => clm3%g%l%c%cps%z0mg
z0qg => clm3%g%l%c%cps%z0qg
zii => clm3%g%l%c%cps%zii
bsw => clm3%g%l%c%cps%bsw
dz => clm3%g%l%c%cps%dz
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
soilalpha => clm3%g%l%c%cws%soilalpha
soilbeta => clm3%g%l%c%cws%soilbeta
soilalpha_u => clm3%g%l%c%cws%soilalpha_u
sucsat => clm3%g%l%c%cps%sucsat
t_soisno => clm3%g%l%c%ces%t_soisno
tssbef => clm3%g%l%c%ces%tssbef
watsat => clm3%g%l%c%cps%watsat
watfc => clm3%g%l%c%cps%watfc
watdry => clm3%g%l%c%cps%watdry
watopt => clm3%g%l%c%cps%watopt
rootfr_road_perv => clm3%g%l%c%cps%rootfr_road_perv
rootr_road_perv => clm3%g%l%c%cps%rootr_road_perv
! Assign local pointers to derived type members (pft-level)
ivt => clm3%g%l%c%p%itype
elai => clm3%g%l%c%p%pps%elai
esai => clm3%g%l%c%p%pps%esai
htop => clm3%g%l%c%p%pps%htop
emv => clm3%g%l%c%p%pps%emv
z0m => clm3%g%l%c%p%pps%z0m
displa => clm3%g%l%c%p%pps%displa
z0mv => clm3%g%l%c%p%pps%z0mv
z0hv => clm3%g%l%c%p%pps%z0hv
z0qv => clm3%g%l%c%p%pps%z0qv
eflx_sh_tot => clm3%g%l%c%p%pef%eflx_sh_tot
eflx_sh_tot_u => clm3%g%l%c%p%pef%eflx_sh_tot_u
eflx_sh_tot_r => clm3%g%l%c%p%pef%eflx_sh_tot_r
eflx_lh_tot => clm3%g%l%c%p%pef%eflx_lh_tot
eflx_lh_tot_u => clm3%g%l%c%p%pef%eflx_lh_tot_u
eflx_lh_tot_r => clm3%g%l%c%p%pef%eflx_lh_tot_r
eflx_sh_veg => clm3%g%l%c%p%pef%eflx_sh_veg
qflx_evap_tot => clm3%g%l%c%p%pwf%qflx_evap_tot
qflx_evap_veg => clm3%g%l%c%p%pwf%qflx_evap_veg
qflx_tran_veg => clm3%g%l%c%p%pwf%qflx_tran_veg
cgrnd => clm3%g%l%c%p%pef%cgrnd
cgrnds => clm3%g%l%c%p%pef%cgrnds
cgrndl => clm3%g%l%c%p%pef%cgrndl
forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft
forc_hgt_t_pft => clm3%g%l%c%p%pps%forc_hgt_t_pft
forc_hgt_q_pft => clm3%g%l%c%p%pps%forc_hgt_q_pft
plandunit => clm3%g%l%c%p%landunit
frac_veg_nosno => clm3%g%l%c%p%pps%frac_veg_nosno
thm => clm3%g%l%c%p%pes%thm
pgridcell => clm3%g%l%c%p%gridcell
pcolumn => clm3%g%l%c%p%column
pwtgcell => clm3%g%l%c%p%wtgcell
! Assign local pointers to derived type members (ecophysiological)
z0mr => pftcon%z0mr
displar => pftcon%displar
do j = -nlevsno+1, nlevgrnd
do fc = 1,num_nolakec
c = filter_nolakec(fc)
tssbef(c,j) = t_soisno(c,j)
end do
end do
do fc = 1,num_nolakec
c = filter_nolakec(fc)
l = clandunit(c)
if (ctype(c) == icol_road_perv) then
hr_road_perv = 0._r8
end if
! begin calculations that relate only to the column level
! Ground and soil temperatures from previous time step
t_grnd(c) = t_soisno(c,snl(c)+1)
! Saturated vapor pressure, specific humidity and their derivatives
! at ground surface
qred = 1._r8
if (ityplun(l)/=istwet .AND. ityplun(l)/=istice &
.AND. ityplun(l)/=istice_mec) then
if (ityplun(l) == istsoil) then
wx = (h2osoi_liq(c,1)/denh2o+h2osoi_ice(c,1)/denice)/dz(c,1)
fac = min(1._r8, wx/watsat(c,1))
fac = max( fac, 0.01_r8 )
psit = -sucsat(c,1) * fac ** (-bsw(c,1))
psit = max(smpmin(c), psit)
hr = exp(psit/roverg/t_grnd(c))
qred = (1.-frac_sno(c))*hr + frac_sno(c)
!! Lee and Pielke 1992 beta, added by K.Sakaguchi
if (wx < watfc(c,1) ) then !when water content of ths top layer is less than that at F.C.
fac_fc = min(1._r8, wx/watfc(c,1)) !eqn5.66 but divided by theta at field capacity
fac_fc = max( fac_fc, 0.01_r8 )
! modifiy soil beta by snow cover. soilbeta for snow surface is one
soilbeta(c) = (1._r8-frac_sno(c))*0.25_r8*(1._r8 - cos(SHR_CONST_PI*fac_fc))**2._r8 &
+ frac_sno(c)
else !when water content of ths top layer is more than that at F.C.
soilbeta(c) = 1._r8
end if
soilalpha(c) = qred
! Pervious road depends on water in total soil column
else if (ctype(c) == icol_road_perv) then
do j = 1, nlevsoi
if (t_soisno(c,j) >= tfrz) then
vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
eff_porosity = watsat(c,j)-vol_ice
vol_liq = min(eff_porosity, h2osoi_liq(c,j)/(dz(c,j)*denh2o))
fac = min( max(vol_liq-watdry(c,j),0._r8) / (watopt(c,j)-watdry(c,j)), 1._r8 )
else
fac = 0._r8
end if
rootr_road_perv(c,j) = rootfr_road_perv(c,j)*fac
hr_road_perv = hr_road_perv + rootr_road_perv(c,j)
end do
! Allows for sublimation of snow or dew on snow
qred = (1.-frac_sno(c))*hr_road_perv + frac_sno(c)
! Normalize root resistances to get layer contribution to total ET
if (hr_road_perv .gt. 0._r8) then
do j = 1, nlevsoi
rootr_road_perv(c,j) = rootr_road_perv(c,j)/hr_road_perv
end do
end if
soilalpha_u(c) = qred
soilbeta(c) = 0._r8
else if (ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall) then
qred = 0._r8
soilbeta(c) = 0._r8
soilalpha_u(c) = spval
else if (ctype(c) == icol_roof .or. ctype(c) == icol_road_imperv) then
qred = 1._r8
soilbeta(c) = 0._r8
soilalpha_u(c) = spval
end if
else
soilalpha(c) = spval
soilbeta(c) = 1._r8
end if
call QSat(t_grnd(c), forc_pbot(c), eg, degdT, qsatg, qsatgdT)
qg(c) = qred*qsatg
dqgdT(c) = qred*qsatgdT
if (qsatg > forc_q(c) .and. forc_q(c) > qred*qsatg) then
qg(c) = forc_q(c)
dqgdT(c) = 0._r8
end if
! Ground emissivity - only calculate for non-urban landunits
! Urban emissivities are currently read in from data file
if (ityplun(l) /= isturb) then
if (ityplun(l)==istice .or. ityplun(l)==istice_mec) then
emg(c) = 0.97_r8
else
emg(c) = (1._r8-frac_sno(c))*0.96_r8 + frac_sno(c)*0.97_r8
end if
end if
! Latent heat. We arbitrarily assume that the sublimation occurs
! only as h2osoi_liq = 0
htvp(c) = hvap
if (h2osoi_liq(c,snl(c)+1) <= 0._r8 .and. h2osoi_ice(c,snl(c)+1) > 0._r8) htvp(c) = hsub
! Switch between vaporization and sublimation causes rapid solution
! separation in perturbation growth test
#if (defined PERGRO)
htvp(c) = hvap
#endif
! Ground roughness lengths over non-lake columns (includes bare ground, ground
! underneath canopy, wetlands, etc.)
if (frac_sno(c) > 0._r8) then
z0mg(c) = zsno
else
z0mg(c) = zlnd
end if
z0hg(c) = z0mg(c) ! initial set only
z0qg(c) = z0mg(c) ! initial set only
! Potential, virtual potential temperature, and wind speed at the
! reference height
beta(c) = 1._r8
zii(c) = 1000._r8
thv(c) = forc_th(c)*(1._r8+0.61_r8*forc_q(c))
end do ! (end of columns loop)
! Initialization
do fp = 1,num_nolakep
p = filter_nolakep(fp)
! Initial set (needed for history tape fields)
eflx_sh_tot(p) = 0._r8
l = plandunit(p)
if (ityplun(l) == isturb) then
eflx_sh_tot_u(p) = 0._r8
else if (ityplun(l) == istsoil) then
eflx_sh_tot_r(p) = 0._r8
end if
eflx_lh_tot(p) = 0._r8
if (ityplun(l) == isturb) then
eflx_lh_tot_u(p) = 0._r8
else if (ityplun(l) == istsoil) then
eflx_lh_tot_r(p) = 0._r8
end if
eflx_sh_veg(p) = 0._r8
qflx_evap_tot(p) = 0._r8
qflx_evap_veg(p) = 0._r8
qflx_tran_veg(p) = 0._r8
! Initial set for calculation
cgrnd(p) = 0._r8
cgrnds(p) = 0._r8
cgrndl(p) = 0._r8
! Vegetation Emissivity
avmuir = 1._r8
emv(p) = 1._r8-exp(-(elai(p)+esai(p))/avmuir)
! Roughness lengths over vegetation
z0m(p) = z0mr(ivt(p)) * htop(p)
displa(p) = displar(ivt(p)) * htop(p)
z0mv(p) = z0m(p)
z0hv(p) = z0mv(p)
z0qv(p) = z0mv(p)
end do
! Make forcing height a pft-level quantity that is the atmospheric forcing
! height plus each pft's z0m+displa
do pi = 1,max_pft_per_gcell
do g = lbg, ubg
if (pi <= npfts(g)) then
p = pfti(g) + pi - 1
l = plandunit(p)
! Note: Some glacier_mec pfts may have zero weight
if (pwtgcell(p) > 0._r8 .or. ityplun(l)==istice_mec) then
c = pcolumn(p)
if (ityplun(l) == istsoil) then
if (frac_veg_nosno(p) == 0) then
forc_hgt_u_pft(p) = forc_hgt_u(g) + z0mg(c) + displa(p)
forc_hgt_t_pft(p) = forc_hgt_t(g) + z0mg(c) + displa(p)
forc_hgt_q_pft(p) = forc_hgt_q(g) + z0mg(c) + displa(p)
else
forc_hgt_u_pft(p) = forc_hgt_u(g) + z0m(p) + displa(p)
forc_hgt_t_pft(p) = forc_hgt_t(g) + z0m(p) + displa(p)
forc_hgt_q_pft(p) = forc_hgt_q(g) + z0m(p) + displa(p)
end if
else if (ityplun(l) == istwet .or. ityplun(l) == istice &
.or. ityplun(l) == istice_mec) then
forc_hgt_u_pft(p) = forc_hgt_u(g) + z0mg(c)
forc_hgt_t_pft(p) = forc_hgt_t(g) + z0mg(c)
forc_hgt_q_pft(p) = forc_hgt_q(g) + z0mg(c)
else if (ityplun(l) == istdlak) then
! Should change the roughness lengths to shared constants
if (t_grnd(c) >= tfrz) then
forc_hgt_u_pft(p) = forc_hgt_u(g) + 0.01_r8
forc_hgt_t_pft(p) = forc_hgt_t(g) + 0.01_r8
forc_hgt_q_pft(p) = forc_hgt_q(g) + 0.01_r8
else
forc_hgt_u_pft(p) = forc_hgt_u(g) + 0.04_r8
forc_hgt_t_pft(p) = forc_hgt_t(g) + 0.04_r8
forc_hgt_q_pft(p) = forc_hgt_q(g) + 0.04_r8
end if
else if (ityplun(l) == isturb) then
forc_hgt_u_pft(p) = forc_hgt_u(g) + z_0_town(l) + z_d_town(l)
forc_hgt_t_pft(p) = forc_hgt_t(g) + z_0_town(l) + z_d_town(l)
forc_hgt_q_pft(p) = forc_hgt_q(g) + z_0_town(l) + z_d_town(l)
end if
end if
end if
end do
end do
do fp = 1,num_nolakep
p = filter_nolakep(fp)
c = pcolumn(p)
thm(p) = forc_t(c) + 0.0098_r8*forc_hgt_t_pft(p)
end do
end subroutine Biogeophysics1
end module Biogeophysics1Mod