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.

Sign up for GitHub

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

Small answer changes in preparation for adding option for bio-mass heat storage #1241

Merged
merged 11 commits into from
Dec 30, 2020
Merged

Small answer changes in preparation for adding option for bio-mass heat storage #1241

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
b9b9c21
Preparation for biomass heat storage by rearrnaging terms a bit for d…
ekluzek Dec 29, 2020
94f8223
Merge commit 'b9b9c216' into bhsprep_roffanschange
ekluzek Dec 29, 2020
49b0e9f
Add checks that differences are within roundoff, get it working corre…
ekluzek Dec 29, 2020
eb7bdfe
Add new clm51 compsets with Sp mode, change one test to use it and ad…
ekluzek Dec 29, 2020
501f5f8
Add some more terms to match the full BHS work, increase tolerance to…
ekluzek Dec 29, 2020
e3797fc
Do relative diffs rather than absolute, lower tolerance to 1.e-14, ra…
ekluzek Dec 30, 2020
22fcf15
Remove temporary tests for terms being only different by relative dif…
ekluzek Dec 30, 2020
5dfd628
Correct number of tests
ekluzek Dec 30, 2020
5bb7e59
Work on ChangeLog
ekluzek Dec 30, 2020
f0abdd9
Update changelog
ekluzek Dec 30, 2020
7410b45
Update change files
ekluzek Dec 30, 2020
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
35 changes: 35 additions & 0 deletions cime_config/config_compsets.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,11 @@

<!-- I single point forcing -->

<compset>
<alias>I1PtClm51SpRs</alias>
<lname>2000_DATM%1PT_CLM51%SP_SICE_SOCN_SROF_SGLC_SWAV</lname>
</compset>

<compset>
<alias>I1PtClm50SpRs</alias>
<lname>2000_DATM%1PT_CLM50%SP_SICE_SOCN_SROF_SGLC_SWAV</lname>
Expand All @@ -59,6 +64,17 @@
<lname>2000_DATM%GSWP3v1_CLM50%SP_SICE_SOCN_SROF_SGLC_SWAV</lname>
</compset>

<compset>
<alias>I2000Clm51Sp</alias>
<lname>2000_DATM%GSWP3v1_CLM51%SP_SICE_SOCN_MOSART_SGLC_SWAV</lname>
</compset>

<!-- Primarily for testing -->
<compset>
<alias>I2000Clm51SpRs</alias>
<lname>2000_DATM%GSWP3v1_CLM51%SP_SICE_SOCN_SROF_SGLC_SWAV</lname>
</compset>

<!-- Primarily for testing -->
<compset>
<alias>I2000Clm50SpRtm</alias>
Expand Down Expand Up @@ -128,6 +144,11 @@
<science_support grid="f19_g17"/>
</compset>

<compset>
<alias>I1850Clm51Sp</alias>
<lname>1850_DATM%GSWP3v1_CLM51%SP_SICE_SOCN_MOSART_SGLC_SWAV</lname>
</compset>

<compset>
<alias>I1850Clm51Bgc</alias>
<lname>1850_DATM%GSWP3v1_CLM51%BGC_SICE_SOCN_MOSART_SGLC_SWAV</lname>
Expand Down Expand Up @@ -177,6 +198,10 @@
</compset>

<!---I FATES compsets -->
<compset>
<alias>I2000Clm51Fates</alias>
<lname>2000_DATM%GSWP3v1_CLM51%FATES_SICE_SOCN_MOSART_SGLC_SWAV</lname>
</compset>
<compset>
<alias>I2000Clm50Fates</alias>
<lname>2000_DATM%GSWP3v1_CLM50%FATES_SICE_SOCN_MOSART_SGLC_SWAV</lname>
Expand Down Expand Up @@ -213,6 +238,16 @@
<science_support grid="f19_g17"/>
</compset>

<compset>
<alias>I1850Clm51SpNoAnthro</alias>
<lname>1850_DATM%GSWP3v1_CLM51%SP-NOANTHRO_SICE_SOCN_MOSART_SGLC_SWAV</lname>
</compset>

<compset>
<alias>IHistClm51Sp</alias>
<lname>HIST_DATM%GSWP3v1_CLM51%SP_SICE_SOCN_MOSART_SGLC_SWAV</lname>
</compset>

