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

Assumption that building width (Wb) equals street width can be relaxed and Wb can be derived from morphology #803

Closed
olyson opened this issue Sep 12, 2019 · 13 comments
Assignees
Labels
bug something is working incorrectly science Enhancement to or bug impacting science

Comments

@olyson
Copy link
Contributor

olyson commented Sep 12, 2019

Brief summary of bug

There is an implicit assumption in the urban building energy model that building width equals street width. However, this assumption can/should be relaxed and building width can be derived from the morphology dataset.

General bug information

CTSM version you are using: release-clm

Does this bug cause significantly incorrect results in the model's science? Yes for urban (greater than roundoff), No for general climate simulations (same gridcell/regional/global climate).

Configurations affected: All configurations with urban landunits

Details of bug

There is an implicit assumption in the urban building energy model that building width equals street width. However, this assumption can be relaxed and building width can be derived from the Jackson morphology dataset. Specifically, the equations which use H/Ws (building height over street width) should use H/Wb (building height over building width). Building width can be derived from (Ws*froof)/(1-froof), where froof is the roof (building) fraction.
A simulation with the new equations indicate that as expected, for froof=0.5, answers are identical (because H/Ws = H/Wb). Differences in other cases depend on the degree of departure from froof=0.5 and the presence and magnitude of space heating and air conditioning. A full assessment of differences by variable, season, density type, and region using ctsm1.0.dev098 can be found in the attached file. This uses the default CLM5 urban datasets.

CTSM98_BUILDENERGY_Analysis_Pub.xlsx

A similar assessment using release-cesm2.0.1 but with the optional new urban datasets (Oleson and Feddema 2019) can be found in this attached file:

CLM50_BUILDENERGY_Analysis_Pub.xlsx

I have SourceMods for this currently, but can put together a branch. It would be nice to get this into the release branch, in particular for the CESM2 large ensemble runs, since it's not going to affect climate and would make the urban output more useful. Along with the new urban properties data (see PR#591).

@olyson olyson self-assigned this Sep 12, 2019
@ekluzek ekluzek self-assigned this Sep 12, 2019
@ekluzek
Copy link
Collaborator

ekluzek commented Sep 12, 2019

@olyson has this been approved to go on the release branch and for cesm2.1.2? I was thinking urban properties update would only go on master. Is that approved to go on the release branch as well?

Also if you could point me to your case where you have the SourceMods that would be helpful. Thanks.

@olyson
Copy link
Contributor Author

olyson commented Sep 12, 2019

No, hasn't been approved, I'm simply proposing that these go on the release branch. I'm mainly thinking about the large ensemble that is coming up. I'm not sure whether a decision has been made about whether "bug" fixes will be allowed for the large ensemble tag. I will ask Dave when I get a chance.
The SourceMods are confined to one module:
/glade/work/oleson/release-cesm2.0.1/cime/scripts/clm50_cesm201R_1deg_GSWP3V1_CON_FIXBUILDENERGY_2000/SourceMods/src.clm/UrbBuildTempOleson2015Mod.F90

@olyson
Copy link
Contributor Author

olyson commented Sep 12, 2019

The word from Dave is that these should only go on master.

@billsacks billsacks added the bug something is working incorrectly label Sep 12, 2019
@ekluzek
Copy link
Collaborator

ekluzek commented Sep 12, 2019

OK, the changes are to increase nlevurb to 10. Which is done because of the updated surface dataset that is used in the case...

fsurdat = '/gpfs/fs1/p/cgd/tss/people/oleson/urban_sfcdata/Feddema_urban_data_080410/HIGH_RES_VERSION2/URBAN_PROPERTIES_THESIS_TOOL/JUN_22_2018/surfdata_0.9x1.25_16pfts_Irrig_CMIP6_simyr2000_Control_c190405.nc'

It also uses an update urbantv dataset.

And then this change in the code...

