Merge pull request #1034 from olyson/fixbuildenergy

Building energy fix (issue #803).

Author: ekluzek

src/biogeophys/UrbBuildTempOleson2015Mod.F90

@@ -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)              
@@ -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
@@ -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
@@ -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)
@@ -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) &
@@ -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) &
@@ -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) &
@@ -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) &
@@ -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) &
@@ -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) &
@@ -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
@@ -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)
@@ -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'
@@ -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'
@@ -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