Allow collapsing to N dominant PFTs

[billsacks]

For the initial NWP configuration, we're pointing to a surface dataset that @barlage made manually. Eventually we want to be able to support this configuration out-of-the-box, both so (1) users can easily create their own datasets and set up runs equivalent to this at their own resolution, and (2) we can easily recreate this surface dataset for testing purposes as input requirements change.

One key aspect of this is allowing the model to just run with the dominant N PFTs in each grid cell (where N may often be 1). This could be done either at runtime or at mksurfdata_map time. Based on discussion yesterday, it seemed like doing this particular piece at runtime might be easiest and most useful, but either way could be done (note that, for NWP applications, it will be typical for users to create their own surface datasets for their own particular grid, rather than relying on out-of-the-box surface datasets -- so setting these kinds of things at mksurfdata_map time becomes more reasonable).

Other things we may want for this are:

• Using higher thresholds for whether special landunits are present. For example, maybe lake is only present if it takes up 10% or more of the grid cell? We may want these thresholds to be selectable via the mksurfdata_map namelist (currently I think they're hardcoded in mksurfdata_map). (Update 2018-11-29: moved to Allow zeroing out special landunits if their area is below some threshold #581 .)

• Turning off some things that take a long time to generate - such as the surface topographic variability and glacier elevation classes. We could have options to just use constant, pre-defined values for these. (Update 2018-11-29: moved to Allow turning off pieces of mksurfdata_esmf that take a long time to generate #582 .)

[billsacks]

Initial target for this is ability to say, at runtime, "I want to just run with the dominant N PFT(s) in each grid cell" (rather than all PFTs). We're not sure if we also want the bare land PFT in each grid cell. Some old CLM4 (or earlier) code used to do something like this (I think in SubgridMod??).

@dlawrenncar asks how this would combine with collapsing crops from 78 to 16? Our initial thought is we would do the crop collapsing first, then pick the dominant N PFTs, treating all crop as one.

Then this could be extended to turn off special landunits if their area is below some threshold.

There may also be subtleties with how this works in transient mode. We might want to start by not accommodating transient.

[slevis-lmwg]

Gathering additional details.

Received by email from @dlawrenncar:
"...we want the ability to reduce to maximum of 2 or 3 PFTs for the purpose of faster simulations for NWP (or any other time we want to increase speed). Seems like there could actually be two settings. One for the # of PFTs and one for minimum PFT weight. E.g., could set max_PFTs = 3 and min_PFT_weight to 5% or something like that. Not sure how to deal with crops under this scenario. Maybe leave the CFTs as they are now so that this only affects natural vegetation."

Received by email from @billsacks:
"I think changes are needed in subgridMod and initGridCellsMod, and possibly elsewhere."

[slevis-lmwg]

Notes that should not clutter this page, I will gather here:
https://docs.google.com/document/d/1XZUXk-1VrvV6Rvl_cQPG9slYXiO24_4ujaiJGoI-ixI/edit

[billsacks]

I'm renaming this issue to be specific to allowing N dominant PFTs; I'll open separate issues for the other desired features.

[billsacks]

@slevisconsulting - I'm responding to some of your initial comments in the google doc here (so that these thoughts are easier for me to find later):

1. My intuition tells me to start with subroutine surfrd_veg_all in surfrdMod.F90.

Right after…
call collapse_crop_types (also called in subroutine dyncrop_interp for transient runs)
add
call collapse_all_pfts (I doubt that NWP will perform transient runs)
...to perform a similar “collapse” for the pfts that are left.

... some text omitted ...

2. Bill expects changes in subgridMod, initGridCellsMod, and possibly elsewhere.

I like your thought. I had been stuck in the mode of thinking about what (I think) was in place a number of years ago, but your suggestion is better. With your suggestion, it's likely that no changes will be needed in subgridMod and initGridCellsMod.

The new subroutine will need to return updated
wt_lunit(begg:endg,istsoil)
wt_lunit(begg:endg,istcrop)
wt_nat_patch(begg:endg,:)
wt_cft(begg:endg,:)
fert_cft(begg:endg,:)

I have been a little uneasy with the fact that fert_cft is updated in this same collapse routine: This design makes it hard to add additional fields that should be collapsed, and makes the current collapse routine harder to understand because it's doing multiple things. This is not just a theoretical issue: we'll need to do this for irrigation method (see #565) (this will use a slightly different method since it's a discrete rather than continuous field, but it will follow similar logic). I haven't given this much thought, but what I think would be cleanest is if the main collapse routine returns the information needed to collapse fert_cft and any similar field, then there is a separate routine that takes that information and applies it to any field – so you can call that separate routine on fert_cft or any other field that needs similar collapsing. However, I recognize that that would require more rework, so I'm okay if you don't want to tackle that change right now, and we can come back to it later, if you prefer.

I doubt that NWP will perform transient runs

Yeah, at least for now let's assume that this new collapsing cannot be combined with transient runs (and add an error-check preventing that), because it seems like it will be a little messy to allow their combination - possibly requiring some rearrangement of dyn_subgrid code.

[dlawrenncar]

I guess it is fine to assume that collapsing cannot be combined with transient runs for now, but I was hoping that this might be a means to speed up CTSM-climate configuration as well by reducing the max number of PFTs represented in any grid cell, which would require having the capability for this to work with transient runs as well.

…

On Thu, Nov 29, 2018 at 3:18 PM Bill Sacks ***@***.***> wrote: @slevisconsulting <https://github.com/slevisconsulting> - I'm responding to some of your initial comments in the google doc here (so that these thoughts are easier for me to find later): 1. My intuition tells me to start with subroutine surfrd_veg_all in surfrdMod.F90. Right after… call collapse_crop_types (also called in subroutine dyncrop_interp for transient runs) add call collapse_all_pfts (I doubt that NWP will perform transient runs) ...to perform a similar “collapse” for the pfts that are left. ... some text omitted ... 1. Bill expects changes in subgridMod, initGridCellsMod, and possibly elsewhere. I like your thought. I had been stuck in the mode of thinking about what (I think) was in place a number of years ago, but your suggestion is better. With your suggestion, it's likely that no changes will be needed in subgridMod and initGridCellsMod. The new subroutine will need to return updated wt_lunit(begg:endg,istsoil) wt_lunit(begg:endg,istcrop) wt_nat_patch(begg:endg,:) wt_cft(begg:endg,:) fert_cft(begg:endg,:) I have been a little uneasy with the fact that fert_cft is updated in this same collapse routine: This design makes it hard to add additional fields that should be collapsed, and makes the current collapse routine harder to understand because it's doing multiple things. This is not just a theoretical issue: we'll need to do this for irrigation method (see #565 <#565>) (this will use a slightly different method since it's a discrete rather than continuous field, but it will follow similar logic). I haven't given this much thought, but what I think would be cleanest is if the main collapse routine returns the information needed to collapse fert_cft and any similar field, then there is a separate routine that takes that information and applies it to any field – so you can call that separate routine on fert_cft or any other field that needs similar collapsing. However, I recognize that that would require more rework, so I'm okay if you don't want to tackle that change right now, and we can come back to it later, if you prefer. I doubt that NWP will perform transient runs Yeah, at least for now let's assume that this new collapsing cannot be combined with transient runs (and add an error-check preventing that), because it seems like it will be a little messy to allow their combination - possibly requiring some rearrangement of dyn_subgrid code. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#457 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AUAcVEuAELGoiPb0h0VqYgXgBUmMLb-Qks5u0F0wgaJpZM4VkJ4Q> .

[billsacks]

@slevisconsulting based on @dlawrenncar 's comment, maybe we should at least consider what it would take to combine this with transient. I don't think it will be terrible, and if it's wanted in the near future anyway, it might be easiest just to go ahead and do it now. Let me know if you want to look at this together or if you want to just take an initial look yourself.

[slevis-lmwg]

@billsacks thank you for validating my intuition about where to start :-) I like your recommendation regarding fert_cft and other variables and will try to tackle the way that you suggested.

I'm happy to consider combining with transient. I do not see that as changing my starting point, ie adding call collapse_all_pfts in surfrdMod.F90. Once the non-transient changes work, I can also add call collapse_all_pfts in subroutine dynpft_interp for transient cases. Does that sound right?

[slevis-lmwg]

@billsacks I've started writing the code for this. In subroutine collapse_all_pfts I need to sort pfts and cfts to determine the dominant ones. This raises the question whether we have a function in the model for sorting. I did not find one when I searched in my branch. Do you know of one? If not, no problem.

[ekluzek]

We do have access to LAPACK in CTSM and there's a subroutine called lasrt in there. There might be other home grown sort routines as well. Since, your lists are going to be short, it doesn't matter too much what algorithm is used. But, good to use something that's well done and already vetted...

[billsacks]

I'm happy to consider combining with transient. I do not see that as changing my starting point, ie adding call collapse_all_pfts in surfrdMod.F90. Once the non-transient changes work, I can also add call collapse_all_pfts in subroutine dynpft_interp for transient cases. Does that sound right?

I agree with your starting point. However, I don't think transient will be quite that simple, because right now dynpft_interp is called before dyncrop_interp. It may be as easy as swapping the order of those two in dynSubgrid_driver. But it's possible that you'll need to rework these so that they first read all of the data for the new year (both pft and cft), then do the collapsing, then do the final bits of code in those routines. I'm not sure....

[billsacks]

Regarding sorting, @ekluzek 's suggestion of using lasrt seems like a good one. I hadn't realized that existed, but seeing http://www.netlib.org/lapack/explore-3.1.1-html/dlasrt.f.html#DLASRT.1 this seems like as good a choice as any.

[billsacks]

Actually, regarding sort: I realized that most out-of-the-box sorting routines (including this lapack one) won't do exactly what you want: Most sorting routines just sort the data, but don't tell you the resulting dominant indices into the data. So I think it's easier to roll our own here rather than trying to dig up a routine that does exactly what we want. I just whipped something up: see billsacks@15bf89c . If this will work for you, you can merge my branch (find_k_max_indices) into yours.

[billsacks]

I added a commit to my branch and turned it into a PR for easier reviewing: see #583

[slevis-lmwg]

@billsacks two questions about merging your branch into mine:

1. I have uncommitted mods in my branch. Do I need to commit them before merging your branch into mine? Not that this would be a problem; just that I would have likely kept developing before committing.
2. The merge itself should look like this?
   git merge 15bf89c

[billsacks]

I have uncommitted mods in my branch. Do I need to commit them before merging your branch into mine? Not that this would be a problem; just that I would have likely kept developing before committing.

I think it's possible to do a merge when you have uncommitted changes, but it's not recommended. From git help merge:

Warning: Running git merge with non-trivial uncommitted changes is discouraged: while possible, it may leave you in a state that is hard to back out of in the case of a conflict.

The general recommendation is to either (1) commit or (2) stash your changes before merging.

Regarding option (2): have you used git stash before? I'll admit that I often get tripped up when trying to use it, but it can be helpful in this case. The main way I get tripped up is when I have a mix of staged and unstaged changes. If all you have are unstaged changes (i.e., you haven't done a git add yet), then stashing can work pretty smoothly. The workflow is:

git stash
git merge ...
git stash pop

git stash saves your uncommitted changes in a stash stack, then backs them out. git stash pop applies the most recently stashed changes, then removes them from the stack.

If you don't want to deal with git stash, then I would manually save a git diff (e.g., git diff > my_diffs) before doing the merge, to help you recover in case of problems.

• The merge itself should look like this?
  git merge 15bf89c

First, if you haven't added my fork as a remote, you'll need

git remote add billsacks https://github.com/billsacks/ctsm.git
git fetch billsacks

Then, note that I have pushed a second commit. The easiest thing is:

git merge billsacks/find_k_max_indices

[slevis-lmwg]

Ok, merge complete. (I was only conceptually aware of git stash before.) Thanks!

[slevis-lmwg]

Nothing committed, yet.
Made first attempt to build; corrected errors.
Second attempt to build thwarted by unscheduled cheyenne outage.

Here's a quick summary of work so far:

1. clm_varpar.F90
   + ! TODO User will set n_dom_soil_patches in namelist
+ ! Namelist scripts or model code will:
+ ! Default n_dom_soil_patches to maxsoil_patches determined from input data
+ integer, public :: n_dom_soil_patches = 2 ! # of dominant soil patches; determines the number of active soil patches

2. controlMod.F90
   Introduced n_dom_soil_patches using maxpatch_pft or similar as template

3. surfrdMod.F90
   After call collapse_crop_types, I have added
   call collapse_all_pfts(...)
   call collapse_crop_var(...)

4. surfrdUtilsMod.F90
   Introduced the new subroutines called from surfrd_veg_all in surfrdMod.F90

Subr. collapse_all_pfts:
a) Normalizes natpft and cft weights by their landunit weights in the gridcell and concatenates into a single array wt_all_patch(natpft_lb:cft_ub) before calling find_k_max_indices.
b) Adjusts wt_nat_patch, wt_cft, and their respective wt_lunit by normalizing the dominant weights to 1 (else they would sum to <=1) and setting the remaining weights to 0. If natpft or crop landunit is now non-existent, merge the two landunit weights.
c) Sets non-dominant pft and cft weights to 0 and is_pft_known_to_model to .false.

Subr. collapse_crop_var sets the passed crop variable to 0 where is_pft_known_to_model has been changed to .false.

[slevis-lmwg]

Updates and discussion will continue in #588.

[slevis-lmwg]

[...] I don't think transient will be quite that simple, because right now dynpft_interp is called before dyncrop_interp. It may be as easy as swapping the order of those two in dynSubgrid_driver. But it's possible that you'll need to rework these so that they first read all of the data for the new year (both pft and cft), then do the collapsing, then do the final bits of code in those routines. I'm not sure....

For the record: I tested gnu standard testing with the calling order reversed and it PASSed.

However, you are correct @billsacks that I need to rework these. So in the end reordering the calls is not relevant.