SourceMods/src.clm> diff -wbc /gpfs/fs1/work/oleson/release-cesm2.0.1/components/clm/src/biogeophys/UrbBuildTempOleson2015Mod.F90 .
*** /gpfs/fs1/work/oleson/release-cesm2.0.1/components/clm/src/biogeophys/UrbBuildTempOleson2015Mod.F90	2018-08-23 11:52:37.943164484 -0600
--- ./UrbBuildTempOleson2015Mod.F90	2019-09-12 14:12:02.211103587 -0600
***************
*** 230,235 ****
--- 230,238 ----
      integer, parameter :: neq = 5          ! number of equation/unknowns
      integer  :: fc,fl,c,l                  ! indices
      real(r8) :: dtime                      ! land model time step (s)
+ !KO
+     real(r8) :: building_hwr(bounds%begl:bounds%endl)      ! building height to building width ratio (-)
+ !KO
      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)              
***************
*** 342,347 ****
--- 345,351 ----
      !    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
***************
*** 374,386 ****
           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
--- 378,397 ----
           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))
+ !KO
+          ! Building height to building width ratio
+          building_hwr(l) = canyon_hwr(l)*(1._r8-wtlunit_roof(l))/wtlunit_roof(l)
+ !KO
         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
! !KO    ! simultaneous equations this must be converted to W m-2 ground area. This is done below when 
! !KO
!     ! simultaneous equations this must be converted to W m-2 floor area. This is done below when 
! !KO
      ! setting up the equation coefficients.
  
      do fc = 1,num_nolakec
***************
*** 414,427 ****
         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)
--- 425,444 ----
         l = filter_urbanl(fl)
         if (urbpoi(l)) then
  
! !KO         vf_rf(l) = sqrt(1._r8 + canyon_hwr(l)**2._r8) - canyon_hwr(l)
! !KO
!          vf_rf(l) = sqrt(1._r8 + building_hwr(l)**2._r8) - building_hwr(l)
! !KO
           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
! !KO         vf_fw(l) = vf_wf(l) / canyon_hwr(l)
! !KO
!          vf_fw(l) = vf_wf(l) / building_hwr(l)
! !KO
  
           ! This view factor implicitly converts from per unit roof area to per unit wall area
           vf_rw(l)  = vf_fw(l)
***************
*** 516,523 ****
                    - 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) &
--- 533,544 ----
                    - 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)
  
! !KO         a(2,2) =   0.5_r8*hcv_sunwi(l)*canyon_hwr(l) &
! !KO                  + 0.5_r8*tk_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*canyon_hwr(l) &
! !KO
!          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) &
! !KO
                    + 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) &
***************
*** 530,540 ****
           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) &
--- 551,569 ----
           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)
! !KO         a(2,5) = - 0.5_r8*hcv_sunwi(l)*canyon_hwr(l)
! !KO
!          a(2,5) = - 0.5_r8*hcv_sunwi(l)*building_hwr(l)
! !KO
! 
! !KO         result(2) =   0.5_r8*tk_sunw_innerl(l)*t_sunw_innerl(l)/(zi_sunw_innerl(l) - z_sunw_innerl(l))*canyon_hwr(l) &
! !KO                     - 0.5_r8*tk_sunw_innerl(l)*(t_sunw_inner_bef(l)-t_sunw_innerl_bef(l))/(zi_sunw_innerl(l) &
! !KO                     - z_sunw_innerl(l))*canyon_hwr(l) &
! !KO
!          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))*building_hwr(l) &
! !KO
                       - 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) &
***************
*** 548,554 ****
                       - 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) &
--- 577,586 ----
                       - 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) &
! !KO                     - 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l)
! !KO
!                      - 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*building_hwr(l)
! !KO
  
           ! SHADEWALL
           a(3,1) = - 4._r8*em_shdwi(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rw(l) &
***************
*** 559,566 ****
                    - 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) &
--- 591,602 ----
                    - 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)
  
! !KO         a(3,3) =   0.5_r8*hcv_shdwi(l)*canyon_hwr(l) &
! !KO                  + 0.5_r8*tk_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*canyon_hwr(l) &
! !KO
!          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) &
! !KO
                    + 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) &