<compset>
<alias>IHistClm51Bgc</alias>
<lname>HIST_DATM%GSWP3v1_CLM51%BGC_SICE_SOCN_MOSART_SGLC_SWAV</lname>
Expand Down
12 changes: 11 additions & 1 deletion cime_config/testdefs/testlist_clm.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -381,12 +381,13 @@
<option name="wallclock">00:20:00</option>
</options>
</test>
<test name="ERP_D_Ld5" grid="f10_f10_musgs" compset="I2000Clm50Sp" testmods="clm/decStart">
<test name="ERP_D_Ld5" grid="f10_f10_musgs" compset="I2000Clm51Sp" testmods="clm/decStart">
<machines>
<machine name="cheyenne" compiler="intel" category="aux_clm"/>
</machines>
<options>
<option name="wallclock">00:20:00</option>
<option name="comment">2000 Sp test for CLM51</option>
</options>
</test>
<test name="ERP_D_Ld5" grid="f10_f10_musgs" compset="I2000Clm50Sp" testmods="clm/reduceOutput">
Expand Down Expand Up @@ -423,6 +424,15 @@
<option name="wallclock">00:20:00</option>
</options>
</test>
<test name="ERP_D_Ld5" grid="f10_f10_musgs" compset="IHistClm51Sp" testmods="clm/default">
<machines>
<machine name="cheyenne" compiler="intel" category="aux_clm"/>
</machines>
<options>
<option name="wallclock">00:20:00</option>
<option name="comment" >Test Hist compset with Sp for CLM5.1</option>
</options>
</test>
<test name="SMS_Ld5" grid="f09_g17" compset="IHistClm50SpCru" testmods="clm/default">
<machines>
<machine name="cheyenne" compiler="intel" category="ctsm_sci"/>
Expand Down
69 changes: 69 additions & 0 deletions src/biogeophys/CanopyFluxesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -402,6 +402,15 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
real(r8) :: dt_veg_temp(bounds%begp:bounds%endp)
integer :: iv
logical :: is_end_day ! is end of current day
real(r8) :: tl_ini ! leaf temperature from beginning of time step [K]
real(r8) :: ts_ini ! stem temperature from beginning of time step [K]
real(r8) :: cp_leaf !heat capacity of leaves
real(r8) :: dt_stem !change in stem temperature
real(r8) :: frac_rad_abs_by_stem !fraction of incoming radiation absorbed by stems
real(r8) :: lw_stem !internal longwave stem
real(r8) :: lw_leaf !internal longwave leaf
real(r8) :: sa_internal !min(sa_stem,sa_leaf)
real(r8) :: temp

integer :: dummy_to_make_pgi_happy
!------------------------------------------------------------------------------
Expand Down Expand Up @@ -627,6 +636,14 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
obuold(p) = 0._r8
btran(p) = btran0
end do
frac_rad_abs_by_stem = 0._r8
lw_leaf = 0._r8
cp_leaf = 0._r8
lw_stem = 0._r8
sa_internal = 0._r8
dt_stem = 0._r8
tl_ini = 0._r8
ts_ini = 0._r8

! calculate daylength control for Vcmax
do f = 1, fn
Expand Down Expand Up @@ -1049,6 +1066,17 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
dt_veg(p) = (sabv(p) + air(p) + bir(p)*t_veg(p)**4 + &
cir(p)*lw_grnd - efsh - efe(p)) / &
(- 4._r8*bir(p)*t_veg(p)**3 +dc1*wtga +dc2*wtgaq*qsatldT(p))
temp = dt_veg(p)
dt_veg(p) = ((1._r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) &
+ bir(p)*t_veg(p)**4 + cir(p)*lw_grnd) &
- efsh - efe(p) - lw_leaf + lw_stem &
- (cp_leaf/dtime)*(t_veg(p) - tl_ini)) &
/ ((1._r8-frac_rad_abs_by_stem)*(- 4._r8*bir(p)*t_veg(p)**3 &
+ 4._r8*sa_internal*emv(p)*sb*t_veg(p)**3 &
+dc1*wtga+dc2*wtgaq*qsatldT(p))+ cp_leaf/dtime)
if ( abs(temp - dt_veg(p) ) > 1.e-12 )then
call endrun( "Difference greater than roundoff" )
end if
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ekluzek - Thanks for all of your work to demonstrate only roundoff-level diffs. I know from my own experience how painful that can sometimes be.

One question: are the terms you're checking here all order 1 (in magnitude)? If not, or if you're not sure, I generally use relative diff checks for things like this:

if (abs(temp - dt_veg(p)) > 1.e-13_r8 * temp) then

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good point. Actually I'm sure they aren't order 1. dt_veg and err will be order 1, but could also be very small. ulrad and dlrad should be order 100. So doing a relative diff would be a better comparison. If I do that I'm likely to be able to show a relative difference that's a bit smaller as well.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to redo things, though, if you have already done it this way and you have a general sense of the magnitude. My main worry with abs diffs is if the terms are very small to begin with - so an abs diff of 1e-12 could be a relative diff of, say, 1e-6. If you're pretty confident that isn't generally the case here, then I think this is good enough.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I went ahead and did a little testing with this to show it's fine. Relative diff's potentially a bit more troublesome if temp is identically zero, or negative. But with a relative diff I was able to show a tolerance of e-14 was fine for two tests.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Relative diff's potentially a bit more troublesome if temp is identically zero, or negative.

Just a thought for the future: If temp is the old (presumed correct) version, then I usually have success with something like if (abs(temp - dt_veg(p)) > 1.e-13_r8 * temp) – or change the last temp to abs(temp) if negative values are acceptable. This avoids divide by 0. If temp (the old version) was exactly 0 at a point, then I typically expect the new version to be exactly zero, too, and I generally want to ensure that's true (which is done by the above) – since exactly 0 vs. slightly different from 0 can sometimes lead to large downstream differences.

t_veg(p) = tlbef(p) + dt_veg(p)
dels = dt_veg(p)
del(p) = abs(dels)
Expand All @@ -1060,6 +1088,18 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
4._r8*dt_veg(p)) + cir(p)*lw_grnd - &
(efsh + dc1*wtga*dt_veg(p)) - (efe(p) + &
dc2*wtgaq*qsatldT(p)*dt_veg(p))
temp = err(p)
err(p) = (1._r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) &
+ bir(p)*tlbef(p)**3*(tlbef(p) + &
4._r8*dt_veg(p)) + cir(p)*lw_grnd) &
-sa_internal*emv(p)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) &
+ lw_stem &
- (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + &
dc2*wtgaq*qsatldT(p)*dt_veg(p)) &
- (cp_leaf/dtime)*(t_veg(p) - tl_ini)
if ( abs(temp - err(p) ) > 1.e-12 )then
call endrun( "Difference greater than roundoff" )
end if
end if

! Fluxes from leaves to canopy space
Expand Down Expand Up @@ -1186,6 +1226,15 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
err(p) = sabv(p) + air(p) + bir(p)*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) &
!+ cir(p)*t_grnd(c)**4 - eflx_sh_veg(p) - hvap*qflx_evap_veg(p)
+ cir(p)*lw_grnd - eflx_sh_veg(p) - hvap*qflx_evap_veg(p)
temp = err(p)
err(p) = (1.0_r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) + bir(p)*tlbef(p)**3 &
*(tlbef(p) + 4._r8*dt_veg(p)) + cir(p)*lw_grnd) &
- lw_leaf + lw_stem - eflx_sh_veg(p) - hvap*qflx_evap_veg(p) &
- ((t_veg(p)-tl_ini)*cp_leaf/dtime)
if ( abs(temp - err(p) ) > 1.e-11 )then
write(iulog,*) 'Difference = ', temp - err(p)
call endrun( "Difference greater than roundoff" )
end if

! Fluxes from ground to canopy space

Expand Down Expand Up @@ -1271,12 +1320,32 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,

dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(c) + &
emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))
temp = dlrad(p)
dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(c) &
+ emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) &
*(1.0_r8-frac_rad_abs_by_stem) &
+ emv(p)*emg(c)*sb*ts_ini**3*(ts_ini + 4._r8*dt_stem) &
*frac_rad_abs_by_stem
if ( abs(temp - dlrad(p) ) > 1.e-12 )then
call endrun( "Difference greater than roundoff" )
end if

! Upward longwave radiation above the canopy

ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(c) &
+ emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + &
4._r8*dt_veg(p)) + emg(c)*(1._r8-emv(p))*sb*lw_grnd)
temp = ulrad(p)
ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(c) &
+ emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb &
*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))*(1._r8-frac_rad_abs_by_stem) &
+ emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb &
*ts_ini**3*(ts_ini+ 4._r8*dt_stem)*frac_rad_abs_by_stem &
+ emg(c)*(1._r8-emv(p))*sb*lw_grnd)
if ( abs(temp - ulrad(p) ) > 1.e-11 )then
write(iulog,*) 'Difference = ', temp - ulrad(p)
call endrun( "Difference greater than roundoff" )
end if

! Calculate the skin temperature as a weighted sum of all the ground and vegetated fraction
! The weight is the so-called vegetation emissivity, but not that emv is actually an attentuation
Expand Down