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

Modularize snow cover fraction method #769

Merged
merged 47 commits into from Aug 6, 2019

Conversation

billsacks
Copy link
Member

@billsacks billsacks commented Jul 27, 2019

Description of changes

This PR moves the calculation of frac_sno - and the related updates of
snow_depth - into a new set of classes, with one class for each
parameterization (Niu & Yang 2007 and the CLM5 parameterization).

Previously, the code always calculated frac_sno the new way, but then
possibly overwrote it if using the older Niu & Yang method. The new code
cleans this up, only doing the calculations that are needed for each
method.

In addition, other code that is specific to one of the two methods is
now moved to a home that makes this dependence on method explicit. This
includes the addition of newsnow to int_snow: previously, int_snow was
always updated using an equation specific to the newer CLM5
parameterization of frac_sno, which was not appropriate if using the
Niu & Yang 2007 parameterization; this doesn't make a difference
currently, since int_snow is only referenced if using the CLM5
parameterization, but this clears up some confusion. Also, time-constant
parameters read from namelist or the netCDF parameter file now reside in
the appropriate class rather than being more global.

This PR also renames two namelist options to increase clarity:

  • subgridflag is renamed to use_subgrid_fluxes, and is now a logical
  • oldfflag is renamed to snow_cover_fraction_method, and is now a
    string

Specific notes

Contributors other than yourself, if any: @swensosc provided significant
feedback

CTSM Issues Fixed (include github issue #):

Are answers expected to change (and if so in what way)? No

Any User Interface Changes (namelist or namelist defaults changes)?
See notes above

Testing performed, if any:
So far, just these two tests on my mac: both pass and are bit-for-bit:
SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm50BgcCropQianRsGs.bishorn_gnu.clm-default
SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm45BgcCropQianRsGs.bishorn_gnu.clm-oldhyd

These options are incompatible.

Partially addresses ESCOMP#502
Also, fix error checks involving use_subgrid_fluxes in CLMBuildNamelist:
previously, subgridflag was not set at the point when the error checks
involving that value were being done; I have fixed that by always
including use_subgrid_fluxes in the default namelist.

Resolves ESCOMP#502
Testing:
SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm50BgcCropQianRsGs.bishorn_gnu.clm-default
passes and bit-for-bit
This will facilitate upcoming refactoring

Testing:
SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm50BgcCropQianRsGs.bishorn_gnu.clm-default
passes and is bit-for-bit
This will facilitate upcoming refactoring
Make oldhyd inherit from default so that it has some output in a one-day
test
I plan to use this for testing on my laptop, specifically with the test:

SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm45BgcCropQianRsGs.bishorn_gnu.clm-oldhyd
This will be done differently for different snow cover methods, so it
will be part of the snow_cover_fraction polymorphic class.
This will be called from SnowCompaction, too
There was a minor difference between the formulation in SnowCompaction
and the one in FracSnowDuringMeltClm5: the latter has an extra min
statement. It looks like this shouldn't change answers (the min
statement shouldn't be necessary, given the min that is already applied
to smr), but it's possible that it will change answers by roundoff in
rare cases.
This isn't actually used yet in the code, and the two most important
methods are just stubs.
I needed to move it out of SnowHydrologyMod because it will now be used
in SnowCoverFractionMod. I thought about reading it in in the factory
method in SnowCoverFractionMod, but I prefer having this (and other
similar control flags that we might need in the future) in clm_varctl
for a few reasons:

- It avoids introducing a whole new namelist group and a new namelist
  read routine just to hold and read this one variable

- It makes the interface of the factory method more clear, in that we'll
  pass in the method directly

- I like the idea of having these big controls (determining relatively
  big parameterizations) all be in one namelist group for the sake of
  build-namelist: we can ensure that we read this namelist group first,
  then do some other logic based on things set there (see also
  ESCOMP#26 (comment))
We still have to move some code from SnowHydrologyMod into
SnowCoverFractionMod for the N&Y07 method
Move the code out of SnowHydrologyMod, into SnowCoverFractionMod
It's clearer to have the error check more closely associated with the
relevant code.
I have confirmed that, in clm4_0_00, the update of snowdp looked like:

          dz_snowf = qflx_snow_grnd_col(c)/bifall
          snowdp(c) = snowdp(c) + dz_snowf*dtime
Also fix build-namelist error checking of compatibility with
use_subgrid_fluxes
@billsacks billsacks added the PR status: ready PR: author feels this is ready to merge in label Jul 27, 2019
@billsacks billsacks requested a review from ekluzek July 27, 2019 13:19
@billsacks billsacks added this to In progress - master in Upcoming tags Jul 27, 2019
This fixes the error on cheyenne_intel related to trying to allocate an
array (n_melt) that is already allocated.
This tests that adding water isotopes does not change answers

This is not yet passing
1. Make 'base' case have water isotopes on. This way the baselines
   include water isotopes, which is nice for having stronger baseline
   comparisons.

2. For now, use enable_water_tracer_consistency_checks rather than
   enable_water_isotopes. I hope to revert this later, but this seems
   needed for now in order for the test to run to completion.
This covers some code that isn't covered by any existing tests (such as
the oldhyd test), though the amount of additional code coverage is
small, so we don't necessarily need to keep this test long-term.
@billsacks
Copy link
Member Author