***************
*** 570,580 ****
                    - 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) &
--- 606,624 ----
                    - 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)
  
! !KO         a(3,5) = - 0.5_r8*hcv_shdwi(l)*canyon_hwr(l)
! !KO
!          a(3,5) = - 0.5_r8*hcv_shdwi(l)*building_hwr(l)
! !KO
! 
! !KO         result(3) =   0.5_r8*tk_shdw_innerl(l)*t_shdw_innerl(l)/(zi_shdw_innerl(l) - z_shdw_innerl(l))*canyon_hwr(l) &
! !KO                     - 0.5_r8*tk_shdw_innerl(l)*(t_shdw_inner_bef(l)-t_shdw_innerl_bef(l))/(zi_shdw_innerl(l) &
! !KO                     - z_shdw_innerl(l))*canyon_hwr(l) &
! !KO
!          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))*building_hwr(l) &
! !KO
                       - 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) &
***************
*** 588,594 ****
                       - 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) &
--- 632,641 ----
                       - 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) &
! !KO                     - 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l)
! !KO
!                      - 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*building_hwr(l)
! !KO
  
           ! FLOOR
           a(4,1) = - 4._r8*em_floori(l)*em_roofi(l)*sb*t_roof_inner_bef(l)**3._r8*vf_rf(l) &
***************
*** 629,652 ****
  
           ! 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
--- 676,713 ----
  
           ! Building air temperature
           a(5,1) = - 0.5_r8*hcv_roofi(l)
! !KO         a(5,2) = - 0.5_r8*hcv_sunwi(l)*canyon_hwr(l)
! !KO
!          a(5,2) = - 0.5_r8*hcv_sunwi(l)*building_hwr(l)
! !KO
! 
! !KO         a(5,3) = - 0.5_r8*hcv_shdwi(l)*canyon_hwr(l)
! !KO
!          a(5,3) = - 0.5_r8*hcv_shdwi(l)*building_hwr(l)
! !KO
  
           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) + &
! !KO                   0.5_r8*hcv_sunwi(l)*canyon_hwr(l) + &
! !KO                   0.5_r8*hcv_shdwi(l)*canyon_hwr(l) + &
! !KO
!                    0.5_r8*hcv_sunwi(l)*building_hwr(l) + &
!                    0.5_r8*hcv_shdwi(l)*building_hwr(l) + &
! !KO
                     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)) &
! !KO                      + 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
! !KO                      + 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
! !KO
!                       + 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) &
! !KO
                        + 0.5_r8*hcv_floori(l)*(t_floor_bef(l) - t_building_bef(l))
  
           ! Solve equations
***************
*** 827,833 ****
                          + 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)
--- 888,894 ----
                          + 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) + 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)
***************
*** 852,858 ****
           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'
--- 913,922 ----
           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))
! !KO         enrgy_bal_sunw(l) = qrd_sunw(l) + qcv_sunw(l)*canyon_hwr(l) + qcd_sunw(l)*canyon_hwr(l)
! !KO
!          enrgy_bal_sunw(l) = qrd_sunw(l) + qcv_sunw(l)*building_hwr(l) + qcd_sunw(l)*building_hwr(l)
! !KO
           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'
***************
*** 864,870 ****
           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'
--- 928,937 ----
           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))
! !KO         enrgy_bal_shdw(l) = qrd_shdw(l) + qcv_shdw(l)*canyon_hwr(l) + qcd_shdw(l)*canyon_hwr(l)
! !KO
!          enrgy_bal_shdw(l) = qrd_shdw(l) + qcv_shdw(l)*building_hwr(l) + qcd_shdw(l)*building_hwr(l)
! !KO
           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'
***************
*** 885,894 ****
                                   - 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
--- 952,967 ----
                                   - 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)) &
