module UrbBuildTempOleson2015Mod
#include "shr_assert.h"
!-----------------------------------------------------------------------
! !DESCRIPTION:
! Calculates internal building air temperature
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_log_mod , only : errMsg => shr_log_errMsg
use decompMod , only : bounds_type
use abortutils , only : endrun
use perf_mod , only : t_startf, t_stopf
use clm_varctl , only : iulog
use UrbanParamsType , only : urbanparams_type
use UrbanTimeVarType , only : urbantv_type
use EnergyFluxType , only : energyflux_type
use TemperatureType , only : temperature_type
use LandunitType , only : lun
use ColumnType , only : col
!
! !PUBLIC TYPES:
implicit none
save
private
!
! !PUBLIC MEMBER FUNCTIONS:
public :: BuildingTemperature ! Calculation of interior building air temperature, inner
! surface temperatures of walls and roof, and floor temperature
character(len=*), parameter, private :: sourcefile = &
__FILE__
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: BuildingTemperature
!
! !INTERFACE:
subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec, &
filter_nolakec, tk, urbanparams_inst, temperature_inst, &
energyflux_inst, urbantv_inst)
!
! !DESCRIPTION:
! Solve for t_building, inner surface temperatures of roof, sunw, shdw, and floor temperature
! Five equations, five unknowns (t_roof_inner,t_sunw_inner,t_shdw_inner,t_floor,t_building at n+1)
! Derived from energy balance equations at each surface and building air
! rd (radiation), cd (conduction), cv (convection)
! qrd_roof + qcd_roof + qcv_roof = 0
! qrd_sunw + qcd_sunw + qcv_sunw = 0
! qrd_shdw + qcd_shdw + qcv_shdw = 0
! qrd_floor + qcd_floor + qcv_floor = 0
! Vbld*rho_dair*cpair*(dt_building/dt) = sum(Asfc*hcv_sfc*(t_sfc - t_building)
! + Vvent*rho_dair*cpair*(taf - t_building)
! where Vlbd is volume of building air,
! rho_dair is density of dry air at t_building (kg m-3),
! cpair is specific heat of dry air (J kg-1 K-1),
! dt_building is change in interior building temperature (K),
! dt is timestep (s),
! Asfc is surface area of roof, sunw, shdw, floor (m2)
! hcv_sfc is convective heat transfer coefficient for roof, sunw, shdw, floor (W m-2 K-1)
! t_sfc is inner surface temperature of roof, sunw, shdw, floor (K)
! t_building is interior building temperature (K)
! Vvent is ventilation airflow rate (m3 s-1)
! taf is urban canyon air temperature (K)
!
! This methodology was introduced as part of CLM5.0.
!
! Conduction fluxes are obtained from terms of soil temperature equations
! Radiation fluxes are obtained from linearizing the longwave radiation equations taking into
! account view factors for each surface.
! qrd is positive away from the surface toward room air, so qrd = emitted - absorbed,
! so positive qrd will result in a decrease in temperature
! qcd_floor is positive away from surface toward room air, so positive
! qcd will result in a decrease in temperature
! qcv is positive toward room air, so positive qcv (t_surface > t_room) will
! result in a decrease in temperature
! The LAPACK routine DGESV is used to compute the solution to the real system of linear equations
! a * x = b,
! where a is an n-by-n matrix and x and b are n-by-nrhs matrices.
!
! The LU decomposition with partial pivoting and row interchanges is
! used to factor a as
! a = P * L * U,
! where P is a permutation matrix, L is unit lower triangular, and U is
! upper triangular. The factored form of a is then used to solve the
! system of equations a * x = b.
! The following is from LAPACK documentation
! DGESV computes the solution to system of linear equations A * X = B for GE matrices
!
! =========== DOCUMENTATION ===========
!
! Online html documentation available at
! http://www.netlib.org/lapack/explore-html/
!
! Download DGESV + dependencies
!
! [TGZ]
!
! [ZIP]
!
! [TXT]
!
! Definition:
! ===========
!
! SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
!
! .. Scalar Arguments ..
! INTEGER INFO, LDA, LDB, N, NRHS
! ..
! .. Array Arguments ..
! INTEGER IPIV( * )
! DOUBLE PRECISION A( LDA, * ), B( LDB, * )
! ..
!
!
! =============
!
!
! DGESV computes the solution to a real system of linear equations
! A * X = B,
! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!
! The LU decomposition with partial pivoting and row interchanges is
! used to factor A as
! A = P * L * U,
! where P is a permutation matrix, L is unit lower triangular, and U is
! upper triangular. The factored form of A is then used to solve the
! system of equations A * X = B.
!
! Arguments:
! ==========
!
! \param[in] N
! N is INTEGER
! The number of linear equations, i.e., the order of the
! matrix A. N >= 0.
!
! \param[in] NRHS
! NRHS is INTEGER
! The number of right hand sides, i.e., the number of columns
! of the matrix B. NRHS >= 0.
!
! \param[in,out] A
! A is DOUBLE PRECISION array, dimension (LDA,N)
! On entry, the N-by-N coefficient matrix A.
! On exit, the factors L and U from the factorization
! A = P*L*U; the unit diagonal elements of L are not stored.
!
! \param[in] LDA
! LDA is INTEGER
! The leading dimension of the array A. LDA >= max(1,N).
!
! \param[out] IPIV
! IPIV is INTEGER array, dimension (N)
! The pivot indices that define the permutation matrix P;
! row i of the matrix was interchanged with row IPIV(i).
!
! \param[in,out] B
! B is DOUBLE PRECISION array, dimension (LDB,NRHS)
! On entry, the N-by-NRHS matrix of right hand side matrix B.
! On exit, if INFO = 0, the N-by-NRHS solution matrix X.
!
! \param[in] LDB
! LDB is INTEGER
! The leading dimension of the array B. LDB >= max(1,N).
!
! \param[out] INFO
! INFO is INTEGER
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
! > 0: if INFO = i, U(i,i) is exactly zero. The factorization
! has been completed, but the factor U is exactly
! singular, so the solution could not be computed.
!
! Authors:
! ========
!
! \author Univ. of Tennessee
! \author Univ. of California Berkeley
! \author Univ. of Colorado Denver
! \author NAG Ltd.
!
! \date November 2011
!
! \ingroup doubleGEsolve
! !CALLED FROM:
! subroutine SoilTemperature in this module
!
! !REVISION HISTORY:
! 08/17/12 Keith Oleson: Initial code
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use clm_time_manager, only : get_step_size
use clm_varcon , only : rair, pstd, cpair, sb, hcv_roof, hcv_roof_enhanced, &
hcv_floor, hcv_floor_enhanced, hcv_sunw, hcv_shdw, &
em_roof_int, em_floor_int, em_sunw_int, em_shdw_int, &
dz_floor, dens_floor, cp_floor, vent_ach
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall
use clm_varctl , only : iulog
use abortutils , only : endrun
use clm_varpar , only : nlevurb, nlevsno, nlevgrnd
use UrbanParamsType , only : urban_hac, urban_hac_off, urban_hac_on, urban_wasteheat_on
!
! !ARGUMENTS:
implicit none
type(bounds_type), intent(in) :: bounds ! bounds
integer , intent(in) :: num_nolakec ! number of column non-lake points in column filter
integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points
integer , intent(in) :: num_urbanl ! number of urban landunits in clump
integer , intent(in) :: filter_urbanl(:) ! urban landunit filter
real(r8), intent(in) :: tk(bounds%begc: , -nlevsno+1: ) ! thermal conductivity (W m-1 K-1) [col, j]
type(urbanparams_type), intent(in) :: urbanparams_inst ! urban parameters
type(temperature_type), intent(inout) :: temperature_inst ! temperature variables
type(energyflux_type) , intent(inout) :: energyflux_inst ! energy flux variables
type(urbantv_type) , intent(in) :: urbantv_inst ! urban time varying variables
!
! !LOCAL VARIABLES:
integer, parameter :: neq = 5 ! number of equation/unknowns
integer :: fc,fl,c,l ! indices
real(r8) :: dtime ! land model time step (s)
real(r8) :: t_roof_inner_bef(bounds%begl:bounds%endl) ! roof inside surface temperature at previous time step (K)
real(r8) :: t_sunw_inner_bef(bounds%begl:bounds%endl) ! sunwall inside surface temperature at previous time step (K)
real(r8) :: t_shdw_inner_bef(bounds%begl:bounds%endl) ! shadewall inside surface temperature at previous time step (K)
real(r8) :: t_floor_bef(bounds%begl:bounds%endl) ! floor temperature at previous time step (K)
real(r8) :: t_building_bef(bounds%begl:bounds%endl) ! internal building air temperature at previous time step [K]
real(r8) :: t_building_bef_hac(bounds%begl:bounds%endl)! internal building air temperature before applying HAC [K]
real(r8) :: hcv_roofi(bounds%begl:bounds%endl) ! roof convective heat transfer coefficient (W m-2 K-1)
real(r8) :: hcv_sunwi(bounds%begl:bounds%endl) ! sunwall convective heat transfer coefficient (W m-2 K-1)
real(r8) :: hcv_shdwi(bounds%begl:bounds%endl) ! shadewall convective heat transfer coefficient (W m-2 K-1)
real(r8) :: hcv_floori(bounds%begl:bounds%endl) ! floor convective heat transfer coefficient (W m-2 K-1)
real(r8) :: em_roofi(bounds%begl:bounds%endl) ! roof inside surface emissivity (-)
real(r8) :: em_sunwi(bounds%begl:bounds%endl) ! sunwall inside surface emissivity (-)
real(r8) :: em_shdwi(bounds%begl:bounds%endl) ! shadewall inside surface emissivity (-)
real(r8) :: em_floori(bounds%begl:bounds%endl) ! floor inside surface emissivity (-)
real(r8) :: dz_floori(bounds%begl:bounds%endl) ! concrete floor thickness (m)
real(r8) :: cp_floori(bounds%begl:bounds%endl) ! concrete floor volumetric heat capacity (J m-3 K-1)
real(r8) :: cv_floori(bounds%begl:bounds%endl) ! intermediate calculation for concrete floor (W m-2 K-1)
real(r8) :: rho_dair(bounds%begl:bounds%endl) ! density of dry air at standard pressure and t_building (kg m-3)
real(r8) :: vf_rf(bounds%begl:bounds%endl) ! view factor of roof for floor (-)
real(r8) :: vf_fr(bounds%begl:bounds%endl) ! view factor of floor for roof (-)
real(r8) :: vf_wf(bounds%begl:bounds%endl) ! view factor of wall for floor (-)
real(r8) :: vf_fw(bounds%begl:bounds%endl) ! view factor of floor for wall (-)
real(r8) :: vf_rw(bounds%begl:bounds%endl) ! view factor of roof for wall (-)
real(r8) :: vf_wr(bounds%begl:bounds%endl) ! view factor of wall for roof (-)
real(r8) :: vf_ww(bounds%begl:bounds%endl) ! view factor of wall for wall (-)
real(r8) :: zi_roof_innerl(bounds%begl:bounds%endl) ! interface depth of nlevurb roof (m)
real(r8) :: z_roof_innerl(bounds%begl:bounds%endl) ! node depth of nlevurb roof (m)
real(r8) :: zi_sunw_innerl(bounds%begl:bounds%endl) ! interface depth of nlevurb sunwall (m)
real(r8) :: z_sunw_innerl(bounds%begl:bounds%endl) ! node depth of nlevurb sunwall (m)
real(r8) :: zi_shdw_innerl(bounds%begl:bounds%endl) ! interface depth of nlevurb shadewall (m)
real(r8) :: z_shdw_innerl(bounds%begl:bounds%endl) ! node depth of nlevurb shadewall (m)
real(r8) :: t_roof_innerl_bef(bounds%begl:bounds%endl) ! roof temperature at nlevurb node depth at previous time step (K)
real(r8) :: t_sunw_innerl_bef(bounds%begl:bounds%endl) ! sunwall temperature at nlevurb node depth at previous time step (K)
real(r8) :: t_shdw_innerl_bef(bounds%begl:bounds%endl) ! shadewall temperature at nlevurb node depth at previous time step (K)
real(r8) :: t_roof_innerl(bounds%begl:bounds%endl) ! roof temperature at nlevurb node depth (K)
real(r8) :: t_sunw_innerl(bounds%begl:bounds%endl) ! sunwall temperature at nlevurb node depth (K)
real(r8) :: t_shdw_innerl(bounds%begl:bounds%endl) ! shadewall temperature at nlevurb node depth (K)
real(r8) :: tk_roof_innerl(bounds%begl:bounds%endl) ! roof thermal conductivity at nlevurb interface depth (W m-1 K-1)
real(r8) :: tk_sunw_innerl(bounds%begl:bounds%endl) ! sunwall thermal conductivity at nlevurb interface depth (W m-1 K-1)
real(r8) :: tk_shdw_innerl(bounds%begl:bounds%endl) ! shadewall thermal conductivity at nlevurb interface depth (W m-1 K-1)
real(r8) :: qrd_roof(bounds%begl:bounds%endl) ! roof inside net longwave for energy balance check (W m-2)
real(r8) :: qrd_sunw(bounds%begl:bounds%endl) ! sunwall inside net longwave for energy balance check (W m-2)
real(r8) :: qrd_shdw(bounds%begl:bounds%endl) ! shadewall inside net longwave for energy balance check (W m-2)
real(r8) :: qrd_floor(bounds%begl:bounds%endl) ! floor inside net longwave for energy balance check (W m-2)
real(r8) :: qrd_building(bounds%begl:bounds%endl) ! building inside net longwave for energy balance check (W m-2)
real(r8) :: qcv_roof(bounds%begl:bounds%endl) ! roof inside convection flux for energy balance check (W m-2)
real(r8) :: qcv_sunw(bounds%begl:bounds%endl) ! sunwall inside convection flux for energy balance check (W m-2)
real(r8) :: qcv_shdw(bounds%begl:bounds%endl) ! shadewall inside convection flux for energy balance check (W m-2)
real(r8) :: qcv_floor(bounds%begl:bounds%endl) ! floor inside convection flux for energy balance check (W m-2)
real(r8) :: qcd_roof(bounds%begl:bounds%endl) ! roof inside conduction flux for energy balance check (W m-2)
real(r8) :: qcd_sunw(bounds%begl:bounds%endl) ! sunwall inside conduction flux for energy balance check (W m-2)
real(r8) :: qcd_shdw(bounds%begl:bounds%endl) ! shadewall inside conduction flux for energy balance check (W m-2)
real(r8) :: qcd_floor(bounds%begl:bounds%endl) ! floor inside conduction flux for energy balance check (W m-2)
real(r8) :: enrgy_bal_roof(bounds%begl:bounds%endl) ! roof inside energy balance (W m-2)
real(r8) :: enrgy_bal_sunw(bounds%begl:bounds%endl) ! sunwall inside energy balance (W m-2)
real(r8) :: enrgy_bal_shdw(bounds%begl:bounds%endl) ! shadewall inside energy balance (W m-2)
real(r8) :: enrgy_bal_floor(bounds%begl:bounds%endl) ! floor inside energy balance (W m-2)
real(r8) :: enrgy_bal_buildair(bounds%begl:bounds%endl)! building air energy balance (W m-2)
real(r8) :: sum ! sum of view factors for floor, wall, roof
integer :: n ! number of linear equations (= neq)
integer :: nrhs ! number of right hand sides (= 1)
real(r8) :: a(neq,neq) ! n-by-n coefficient matrix a
integer :: lda ! leading dimension of the matrix a
integer :: ldb ! leading dimension of the matrix b
real(r8) :: result(neq) ! on entry, the right hand side of matrix b
! on exit, if info = 0, the n-by-nrhs solution matrix x
integer :: info ! exit information for LAPACK routine dgesv
integer :: ipiv(neq) ! the pivot indices that define the permutation matrix P
!EOP
!-----------------------------------------------------------------------
! Enforce expected array sizes
SHR_ASSERT_ALL((ubound(tk) == (/bounds%endc, nlevgrnd/)), errMsg(sourcefile, __LINE__))
associate(&
clandunit => col%landunit , & ! Input: [integer (:)] column's landunit
ctype => col%itype , & ! Input: [integer (:)] column type
zi => col%zi , & ! Input: [real(r8) (:,:)] interface level below a "z" level (m)
z => col%z , & ! Input: [real(r8) (:,:)] layer thickness (m)
ht_roof => lun%ht_roof , & ! Input: [real(r8) (:)] height of urban roof (m)
canyon_hwr => lun%canyon_hwr , & ! Input: [real(r8) (:)] ratio of building height to street hwidth (-)
wtlunit_roof => lun%wtlunit_roof , & ! Input: [real(r8) (:)] weight of roof with respect to landunit
urbpoi => lun%urbpoi , & ! Input: [logical (:)] true => landunit is an urban point
taf => temperature_inst%taf_lun , & ! Input: [real(r8) (:)] urban canopy air temperature (K)
tssbef => temperature_inst%t_ssbef_col , & ! Input: [real(r8) (:,:)] temperature at previous time step (K)
t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:)] soil temperature (K)
t_roof_inner => temperature_inst%t_roof_inner_lun , & ! InOut: [real(r8) (:)] roof inside surface temperature (K)
t_sunw_inner => temperature_inst%t_sunw_inner_lun , & ! InOut: [real(r8) (:)] sunwall inside surface temperature (K)
t_shdw_inner => temperature_inst%t_shdw_inner_lun , & ! InOut: [real(r8) (:)] shadewall inside surface temperature (K)
t_floor => temperature_inst%t_floor_lun , & ! InOut: [real(r8) (:)] floor temperature (K)
t_building => temperature_inst%t_building_lun , & ! InOut: [real(r8) (:)] internal building air temperature (K)
t_building_max => urbantv_inst%t_building_max , & ! Input: [real(r8) (:)] maximum internal building air temperature (K)
t_building_min => urbanparams_inst%t_building_min , & ! Input: [real(r8) (:)] minimum internal building air temperature (K)
eflx_building => energyflux_inst%eflx_building_lun , & ! Output: [real(r8) (:)] building heat flux from change in interior building air temperature (W/m**2)
eflx_urban_ac => energyflux_inst%eflx_urban_ac_lun , & ! Output: [real(r8) (:)] urban air conditioning flux (W/m**2)
eflx_urban_heat => energyflux_inst%eflx_urban_heat_lun & ! Output: [real(r8) (:)] urban heating flux (W/m**2)
)
! Get step size
dtime = get_step_size()
! 1. Save t_* at previous time step
! 2. Set convective heat transfer coefficients (Bueno et al. 2012, GMD).
! An alternative is Salamanca et al. 2010, TAC, where they are all set to 8 W m-2 K-1.
! See clm_varcon.F90
! 3. Set inner surface emissivities (Bueno et al. 2012, GMD).
! 4. Set concrete floor properties (Salamanca et al. 2010, TAC).
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
t_roof_inner_bef(l) = t_roof_inner(l)
t_sunw_inner_bef(l) = t_sunw_inner(l)
t_shdw_inner_bef(l) = t_shdw_inner(l)
t_floor_bef(l) = t_floor(l)
t_building_bef(l) = t_building(l)
if (t_roof_inner_bef(l) .le. t_building_bef(l)) then
hcv_roofi(l) = hcv_roof_enhanced
else
hcv_roofi(l) = hcv_roof
end if
if (t_floor_bef(l) .ge. t_building_bef(l)) then
hcv_floori(l) = hcv_floor_enhanced
else
hcv_floori(l) = hcv_floor
end if
hcv_sunwi(l) = hcv_sunw
hcv_shdwi(l) = hcv_shdw
em_roofi(l) = em_roof_int
em_sunwi(l) = em_sunw_int
em_shdwi(l) = em_shdw_int
em_floori(l) = em_floor_int
! Concrete floor thickness (m)
dz_floori(l) = dz_floor
! Concrete floor volumetric heat capacity (J m-3 K-1)
cp_floori(l) = cp_floor
! Intermediate calculation for concrete floor (W m-2 K-1)
cv_floori(l) = (dz_floori(l) * cp_floori(l)) / dtime
! Density of dry air at standard pressure and t_building (kg m-3)
rho_dair(l) = pstd / (rair*t_building_bef(l))
end if
end do
! Get terms from soil temperature equations to compute conduction flux
! Negative is toward surface - heat added
! Note that the conduction flux here is in W m-2 wall area but for purposes of solving the set of
! simultaneous equations this must be converted to W m-2 ground area. This is done below when
! setting up the equation coefficients.
do fc = 1,num_nolakec
c = filter_nolakec(fc)
l = clandunit(c)
if (urbpoi(l)) then
if (ctype(c) == icol_roof) then
zi_roof_innerl(l) = zi(c,nlevurb)
z_roof_innerl(l) = z(c,nlevurb)
t_roof_innerl_bef(l) = tssbef(c,nlevurb)
t_roof_innerl(l) = t_soisno(c,nlevurb)
tk_roof_innerl(l) = tk(c,nlevurb)
else if (ctype(c) == icol_sunwall) then
zi_sunw_innerl(l) = zi(c,nlevurb)
z_sunw_innerl(l) = z(c,nlevurb)
t_sunw_innerl_bef(l) = tssbef(c,nlevurb)
t_sunw_innerl(l) = t_soisno(c,nlevurb)
tk_sunw_innerl(l) = tk(c,nlevurb)
else if (ctype(c) == icol_shadewall) then
zi_shdw_innerl(l) = zi(c,nlevurb)
z_shdw_innerl(l) = z(c,nlevurb)
t_shdw_innerl_bef(l) = tssbef(c,nlevurb)
t_shdw_innerl(l) = t_soisno(c,nlevurb)
tk_shdw_innerl(l) = tk(c,nlevurb)
end if
end if
end do
! Calculate view factors
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
vf_rf(l) = sqrt(1._r8 + canyon_hwr(l)**2._r8) - canyon_hwr(l)
vf_fr(l) = vf_rf(l)
! This view factor implicitly converts from per unit wall area to per unit floor area
vf_wf(l) = 0.5_r8*(1._r8 - vf_rf(l))
! This view factor implicitly converts from per unit floor area to per unit wall area
vf_fw(l) = vf_wf(l) / canyon_hwr(l)
! This view factor implicitly converts from per unit roof area to per unit wall area
vf_rw(l) = vf_fw(l)
! This view factor implicitly converts from per unit wall area to per unit roof area
vf_wr(l) = vf_wf(l)
vf_ww(l) = 1._r8 - vf_rw(l) - vf_fw(l)
end if
end do
! error check -- make sure view factor sums to one for floor, wall, and roof
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
sum = vf_rf(l) + 2._r8*vf_wf(l)
if (abs(sum-1._r8) > 1.e-06_r8 ) then
write (iulog,*) 'urban floor view factor error',sum
write (iulog,*) 'clm model is stopping'
call endrun()
endif
sum = vf_rw(l) + vf_fw(l) + vf_ww(l)
if (abs(sum-1._r8) > 1.e-06_r8 ) then
write (iulog,*) 'urban wall view factor error',sum
write (iulog,*) 'clm model is stopping'
call endrun()
endif
sum = vf_fr(l) + vf_wr(l) + vf_wr(l)
if (abs(sum-1._r8) > 1.e-06_r8 ) then
write (iulog,*) 'urban roof view factor error',sum
write (iulog,*) 'clm model is stopping'
call endrun()
endif
endif
end do
n = neq
nrhs = 1
lda = neq
ldb = neq
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
! ROOF
a(1,1) = 0.5_r8*hcv_roofi(l) &
+ 0.5_r8*tk_roof_innerl(l)/(zi_roof_innerl(l) - z_roof_innerl(l)) &
+ 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8 &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fr(l)
a(1,2) = - 4._r8*em_roofi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wr(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l)
a(1,3) = - 4._r8*em_roofi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wr(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l)
a(1,4) = - 4._r8*em_roofi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fr(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wr(l)
a(1,5) = - 0.5_r8*hcv_roofi(l)
result(1) = 0.5_r8*tk_roof_innerl(l)*t_roof_innerl(l)/(zi_roof_innerl(l) - z_roof_innerl(l)) &
- 0.5_r8*tk_roof_innerl(l)*(t_roof_inner_bef(l)-t_roof_innerl_bef(l))/(zi_roof_innerl(l) &
- z_roof_innerl(l)) &
- 3._r8*em_roofi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wr(l) &
- 3._r8*em_roofi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wr(l) &
- 3._r8*em_roofi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fr(l) &
+ 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8 &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 0.5_r8*hcv_roofi(l)*(t_roof_inner_bef(l) - t_building_bef(l))
! SUNWALL
a(2,1) = - 4._r8*em_sunwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l)
a(2,2) = 0.5_r8*hcv_sunwi(l)*canyon_hwr(l) &
+ 0.5_r8*tk_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*canyon_hwr(l) &
+ 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8 &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)
a(2,3) = - 4._r8*em_sunwi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)
a(2,4) = - 4._r8*em_sunwi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_ww(l)
a(2,5) = - 0.5_r8*hcv_sunwi(l)*canyon_hwr(l)
result(2) = 0.5_r8*tk_sunw_innerl(l)*t_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*canyon_hwr(l) &
- 0.5_r8*tk_sunw_innerl(l)*(t_sunw_inner_bef(l)-t_sunw_innerl_bef(l))/(zi_sunw_innerl(l) &
- z_sunw_innerl(l))*canyon_hwr(l) &
- 3._r8*em_sunwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l) &
- 3._r8*em_sunwi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_ww(l) &
- 3._r8*em_sunwi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l) &
+ 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8 &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l)
! SHADEWALL
a(3,1) = - 4._r8*em_shdwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l)
a(3,2) = - 4._r8*em_shdwi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)
a(3,3) = 0.5_r8*hcv_shdwi(l)*canyon_hwr(l) &
+ 0.5_r8*tk_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*canyon_hwr(l) &
+ 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8 &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)
a(3,4) = - 4._r8*em_shdwi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_ww(l)
a(3,5) = - 0.5_r8*hcv_shdwi(l)*canyon_hwr(l)
result(3) = 0.5_r8*tk_shdw_innerl(l)*t_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*canyon_hwr(l) &
- 0.5_r8*tk_shdw_innerl(l)*(t_shdw_inner_bef(l)-t_shdw_innerl_bef(l))/(zi_shdw_innerl(l) &
- z_shdw_innerl(l))*canyon_hwr(l) &
- 3._r8*em_shdwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l) &
- 3._r8*em_shdwi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_ww(l) &
- 3._r8*em_shdwi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l) &
+ 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8 &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l)
! FLOOR
a(4,1) = - 4._r8*em_floori(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wf(l)
a(4,2) = - 4._r8*em_floori(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l)
a(4,3) = - 4._r8*em_floori(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wf(l)
a(4,4) = (cv_floori(l) + 0.5_r8*hcv_floori(l)) &
+ 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8 &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wf(l)
a(4,5) = - 0.5_r8*hcv_floori(l)
result(4) = cv_floori(l)*t_floor_bef(l) &
- 3._r8*em_floori(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rf(l) &
- 3._r8*em_floori(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wf(l) &
- 3._r8*em_floori(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wf(l) &
+ 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8 &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 3._r8*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 3._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 3._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 3._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 0.5_r8*hcv_floori(l)*(t_floor_bef(l) - t_building_bef(l))
! Building air temperature
a(5,1) = - 0.5_r8*hcv_roofi(l)
a(5,2) = - 0.5_r8*hcv_sunwi(l)*canyon_hwr(l)
a(5,3) = - 0.5_r8*hcv_shdwi(l)*canyon_hwr(l)
a(5,4) = - 0.5_r8*hcv_floori(l)
a(5,5) = ((ht_roof(l)*rho_dair(l)*cpair)/dtime) + &
((ht_roof(l)*vent_ach)/3600._r8)*rho_dair(l)*cpair + &
0.5_r8*hcv_roofi(l) + &
0.5_r8*hcv_sunwi(l)*canyon_hwr(l) + &
0.5_r8*hcv_shdwi(l)*canyon_hwr(l) + &
0.5_r8*hcv_floori(l)
result(5) = (ht_roof(l)*rho_dair(l)*cpair/dtime)*t_building_bef(l) &
+ ((ht_roof(l)*vent_ach)/3600._r8)*rho_dair(l)*cpair*taf(l) &
+ 0.5_r8*hcv_roofi(l)*(t_roof_inner_bef(l) - t_building_bef(l)) &
+ 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
+ 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
+ 0.5_r8*hcv_floori(l)*(t_floor_bef(l) - t_building_bef(l))
! Solve equations
call dgesv(n, nrhs, a, lda, ipiv, result, ldb, info)
! If dgesv fails, abort
if (info /= 0) then
write(iulog,*)'fl: ',fl
write(iulog,*)'l: ',l
write(iulog,*)'dgesv info: ',info
write (iulog,*) 'dgesv error'
write (iulog,*) 'clm model is stopping'
call endrun()
end if
! Assign new temperatures
t_roof_inner(l) = result(1)
t_sunw_inner(l) = result(2)
t_shdw_inner(l) = result(3)
t_floor(l) = result(4)
t_building(l) = result(5)
end if
end do
! Energy balance checks
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
qrd_roof(l) = - em_roofi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wr(l) &
- 4._r8*em_roofi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wr(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- em_roofi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wr(l) &
- 4._r8*em_roofi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wr(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- em_roofi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fr(l) &
- 4._r8*em_roofi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fr(l)*(t_floor(l) - t_floor_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wr(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wr(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l)*(1._r8-em_floori(l))*vf_fr(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wr(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wr(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l)*(1._r8-em_floori(l))*vf_fr(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wr(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wr(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wr(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wr(l)*(t_floor(l) &
- t_floor_bef(l)) &
+ em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8 &
+ 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*(t_roof_inner(l) - t_roof_inner_bef(l))
qrd_sunw(l) = - em_sunwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l) &
- 4._r8*em_sunwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- em_sunwi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_ww(l) &
- 4._r8*em_sunwi(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_ww(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- em_sunwi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l) &
- 4._r8*em_sunwi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(t_floor(l) - t_floor_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_ww(l)*(1._r8-em_shdwi(l))*vf_ww(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rw(l)*(1._r8-em_shdwi(l))*vf_ww(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_shdwi(l))*vf_ww(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fw(l)*(1._r8-em_shdwi(l))*vf_ww(l)*(t_floor(l) &
- t_floor_bef(l)) &
+ em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8 &
+ 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*(t_sunw_inner(l) - t_sunw_inner_bef(l))
qrd_shdw(l) = - em_shdwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rw(l) &
- 4._r8*em_shdwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- em_shdwi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_ww(l) &
- 4._r8*em_shdwi(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_ww(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- em_shdwi(l)*em_floori(l)*sb*t_floor_bef(l)**4._r8*vf_fw(l) &
- 4._r8*em_shdwi(l)*em_floori(l)*sb*t_floor_bef(l)**3._r8*vf_fw(l)*(t_floor(l) - t_floor_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_ww(l)*(1._r8-em_sunwi(l))*vf_ww(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_wf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rw(l)*(1._r8-em_sunwi(l))*vf_ww(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rf(l)*(1._r8-em_floori(l))*vf_fw(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fr(l)*(1._r8-em_roofi(l))*vf_rw(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_sunwi(l))*vf_ww(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fw(l)*(1._r8-em_sunwi(l))*vf_ww(l)*(t_floor(l) &
- t_floor_bef(l)) &
+ em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8 &
+ 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*(t_shdw_inner(l) - t_shdw_inner_bef(l))
qrd_floor(l) = - em_floori(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8*vf_rf(l) &
- 4._r8*em_floori(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- em_floori(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8*vf_wf(l) &
- 4._r8*em_floori(l)*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3._r8*vf_wf(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- em_floori(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8*vf_wf(l) &
- 4._r8*em_floori(l)*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3._r8*vf_wf(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fr(l)*(1._r8-em_roofi(l))*vf_rf(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fw(l)*(1._r8-em_sunwi(l))*vf_wf(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_floori(l)*sb*t_floor_bef(l)**4._r8)*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*vf_fw(l)*(1._r8-em_shdwi(l))*vf_wf(l)*(t_floor(l) &
- t_floor_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_ww(l)*(1._r8-em_shdwi(l))*vf_wf(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_sunwi(l)*sb*t_sunw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 4._r8*em_sunwi(l)*sb*t_sunw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l)*(t_sunw_inner(l) &
- t_sunw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_wr(l)*(1._r8-em_roofi(l))*vf_rf(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_shdwi(l)*sb*t_shdw_inner_bef(l)**4._r8)*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 4._r8*em_shdwi(l)*sb*t_shdw_inner_bef(l)**3.*vf_ww(l)*(1._r8-em_sunwi(l))*vf_wf(l)*(t_shdw_inner(l) &
- t_shdw_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wf(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rw(l)*(1._r8-em_sunwi(l))*vf_wf(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
- (em_roofi(l)*sb*t_roof_inner_bef(l)**4._r8)*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wf(l) &
- 4._r8*em_roofi(l)*sb*t_roof_inner_bef(l)**3.*vf_rw(l)*(1._r8-em_shdwi(l))*vf_wf(l)*(t_roof_inner(l) &
- t_roof_inner_bef(l)) &
+ em_floori(l)*sb*t_floor_bef(l)**4._r8 &
+ 4._r8*em_floori(l)*sb*t_floor_bef(l)**3.*(t_floor(l) - t_floor_bef(l))
qrd_building(l) = qrd_roof(l) + canyon_hwr(l)*(qrd_sunw(l) + qrd_shdw(l)) + qrd_floor(l)
if (abs(qrd_building(l)) > .10_r8 ) then
write (iulog,*) 'urban inside building net longwave radiation balance error ',qrd_building(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
qcv_roof(l) = 0.5_r8*hcv_roofi(l)*(t_roof_inner(l) - t_building(l)) + 0.5_r8*hcv_roofi(l)*(t_roof_inner_bef(l) &
- t_building_bef(l))
qcd_roof(l) = 0.5_r8*tk_roof_innerl(l)*(t_roof_inner(l) - t_roof_innerl(l))/(zi_roof_innerl(l) - z_roof_innerl(l)) &
+ 0.5_r8*tk_roof_innerl(l)*(t_roof_inner_bef(l) - t_roof_innerl_bef(l))/(zi_roof_innerl(l) &
- z_roof_innerl(l))
enrgy_bal_roof(l) = qrd_roof(l) + qcv_roof(l) + qcd_roof(l)
if (abs(enrgy_bal_roof(l)) > .10_r8 ) then
write (iulog,*) 'urban inside roof energy balance error ',enrgy_bal_roof(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
qcv_sunw(l) = 0.5_r8*hcv_sunwi(l)*(t_sunw_inner(l) - t_building(l)) + 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) &
- t_building_bef(l))
qcd_sunw(l) = 0.5_r8*tk_sunw_innerl(l)*(t_sunw_inner(l) - t_sunw_innerl(l))/(zi_sunw_innerl(l) - z_sunw_innerl(l)) &
+ 0.5_r8*tk_sunw_innerl(l)*(t_sunw_inner_bef(l) - t_sunw_innerl_bef(l))/(zi_sunw_innerl(l) &
- z_sunw_innerl(l))
enrgy_bal_sunw(l) = qrd_sunw(l) + qcv_sunw(l)*canyon_hwr(l) + qcd_sunw(l)*canyon_hwr(l)
if (abs(enrgy_bal_sunw(l)) > .10_r8 ) then
write (iulog,*) 'urban inside sunwall energy balance error ',enrgy_bal_sunw(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
qcv_shdw(l) = 0.5_r8*hcv_shdwi(l)*(t_shdw_inner(l) - t_building(l)) + 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) &
- t_building_bef(l))
qcd_shdw(l) = 0.5_r8*tk_shdw_innerl(l)*(t_shdw_inner(l) - t_shdw_innerl(l))/(zi_shdw_innerl(l) - z_shdw_innerl(l)) &
+ 0.5_r8*tk_shdw_innerl(l)*(t_shdw_inner_bef(l) - t_shdw_innerl_bef(l))/(zi_shdw_innerl(l) &
- z_shdw_innerl(l))
enrgy_bal_shdw(l) = qrd_shdw(l) + qcv_shdw(l)*canyon_hwr(l) + qcd_shdw(l)*canyon_hwr(l)
if (abs(enrgy_bal_shdw(l)) > .10_r8 ) then
write (iulog,*) 'urban inside shadewall energy balance error ',enrgy_bal_shdw(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
qcv_floor(l) = 0.5_r8*hcv_floori(l)*(t_floor(l) - t_building(l)) + 0.5_r8*hcv_floori(l)*(t_floor_bef(l) &
- t_building_bef(l))
qcd_floor(l) = cv_floori(l)*(t_floor(l) - t_floor_bef(l))
enrgy_bal_floor(l) = qrd_floor(l) + qcv_floor(l) + qcd_floor(l)
if (abs(enrgy_bal_floor(l)) > .10_r8 ) then
write (iulog,*) 'urban inside floor energy balance error ',enrgy_bal_floor(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
enrgy_bal_buildair(l) = (ht_roof(l)*rho_dair(l)*cpair/dtime)*(t_building(l) - t_building_bef(l)) &
- ht_roof(l)*(vent_ach/3600._r8)*rho_dair(l)*cpair*(taf(l) - t_building(l)) &
- 0.5_r8*hcv_roofi(l)*(t_roof_inner(l) - t_building(l)) &
- 0.5_r8*hcv_roofi(l)*(t_roof_inner_bef(l) - t_building_bef(l)) &
- 0.5_r8*hcv_sunwi(l)*(t_sunw_inner(l) - t_building(l))*canyon_hwr(l) &
- 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner(l) - t_building(l))*canyon_hwr(l) &
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
- 0.5_r8*hcv_floori(l)*(t_floor(l) - t_building(l)) &
- 0.5_r8*hcv_floori(l)*(t_floor_bef(l) - t_building_bef(l))
if (abs(enrgy_bal_buildair(l)) > .10_r8 ) then
write (iulog,*) 'urban building air energy balance error ',enrgy_bal_buildair(l)
write (iulog,*) 'clm model is stopping'
call endrun()
end if
end if
end do
! Restrict internal building air temperature to between min and max
! Calculate heating or air conditioning flux from energy required to change
! internal building air temperature to t_building_min or t_building_max.
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then
t_building_bef_hac(l) = t_building(l)
! rho_dair(l) = pstd / (rair*t_building(l))
if (t_building_bef_hac(l) > t_building_max(l)) then
t_building(l) = t_building_max(l)
eflx_urban_ac(l) = wtlunit_roof(l) * abs( (ht_roof(l) * rho_dair(l) * cpair / dtime) * t_building(l) &
- (ht_roof(l) * rho_dair(l) * cpair / dtime) * t_building_bef_hac(l) )
else if (t_building_bef_hac(l) < t_building_min(l)) then
t_building(l) = t_building_min(l)
eflx_urban_heat(l) = wtlunit_roof(l) * abs( (ht_roof(l) * rho_dair(l) * cpair / dtime) * t_building(l) &
- (ht_roof(l) * rho_dair(l) * cpair / dtime) * t_building_bef_hac(l) )
else
eflx_urban_ac(l) = 0._r8
eflx_urban_heat(l) = 0._r8
end if
else
eflx_urban_ac(l) = 0._r8
eflx_urban_heat(l) = 0._r8
end if
eflx_building(l) = wtlunit_roof(l) * (ht_roof(l) * rho_dair(l)*cpair/dtime) * (t_building(l) - t_building_bef(l))
end if
end do
end associate
end subroutine BuildingTemperature
!-----------------------------------------------------------------------
end module UrbBuildTempOleson2015Mod