With the latest changes, the full test suite is passing & bit-for-bit.

In addition, I compared against master (with changes on master to support
the new tests, including changing oldhyd to inherit from default, and
having an oldhyd_monthly that inherits from monthly rather than default):

  • ERP_D_P36x2_Ld3.f10_f10_musgs.I2000Clm45BgcCrop.cheyenne_gnu.clm-no_subgrid_fluxes:
    bit-for-bit
  • SMS_D_Ld1_P48x1.f10_f10_musgs.I2000Clm45BgcCrop.hobart_nag.clm-oldhyd:
    bit-for-bit
  • SMS_Ly3.f10_f10_musgs.I2000Clm45BgcCrop.cheyenne_intel.clm-oldhyd_monthly:
    bit-for-bit

So I feel this branch is ready after it has been reviewed.

frac_sno(c) = 0._r8
end if
end if ! end of h2osno_total > 0
end do
Copy link
Contributor

@swensosc swensosc Aug 2, 2019

Choose a reason for hiding this comment

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

I wonder if this could be written in such a way to reduce redundancy? Like this maybe:

  do fc = 1, num_c
       c = filter_c(fc)

! initialize frac_sno
    frac_sno(c) = 0._r8
! reset snow_depth if no h2osno    
    if (h2osno_total(c) == 0.0) snow_depth(c) = 0._r8
! update snow depth with new snowfall
   snow_depth(c) = snow_depth(c) + newsnow(c) / bifall(c)
! calculate frac_sno

          if(h2osno_total(c) < 1.0_r8)  then
             frac_sno(c)=min(frac_sno(c),h2osno_total(c))
	  else
	     if(snow_depth(c) > 0.0_r8)  then
             frac_sno(c) = tanh(snow_depth(c) / (2.5_r8 * this%zlnd * &
                  (min(800._r8,(h2osno_total(c)+ newsnow(c))/snow_depth(c))/100._r8)**1._r8) )
             end if
          end if
enddo

one thing I don't like is this though: if (h2osno_total(c) == 0.0) snow_depth(c) = 0._r8
it seems like a function doing the snow_depth reset should be called any time h2osno is
updated, and then it wouldn't need to be done here.

Copy link
Member Author

Choose a reason for hiding this comment

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

Good point. I'll take a stab at something like this. I'm not sure that the code you wrote is exactly right as is (particularly for the initial setting of frac_sno in the h2osno_total < 1 condition), but I see what you're getting at, and I can try for something similar.

The h2osno_total < 1 block makes this messier, and in looking more closely, I'm not sure it's really correct: although I don't understand what it's trying to do, it feels to me like this should be looking at (h2osno_total+newsnow), both in the condition and in the min, and that a similar min should be applied in the case where h2osno_total==0 but newsnow > 0. But I'm reluctant to make that change without really understanding it, so unless you give me guidance otherwise, I'll plan to keep the current code logic when doing the restructuring, even if this makes things less clean and consistent than they could be.

I agree that it feels like the snow_depth reset should happen somewhere else, but determining all the places that would need to happen would take some time; for now, I'll probably just add a note to look back at that.

Copy link
Member Author

Choose a reason for hiding this comment

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

Fixed in 24ffe4e. I have run a 3-year test (SMS_Ly3.f10_f10_musgs.I2000Clm45BgcCrop.cheyenne_intel.clm-oldhyd_monthly) and verified that it is bit-for-bit with master (despite the theoretical potential for this commit to change answers, as noted in the log message for the commit).

Copy link
Member Author

Choose a reason for hiding this comment

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

@swensosc it might be good if you could look at the end result of 24ffe4e to make sure it (a) looks correct, and (b) accomplishes what you want.

