Bug fixes for energy imbalance associated with surface water and lakes #307

olyson · 2018-02-27T16:43:17Z

These are bug fixes for the land energy imbalance over land as determined by coupler diagnostics. They include bug fixes for surface water phase change and lake/snow interactions developed by Sean Swenson and Keith Oleson.
The following modules have been modified:
src/biogeophys/EnergyFluxType.F90
src/biogeophys/SoilFluxesMod.F90
src/biogeophys/SoilTemperatureMod.F90
src/biogeophys/SurfaceAlbedoMod.F90
src/biogeophys/SurfaceRadiationMod.F90
Diagnostics compared to an offline control are here:
http://webext.cgd.ucar.edu/I1850/clm50_r272_1deg_GSWP3V1_iso_h2osfclakefix_1850/lnd/clm50_r272_1deg_GSWP3V1_iso_h2osfclakefix_1850.1443_1461-clm50_r272_1deg_GSWP3V1_iso_1850.1443_1461/setsIndex.html

billsacks · 2018-02-27T19:25:50Z

@olyson - I know the distinction between assignee and reviewer is a bit unclear. Typically it works best to have a single assignee and multiple reviewers.

billsacks

See a few inline comments

billsacks · 2018-02-27T21:27:13Z

src/biogeophys/EnergyFluxType.F90

@@ -21,6 +21,9 @@ module EnergyFluxType
  type, public :: energyflux_type

     ! Fluxes
+!scs


Need to remove comment lines like this, and other similar comments (scs comments and other commented-out lines).

@swensosc and @olyson - for the future, it makes code reviews easier if you can remove these commented-out things and commented initials before submitting the PR.

billsacks · 2018-02-27T21:34:09Z

src/biogeophys/SurfaceRadiationMod.F90

-             if (subgridflag == 0) then 
+!             if (subgridflag == 0) then
+!scs
+             if (subgridflag == 0 .or. lun%itype(l) == istdlak) then


@swensosc There is one other occurrence of subgridflag, on line 763 of this module. Should that one also check if the landunit type is istdlak?

hi, bill, I think it is okay to add the 'or lake' condition there as well. Then it will be consistent. thanks, sean
…

On Tue, Feb 27, 2018 at 2:38 PM, Bill Sacks ***@***.***> wrote: In src/biogeophys/SurfaceRadiationMod.F90 <#307 (comment)>: > @@ -664,7 +666,9 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & > sabg_soil(p) = sabg(p) > endif > ! if no subgrid fluxes, make sure to set both components equal to weighted average > - if (subgridflag == 0) then > +! if (subgridflag == 0) then > +!scs > + if (subgridflag == 0 .or. lun%itype(l) == istdlak) then @swensosc <https://github.com/swensosc> There is one other occurrence of subgridflag, on line 763 of this module. Should that one also check if the landunit type is istdlak? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://github.com/ESCOMP/ctsm/pull/307/files/248307d5bd096a23ff6f85efaf255c0254e0e348#r171073670>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AaLe_HJr-VWfkTh0OQDdWZvuwvJdKwo8ks5tZHW9gaJpZM4SVQzI> .

I've pushed this change to the branch

billsacks · 2018-02-27T21:37:44Z

src/biogeophys/SoilTemperatureMod.F90


            ! mass of water converted from liquid to ice
            xm(c) = hm(c)*dtime/hfus  
            temp1 = h2osfc(c) + xm(c)

-            ! compute change in cv due to additional ice
-            dcv(c)=cpice*min(abs(xm(c)),h2osfc(c))


Should also remove the declaration of this no-longer-used local variable, dcv

billsacks · 2018-02-27T21:55:30Z

src/biogeophys/SoilTemperatureMod.F90

+                       /(c1 + c2)             
+
+                  eflx_h2osfc_to_snow_col(c) =(t_h2osfc(c)-t_soisno(c,0))*c2/dtime
+               else


It looks like this else block is exactly the same as the above else if block, except that here there is no dhsdT in the calculation of c1. To reduce this code duplication, can you do something like this?:

if (snl(c) == 0) then !initialize for next time step t_soisno(c,0) = t_h2osfc(c) eflx_h2osfc_to_snow_col(c) = 0. else if (snl(c) == -1) then somevar = dhsdT(c) else somevar = 0._r8 end if c1=frac_sno(c)*(dtime/fact(c,0) - somevar*dtime) ! Rest of the code here

swensosc · 2018-02-27T22:04:42Z

or: c1=frac_sno(c)*(dtime/fact(c,0) if (snl(c) == -1) c1 = c1 - frac_sno(c)*dhsdT(c)*dtime then no need for additional variable (?)

…

On Tue, Feb 27, 2018 at 3:00 PM, Bill Sacks ***@***.***> wrote: ***@***.**** requested changes on this pull request. See a few inline comments ------------------------------ In src/biogeophys/EnergyFluxType.F90 <#307 (comment)>: > @@ -21,6 +21,9 @@ module EnergyFluxType type, public :: energyflux_type ! Fluxes +!scs Need to remove comment lines like this, and other similar comments (scs comments and other commented-out lines). @swensosc <https://github.com/swensosc> and @olyson <https://github.com/olyson> - for the future, it makes code reviews easier if you can remove these commented-out things and commented initials before submitting the PR. ------------------------------ In src/biogeophys/SurfaceRadiationMod.F90 <#307 (comment)>: > @@ -664,7 +666,9 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & sabg_soil(p) = sabg(p) endif ! if no subgrid fluxes, make sure to set both components equal to weighted average - if (subgridflag == 0) then +! if (subgridflag == 0) then +!scs + if (subgridflag == 0 .or. lun%itype(l) == istdlak) then @swensosc <https://github.com/swensosc> There is one other occurrence of subgridflag, on line 763 of this module. Should that one also check if the landunit type is istdlak? ------------------------------ In src/biogeophys/SoilTemperatureMod.F90 <#307 (comment)>: > ! mass of water converted from liquid to ice xm(c) = hm(c)*dtime/hfus temp1 = h2osfc(c) + xm(c) - ! compute change in cv due to additional ice - dcv(c)=cpice*min(abs(xm(c)),h2osfc(c)) Should also remove the declaration of this no-longer-used local variable, dcv ------------------------------ In src/biogeophys/SoilTemperatureMod.F90 <#307 (comment)>: > + !initialize for next time step + t_soisno(c,0) = t_h2osfc(c) + eflx_h2osfc_to_snow_col(c) = 0. + else if (snl(c) == -1) then + c1=frac_sno(c)*(dtime/fact(c,0) - dhsdT(c)*dtime) + if ( frac_h2osfc(c) /= 0.0_r8 )then +! c2=frac_h2osfc(c)*(-cpliq*xm(c)) + c2=(-cpliq*xm(c)) + else + c2=0.0_r8 + end if + t_soisno(c,0) = (c1*t_soisno(c,0)+ c2*t_h2osfc(c)) & + /(c1 + c2) + + eflx_h2osfc_to_snow_col(c) =(t_h2osfc(c)-t_soisno(c,0))*c2/dtime + else It looks like this else block is exactly the same as the above else if block, except that here there is no dhsdT in the calculation of c1. To reduce this code duplication, can you do something like this?: if (snl(c) == 0) then !initialize for next time step t_soisno(c,0) = t_h2osfc(c) eflx_h2osfc_to_snow_col(c) = 0. else if (snl(c) == -1) then somevar = dhsdT(c) else somevar = 0._r8 end if c1=frac_sno(c)*(dtime/fact(c,0) - somevar*dtime) ! Rest of the code here — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#307 (review)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AaLe_LdVsmyxKzoj7fzOp1193VJ2ZKDJks5tZHsRgaJpZM4SVQzI> .

billsacks · 2018-02-27T22:22:03Z

or:

c1=frac_sno(c)*(dtime/fact(c,0)
if (snl(c) == -1) c1 = c1 - frac_sno(c)*dhsdT(c)*dtime

then no need for additional variable (?)

Yes, that's fine, too

Add some land ice diagnostic vars needed for CMIP6 Add some diagnostic variables needed for analyzing land ice that have been requested by some of the MIPs in CMIP6 (especially ISMIP). Also, fixes c2l_scale_type to fix urban scaling for SNOWICE, SNOWLIQ Some changes from Leo van Kampenhout.

ekluzek · 2018-03-02T19:07:46Z

src/biogeophys/SoilTemperatureMod.F90

+                       /(c1 + c2)             
+
+                  eflx_h2osfc_to_snow_col(c) =(t_h2osfc(c)-t_soisno(c,0))*c2/dtime
+               else


Also note that there is a division by "fact" that isn't being checked for. Can fact be small or zero?

OK, I did more looking at the code. It looks like based on ComputeHeatDiffFluxAndFactor, fact shouldn't ever be zero. But, it could possibly be small. We possibly should do an assert check in ComputeHeatDiffFluxAndFactor that fact isn't identically zero (which would only be done in DEBUG mode). But, I also wonder about checking for it being small.

ekluzek · 2018-03-02T19:13:39Z

src/biogeophys/SoilTemperatureMod.F90

+                        else
+                           hm(c,j) = frac_sno_eff(c)*(dhsdT(c)*tinc(c,j) - tinc(c,j)/fact(c,j))
+                        endif
+!scs



More divides by fact without checking.

See the comment above.

ekluzek · 2018-03-02T19:16:31Z

src/biogeophys/SoilTemperatureMod.F90

+               else
+                  cv(c,j) = 1.0e-6_r8
+               endif
+!scs


Probably should check for frac_sno > small value rather than greater than zero.

You are probably right but I think Sean should weigh in on this.

Yes, that's fine.

ekluzek · 2018-03-03T00:06:28Z

The following change allows the model to get past the divide by zero issue:

diff --git a/src/biogeophys/SoilTemperatureMod.F90 b/src/biogeophys/SoilTemperatureMod.F90
index ef650abb..582a0312 100644
--- a/src/biogeophys/SoilTemperatureMod.F90
+++ b/src/biogeophys/SoilTemperatureMod.F90
@@ -342,7 +342,7 @@ subroutine SoilTemperature(bounds, num_urbanl, filter_urbanl, num_nolakec, filte
 !         dz_h2osfc(c) = max(1.0e-6_r8,1.0e-3*h2osfc(c))
 !         c_h2osfc(c)  = cpliq*denh2o*dz_h2osfc(c) !"areametric" heat capacity [J/K/m^2]
 !scs
-         if (h2osfc(c) > 0._r8) then 
+         if ( (h2osfc(c) > 1.0e-6_r8) .and. (frac_h2osfc(c) > 1.e-6_r8) ) then 
             c_h2osfc(c)  = max(1.0e-6_r8,cpliq*h2osfc(c)/frac_h2osfc(c))
             dz_h2osfc(c) = max(1.0e-6_r8,1.0e-3*h2osfc(c)/frac_h2osfc(c))
          else

swensosc · 2018-03-05T15:13:11Z

I think fact will always be nonzero. It is basically the soil moisture, multiplied by a couple of constants.

…

On Fri, Mar 2, 2018 at 12:07 PM, Erik Kluzek ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In src/biogeophys/SoilTemperatureMod.F90 <#307 (comment)>: > + !initialize for next time step + t_soisno(c,0) = t_h2osfc(c) + eflx_h2osfc_to_snow_col(c) = 0. + else if (snl(c) == -1) then + c1=frac_sno(c)*(dtime/fact(c,0) - dhsdT(c)*dtime) + if ( frac_h2osfc(c) /= 0.0_r8 )then +! c2=frac_h2osfc(c)*(-cpliq*xm(c)) + c2=(-cpliq*xm(c) - frac_h2osfc(c)*dhsdT(c)*dtime) + else + c2=0.0_r8 + end if + t_soisno(c,0) = (c1*t_soisno(c,0)+ c2*t_h2osfc(c)) & + /(c1 + c2) + + eflx_h2osfc_to_snow_col(c) =(t_h2osfc(c)-t_soisno(c,0))*c2/dtime + else Also note that there is a division by "fact" that isn't being checked for. Can fact be small or zero? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#307 (review)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AaLe_Bb3FvT6h7c9HSUjXF0sXAeS7qVnks5taZiEgaJpZM4SVQzI> .

ekluzek · 2018-03-05T23:01:11Z

The change I made to get past the divide by zero, does change answers. So I will need to rerun testing at minimum.

olyson · 2018-03-06T00:14:51Z

The divide by zero fix looks ok to me, Erik.

…into energy_imbalance_fix

…eclaration of variables that aren't used

billsacks · 2018-03-07T13:50:45Z

Thanks for doing this cleanup, @ekluzek . It looks like most of my comments are addressed now; the one remaining one is #307 (comment) . However, given the time-criticality of this, I'll go ahead and mark this as accepted - and if it doesn't get resolved this time around, we can come back to it later.

…into energy_imbalance_fix

ekluzek · 2018-03-02T19:49:07Z

src/biogeophys/SoilTemperatureMod.F90

+                       /(c1 + c2)             
+
+                  eflx_h2osfc_to_snow_col(c) =(t_h2osfc(c)-t_soisno(c,0))*c2/dtime
+               else


OK, I did more looking at the code. It looks like based on ComputeHeatDiffFluxAndFactor, fact shouldn't ever be zero. But, it could possibly be small. We possibly should do an assert check in ComputeHeatDiffFluxAndFactor that fact isn't identically zero (which would only be done in DEBUG mode). But, I also wonder about checking for it being small.

ekluzek · 2018-03-02T19:49:40Z

src/biogeophys/SoilTemperatureMod.F90

+                        else
+                           hm(c,j) = frac_sno_eff(c)*(dhsdT(c)*tinc(c,j) - tinc(c,j)/fact(c,j))
+                        endif
+!scs



See the comment above.

ekluzek · 2018-03-08T22:37:38Z

I'm having a bunch of problems fully completing testing on cheyenne. But, I'm only see expected fails, and testing is complete on hobart. So going to bring this to master and complete testing afterwards.

Incorporates bug fixes from CLM5 associated with surface water and lakes. The bug fixes incoporated here includes: - Bug fix incoporated in clm4_5_1_r087, and - ESCOMP/CTSM#307 [non-BFB]

Incorporates bug fixes from CLM5 associated with surface water and lakes. The bug fixes incoporated here from CLM include: - clm4_5_1_r087, and - ESCOMP/CTSM#307 [non-BFB]

This is based on bishtgautam/lnd/radiation-fix, with CPP directive APPLY_POST_DECK_BUGFIXES added to allow using CONFIG options to turn on/off the fixes. It incorporates bug fixes from CLM5 associated with surface water and lakes. The bug fixes incoporated here from CLM include: - clm4_5_1_r087, and - ESCOMP/CTSM#307 modified: components/clm/src/biogeophys/BalanceCheckMod.F90 modified: components/clm/src/biogeophys/CanopyHydrologyMod.F90 modified: components/clm/src/biogeophys/EnergyFluxType.F90 modified: components/clm/src/biogeophys/SnowHydrologyMod.F90 modified: components/clm/src/biogeophys/SoilFluxesMod.F90 modified: components/clm/src/biogeophys/SoilTemperatureMod.F90 modified: components/clm/src/biogeophys/SurfaceAlbedoMod.F90 modified: components/clm/src/biogeophys/SurfaceRadiationMod.F90 [non-BFB]

Incorporates bug fixes from CLM5 associated with surface water and lakes. The bug fixes incorporated here are from CLM include: - clm4_5_1_r087, and - ESCOMP/CTSM#307 The bug fixes are added with CPP directive APPLY_POST_DECK_BUGFIXES [BFB]

Bug fixes for energy imbalance associated with surface water and lakes

Bug fixes for energy imbalance associated with surface water and lakes

248307d

olyson added tag: bug - critical big problems in important configurations priority: high High priority task to fix soon, e.g., because it is a problem in important configurations labels Feb 27, 2018

olyson assigned billsacks, ekluzek, dlawrenncar and swensosc Feb 27, 2018

olyson requested a review from ekluzek February 27, 2018 16:43

ekluzek requested review from billsacks, dlawrenncar and swensosc February 27, 2018 19:37

ekluzek unassigned billsacks, dlawrenncar and swensosc Feb 27, 2018

ekluzek added this to the cesm2 milestone Feb 27, 2018

billsacks requested changes Feb 27, 2018

View reviewed changes

olyson and others added 4 commits February 28, 2018 12:52

Another bug fix from Sean and also no snowfall on h2osfc

de46c5d

Add lake landunit check to SurfaceRadiationMod.F90

4a390a7

Fix copy and paste error from Sean, that was setting the wrong variable

41afa17

ekluzek reviewed Mar 2, 2018

View reviewed changes

ekluzek added 4 commits March 5, 2018 22:11

Prevent a divide by zero in new code added in

fd97f39

Merge branch 'energy_imbalance_fix' of https://github.com//olyson/ctsm …

594318b

…into energy_imbalance_fix

Get rid of scs comments and commented out code

96872f2

Get rid of more scs comments, add comments to new variables, remove d…

69cf2c1

…eclaration of variables that aren't used

billsacks approved these changes Mar 7, 2018

View reviewed changes

ekluzek added 3 commits March 7, 2018 11:34

Make 1.e-6 a parameter and fix some formating alignment

41c79f0

Reduce the duplication for c2, eflx_h2osfc_to_snow_col, and t_soisno

0084183

Merge branch 'energy_imbalance_fix' of https://github.com//olyson/ctsm …

e14431d

…into energy_imbalance_fix

ekluzek approved these changes Mar 8, 2018

View reviewed changes

ekluzek merged commit 97cdc7d into ESCOMP:master Mar 8, 2018

ekluzek mentioned this pull request Mar 14, 2018

Energy imbalance over land per coupler diagnostics #304

Closed

wlin7 mentioned this pull request Mar 28, 2018

Fixes land energy imbalance and applies post deck CPP directive E3SM-Project/E3SM#2204

Merged

billsacks pushed a commit to billsacks/ctsm that referenced this pull request Feb 22, 2019

Merge pull request ESCOMP#307 from olyson/energy_imbalance_fix

4bc1bc6

Bug fixes for energy imbalance associated with surface water and lakes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bug fixes for energy imbalance associated with surface water and lakes #307

Bug fixes for energy imbalance associated with surface water and lakes #307

olyson commented Feb 27, 2018

billsacks commented Feb 27, 2018

billsacks left a comment

billsacks Feb 27, 2018

billsacks Feb 27, 2018

swensosc Feb 27, 2018 via email

olyson Feb 28, 2018

billsacks Feb 27, 2018

billsacks Feb 27, 2018

swensosc commented Feb 27, 2018 via email

billsacks commented Feb 27, 2018

ekluzek Mar 2, 2018

ekluzek Mar 2, 2018

ekluzek Mar 2, 2018

ekluzek Mar 2, 2018

ekluzek Mar 2, 2018

olyson Mar 6, 2018

swensosc Mar 6, 2018

ekluzek commented Mar 3, 2018

swensosc commented Mar 5, 2018 via email

ekluzek commented Mar 5, 2018

olyson commented Mar 6, 2018

billsacks commented Mar 7, 2018

ekluzek Mar 2, 2018

ekluzek Mar 2, 2018

ekluzek commented Mar 8, 2018

Bug fixes for energy imbalance associated with surface water and lakes #307

Bug fixes for energy imbalance associated with surface water and lakes #307

Conversation

olyson commented Feb 27, 2018

billsacks commented Feb 27, 2018

billsacks left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

swensosc Feb 27, 2018 via email

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

swensosc commented Feb 27, 2018 via email

billsacks commented Feb 27, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ekluzek commented Mar 3, 2018

swensosc commented Mar 5, 2018 via email

ekluzek commented Mar 5, 2018

olyson commented Mar 6, 2018

billsacks commented Mar 7, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ekluzek commented Mar 8, 2018