! !KO                                 - 0.5_r8*hcv_sunwi(l)*(t_sunw_inner(l) - t_building(l))*canyon_hwr(l) &
! !KO                                 - 0.5_r8*hcv_sunwi(l)*(t_sunw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
! !KO                                 - 0.5_r8*hcv_shdwi(l)*(t_shdw_inner(l) - t_building(l))*canyon_hwr(l) &
! !KO                                 - 0.5_r8*hcv_shdwi(l)*(t_shdw_inner_bef(l) - t_building_bef(l))*canyon_hwr(l) &
! !KO
!                                  - 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) &
! !KO
                                   - 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

@olyson
Copy link
Contributor Author

olyson commented Sep 12, 2019

Note that the change to nlevurb=10, the updated surface dataset, and the updated urbantv file are part of PR #591.

@ekluzek ekluzek added the next this should get some attention in the next week or two. Normally each Thursday SE meeting. label Sep 23, 2019
@billsacks billsacks removed the next this should get some attention in the next week or two. Normally each Thursday SE meeting. label Sep 23, 2019
olyson added a commit to olyson/ctsm that referenced this issue Jun 8, 2020
Performance verified through spreadsheet CTSM98_BUILDENERGY_Analysis_Pub
@olyson
Copy link
Contributor Author

olyson commented Jun 8, 2020

I created a branch for the code changes:

https://github.com/olyson/ctsm/tree/fixbuildenergy

@ekluzek
Copy link
Collaborator

ekluzek commented Jun 8, 2020

@olyson I looked at your branch, the changes are pretty small and straightforward. And my understanding is that they won't change answers for the current datasets? They just allow different behavior for the future? Is that right?

Do you want to turn the branch into a PR, and we can start moving it to master?

@olyson
Copy link
Contributor Author

olyson commented Jun 8, 2020

No, they will change answers for the current datasets. The code changes could come to master either by themselves or as part of PR #591 (the change to nlevurb=10, the updated surface dataset, and the updated urbantv file). It would be nice to have the code changes on master sooner rather than later. So I will turn the branch into a PR if that is ok.

@ekluzek
Copy link
Collaborator

ekluzek commented Jun 8, 2020

@olyson it looks like PR #591 changes answers, so it would be good to have this separated from that with the bit-for-bit part of it as a separate tag. If there's other bit-for-bit parts of #591, we could lift them out into bit-for-bit PR's as well.

@olyson
Copy link
Contributor Author

olyson commented Jun 8, 2020

There actually aren't any bit-for-bit parts of #591. And this branch with the code changes isn't bit-for-bit either, it's a bug fix.

@ekluzek
Copy link
Collaborator

ekluzek commented Jun 8, 2020

Ahh, OK. Thanks for that clarification. As I read #591 I see it as mostly a new feature for CTSM5_1. The new datasets would be triggered with CTSM5_1, but the old used with CLM5_0. It's good to keep things that are a bug fix that change other versions (like clm4_5, and clm5_0) outside of the new feature for a new version. Does that sound right to you? I may need to discuss this with your more to figure this out...

If this changes answers as a bug fix, do we know that the changes are reasonably small?

@olyson
Copy link
Contributor Author

olyson commented Jun 8, 2020

Let's talk at your convenience (I know you are busy with the release and other things this week).

@olyson
Copy link
Contributor Author

olyson commented Jun 8, 2020

The answer changes due to the bug fix are very small, nothing statistically significant. So it seems like the bug fix itself could go onto master.

http://webext.cgd.ucar.edu/I2000/clm50_ctsm10d098_1deg_GSWP3V1_CON_FIXBUILDENERGY_2000/lnd/clm50_ctsm10d098_1deg_GSWP3V1_CON_FIXBUILDENERGY_2000.1991_2010-clm50_ctsm10d098_1deg_GSWP3V1_CON_2000.1991_2010/setsIndex.html

ekluzek added a commit that referenced this issue Jun 25, 2020
@ekluzek ekluzek closed this as completed Jun 26, 2020
@samsrabin samsrabin added the science Enhancement to or bug impacting science label Aug 8, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug something is working incorrectly science Enhancement to or bug impacting science
Projects
None yet
Development

No branches or pull requests

4 participants