Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Building energy fix (issue #803). #1034

Merged
merged 2 commits into from
Jun 25, 2020
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 32 additions & 28 deletions src/biogeophys/UrbBuildTempOleson2015Mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -229,6 +229,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
integer, parameter :: neq = 5 ! number of equation/unknowns
integer :: fc,fl,c,l ! indices
real(r8) :: dtime ! land model time step (s)
real(r8) :: building_hwr(bounds%begl:bounds%endl) ! building height to building width ratio (-)
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)
Expand Down Expand Up @@ -341,6 +342,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
! See clm_varcon.F90
! 3. Set inner surface emissivities (Bueno et al. 2012, GMD).
! 4. Set concrete floor properties (Salamanca et al. 2010, TAC).
! 5. Calculate building height to building width ratio
do fl = 1,num_urbanl
l = filter_urbanl(fl)
if (urbpoi(l)) then
Expand Down Expand Up @@ -373,13 +375,15 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
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))
! Building height to building width ratio
building_hwr(l) = canyon_hwr(l)*(1._r8-wtlunit_roof(l))/wtlunit_roof(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
! simultaneous equations this must be converted to W m-2 floor area. This is done below when
! setting up the equation coefficients.

do fc = 1,num_nolakec
Expand Down Expand Up @@ -413,14 +417,14 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
l = filter_urbanl(fl)
if (urbpoi(l)) then

vf_rf(l) = sqrt(1._r8 + canyon_hwr(l)**2._r8) - canyon_hwr(l)
vf_rf(l) = sqrt(1._r8 + building_hwr(l)**2._r8) - building_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)
vf_fw(l) = vf_wf(l) / building_hwr(l)

! This view factor implicitly converts from per unit roof area to per unit wall area
vf_rw(l) = vf_fw(l)
Expand Down Expand Up @@ -515,8 +519,8 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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) &
a(2,2) = 0.5_r8*hcv_sunwi(l)*building_hwr(l) &
+ 0.5_r8*tk_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*building_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) &
Expand All @@ -529,11 +533,11 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
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)
a(2,5) = - 0.5_r8*hcv_sunwi(l)*building_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) &
result(2) = 0.5_r8*tk_sunw_innerl(l)*t_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*building_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) &
- z_sunw_innerl(l))*building_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) &
Expand All @@ -547,7 +551,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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)
- 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*building_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) &
Expand All @@ -558,8 +562,8 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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) &
a(3,3) = 0.5_r8*hcv_shdwi(l)*building_hwr(l) &
+ 0.5_r8*tk_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*building_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) &
Expand All @@ -569,11 +573,11 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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)
a(3,5) = - 0.5_r8*hcv_shdwi(l)*building_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) &
result(3) = 0.5_r8*tk_shdw_innerl(l)*t_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*building_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) &
- z_shdw_innerl(l))*building_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) &
Expand All @@ -587,7 +591,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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)
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*building_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) &
Expand Down Expand Up @@ -628,24 +632,24 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,

! 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,2) = - 0.5_r8*hcv_sunwi(l)*building_hwr(l)

a(5,3) = - 0.5_r8*hcv_shdwi(l)*canyon_hwr(l)
a(5,3) = - 0.5_r8*hcv_shdwi(l)*building_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_sunwi(l)*building_hwr(l) + &
0.5_r8*hcv_shdwi(l)*building_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_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*building_hwr(l) &
+ 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*building_hwr(l) &
+ 0.5_r8*hcv_floori(l)*(t_floor_bef(l) - t_building_bef(l))

! Solve equations
Expand Down Expand Up @@ -826,7 +830,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
+ 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)
qrd_building(l) = qrd_roof(l) + building_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)
Expand All @@ -851,7 +855,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
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)
enrgy_bal_sunw(l) = qrd_sunw(l) + qcv_sunw(l)*building_hwr(l) + qcd_sunw(l)*building_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'
Expand All @@ -863,7 +867,7 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
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)
enrgy_bal_shdw(l) = qrd_shdw(l) + qcv_shdw(l)*building_hwr(l) + qcd_shdw(l)*building_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'
Expand All @@ -884,10 +888,10 @@ subroutine BuildingTemperature (bounds, num_urbanl, filter_urbanl, num_nolakec,
- 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_sunwi(l)*(t_sunw_inner(l) - t_building(l))*building_hwr(l) &
- 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*building_hwr(l) &
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner(l) - t_building(l))*building_hwr(l) &
- 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*building_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
Expand Down