! update fsca by new snow event, add to previous fsca
if (newsnow(c) > 0._r8) then
fsno_new = 1._r8 - (1._r8 - tanh(this%accum_factor * newsnow(c))) * (1._r8 - frac_sno(c))
frac_sno(c) = fsno_new
Copy link
Contributor

Choose a reason for hiding this comment

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

is there a benefit to keeeping this fsno_new variable?

Copy link
Member Author

Choose a reason for hiding this comment

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

No benefit; I'll remove it: thanks.

! initialize frac_sno and snow_depth when no snow present initially
if (newsnow(c) > 0._r8) then
z_avg = newsnow(c)/bifall(c)
fmelt=newsnow(c)
Copy link
Contributor

Choose a reason for hiding this comment

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

can fmelt be removed?

Copy link
Member Author

Choose a reason for hiding this comment

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

Good point; I'll remove it.

@@ -0,0 +1,785 @@
module SnowCoverFractionMod
Copy link
Contributor

Choose a reason for hiding this comment

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

It feels like most of this module is infrastructure. Is there any way to make the 'science' code more front-and-center?

Copy link
Member Author

Choose a reason for hiding this comment

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

I have a few ideas, and would be interested to hear what you think of them, or if you have other thoughts:

  1. I could move the definition of the abstract snow_cover_fraction_type and the whole abstract interface block near the top of this module into its own module. This is what I have done for other object-oriented parts of the code; the downside is just having more files. That would remove about 85 lines from the start of the file.

  2. I could also move the factory method (CreateAndInitScfMethod) into its own module. Again, this is what I have done for other object-oriented parts of the code; the downside is again having more files. (Unfortunately, this cannot be combined with the module for the abstract type, because that introduces a circular dependency.) This would remove about 50 lines of code. Alternatively, I could move this method to the bottom of the file.

  3. I could move the various science routines to above the infrastructurey routines that are just used in initialization. If I do this, I'd like to know if you'd prefer:

(a) The different parameterizations each get their own file

(b) The different parameterizations stay together in a single file, with all subroutines for a given parameterization grouped together (as they are now), but just change the order of these routines within each grouping so that the science routines come first.

(c) The different parameterizations stay together in a single file, with two sets of groupings: first we have all of the science routines for NY07 followed by all of the science routines for SL12; then we have all of the initialization routines for NY07 followed by all of the initialization routines for SL12.

(d) The different parameterizations stay together in a single file, with the groupings flipped around: first we have UpdateSnowDepthAndFracNY07 followed by UpdateSnowDepthAndFracSL12; then we have AddNewsnowToIntsnowNY07 followed by AddNewsnowToIntsnowSL12; etc.

I am totally happy to do any of these, so please let me know what you'd find most helpful.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think personally I like the idea of an abstract interface for the infrastructure, and then a ny07 implementation and clm5 (still not sure I like calling it the clm5 version) implementation. I think it's OK to have the different implementations in different files, and maybe better actually. The science is then pretty much all front and center in either of the implementations. I picture it can be useful to compare the different implementations, but then having them in different files is actually helpful. Since, the implementations are quite different it seems unlikely that you'd have shared code. And if you do you could probably put in in a base class.

Copy link
Member Author

Choose a reason for hiding this comment

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

Thanks for your thoughts @ekluzek . You raise good points about the advantages of having separate files. So if I understand correctly, it sounds like your preferences on (3) would be (a) (which I think would also imply doing (1) & (2)) then (b). I'll talk more with @swensosc about this next week to see if he agrees.

Copy link
Contributor

@ekluzek ekluzek left a comment

Choose a reason for hiding this comment

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

I have a few changes that I ask for. The refactoring is much better than before, and overall I'm very happy with it though so I could be convinced to let some of the requests slide. For example, I ask if "clm5" is the best name for the new version. You might have already thought about this and didn't come up with anything better, and if that's the case I'm fine with that. I wanted to make sure you had thought about it though. When we have clm5_5, and clm6_0 as options, it'll be confusing to refer to something as "clm5".

#######################################################################
# namelist group: scf_clm5_inparm #
#######################################################################
setup_logic_scf_clm5($opts, $nl_flags, $definition, $defaults, $nl);
Copy link
Contributor

Choose a reason for hiding this comment

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

Should "clm5" really be in the name here?

Copy link
Member Author

Choose a reason for hiding this comment

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

Sorry that I wasn't more clear about this: I used "clm5" somewhat as a placeholder just because that was the tentative name used in the discussion in #503 , and I was waiting for a better name. @swensosc has now suggested the name "sl12" (for Swenson & Lawrence 2012), and I will switch the references to clm5 to sl12 throughout this PR.

@@ -0,0 +1,44 @@
"""
Copy link
Contributor

Choose a reason for hiding this comment

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

Just checking, this seems unrelated to the snow fraction refactor. I'm OK with having it in place, I'm just checking to make sure I understand.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, this is unrelated.

call endrun('N&Y07 snow cover fraction parameterization incompatible with use_subgrid_fluxes')
end if

call this%ReadParamsNY07( &
Copy link
Contributor

Choose a reason for hiding this comment

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

I think you are right about having readparams called from the init method rather than readParams. Stefan had setup the readParams configuration, and I have been keeping that design. But, it makes more sense that there will be some sort of Init method and the local init method will know if it should read parameters or not.

@@ -0,0 +1,785 @@
module SnowCoverFractionMod
Copy link
Contributor

Choose a reason for hiding this comment

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

I think personally I like the idea of an abstract interface for the infrastructure, and then a ny07 implementation and clm5 (still not sure I like calling it the clm5 version) implementation. I think it's OK to have the different implementations in different files, and maybe better actually. The science is then pretty much all front and center in either of the implementations. I picture it can be useful to compare the different implementations, but then having them in different files is actually helpful. Since, the implementations are quite different it seems unlikely that you'd have shared code. And if you do you could probably put in in a base class.

type(file_desc_t), intent(inout) :: params_ncid ! pio netCDF file id for parameter file
end subroutine Init_Interface

subroutine UpdateSnowDepthAndFrac_Interface(this, bounds, num_c, filter_c, &
Copy link
Contributor

Choose a reason for hiding this comment

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

Personally I still like to have at least a short description of the subroutine

Copy link
Member Author

Choose a reason for hiding this comment

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

Done in 3daf306

if (newsnow(c) > 0._r8) then
snow_depth(c) = newsnow(c) / bifall(c)
if(snow_depth(c) > 0.0_r8) then
frac_sno(c) = tanh(snow_depth(c) / (2.5_r8 * this%zlnd * &
Copy link
Contributor

Choose a reason for hiding this comment

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

There are a couple magic numbers that could be made parameters here.

Copy link
Member Author

Choose a reason for hiding this comment

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

I'd rely on scientists (@swensosc or others) to provide meaningful names for the parameters.

Resolved Conflicts:
	cime_config/testdefs/testlist_clm.xml
	src/main/controlMod.F90
This has the theoretical potential to change answers under the following
conditions:

- Now frac_sno is being explicitly reset to 0; before, under some
  conditions (h2osno_total > 0, but snow_depth <= 0), frac_sno would
  remain at its old value.

- If h2osno_total or newsnow are ever negative, we could have different
  answers

A one-day
test (SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm45BgcCropQianRsGs.bishorn_gnu.clm-oldhyd)
is still bit-for-bit; I'll run a 3-year test next.
This is from Swenson & Lawrence 2012
@swensosc
Copy link
Contributor

swensosc commented Aug 5, 2019 via email

Sean Swenson and others suggested using these longer names rather than
the more cryptic short abbreviations. Part of the motivation is once we
have 4 separate files, this will make the science modules stand out more
because of their longer names.

I wanted to name the namelist group scf_SwensonLawrence2012_inparm for
consistency with the namelist option, but this didn't work with the
namelist generation scripts, presumably because of the capitalization;
so I'm using the all-lowercase form, scf_swenson_lawrence_2012_inparm.
Since the init method is only called from the factory method, it's
fairly straightforward to allow the different init methods to have
different interfaces, rather than having a bunch of unused arguments in
some of the methods. The big motivation here is avoiding the need to
change all parameterizations' Init methods if you add a new
parameterization with a new input requirement.

This ended up being a little non-trivial because of the need for the
select type statement. I think we could have avoided that if we did all
initialization in the constructor (avoiding an Init method
entirely). But we have had some compilers (e.g., intel17) that don't
seem to work right with function constructors, so for now we're sticking
with an Init subroutine.
Sean Swenson pointed out that it would be better to put the science up
front.
@billsacks billsacks merged commit 6c91dff into ESCOMP:master Aug 6, 2019
Upcoming tags automation moved this from In progress - master to Master Tags/Issues Done Aug 6, 2019
@billsacks billsacks deleted the oldfflag_cleanup branch August 6, 2019 20:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
PR status: ready PR: author feels this is ready to merge in
Projects
Upcoming tags
Done (non release/external)
Development

Successfully merging this pull request may close these issues.

Clean up CanopyHydrology oldfflag can NOT be used with subgridflag==1, and rename subgridflag
3 participants