mksurfdata: Rework mapping of landunits to handle coastal areas more rigorously

[billsacks]

Description of changes

This PR reworks the mapping of landunit areas and the determination of wetland areas according to the plan laid out in #1716. (Note that the final implementation differs from the initial comment in #1716, as discussed in follow-up comments in that issue.) See the discussion there as well as the design document in this PR for details.

Specific notes

Contributors other than yourself, if any: Various discussions with others, especially @ekluzek and @slevisconsulting

CTSM Issues Fixed (include github issue #):

• Resolves Wetland area and areas of landunits in coastal grid cells should be determined more rigorously #1716
• Resolves North and South pole set to glacier on surface datasets regardless of raw data specification #1930

Are answers expected to change (and if so in what way)? Areas will change on surface datasets

Any User Interface Changes (namelist or namelist defaults changes)? No

Testing performed, if any:

(1) Compared 1850 f09 surface dataset before and after these changes (and after each individual change from the individual commits on this branch)

(2) By hard-coding landunit areas before the call to normalize_and_check_landuse, I performed various "manual unit tests" of the normalization code, using diffs like this:

diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90
index eb91be30c..0f347f8c0 100644
--- a/tools/mksurfdata_esmf/src/mksurfdata.F90
+++ b/tools/mksurfdata_esmf/src/mksurfdata.F90
@@ -593,6 +593,17 @@ program mksurfdata
   ! Perform other normalizations
   ! -----------------------------------
 
+  do n = 1, lsize_o
+     pctlnd_pft(n) = 0._r8
+
+     call pctnatpft(n)%set_pct_l2g(0._r8)
+     call pctcft(n)%set_pct_l2g(0._r8)
+     pctlak(n) = 0._r8
+     pcturb(n) = 0._r8
+     pctgla(n) = 0._r8
+     pctwet(n) = 0._r8
+  end do
+
   ! Normalize land use and make sure things add up to 100% as well as
   ! checking that things are as they should be.
   call normalize_and_check_landuse(lsize_o)

In the following, any landunit not explicitly mentioned had 0% area:

(a) With pctlnd_pft = 0% and all initial landunit areas 0%: wetland area = 100%

(b) With pctlnd_pft = 0% and initial glacier area 1%: glacier area = 100%

(c) With pctlnd_pft > 0% and all initial landunit areas 0%: natural vegetated area = 100%

(d) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%: crop area = 50%, natural vegetated area = 50%

(e) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 10%: crop area = 50%, natural vegetated area = 25%, glacier area = 25%

(f) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 15%: crop area = 50%, natural vegetated area = 12.5%, glacier area = 37.5%

(g) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 20%: crop area = 50%, glacier area = 50%

(h) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 30%: crop area = 40%, glacier area = 60%

(i) With pctlnd_pft = 40%, initial crop area 0%, natural vegetated area 40%, glacier area 40%: glacier area = 100%

(j) With pctlnd_pft = 2%, initial natural vegetated area 1%, glacier area 1%: natural vegetated area = 50%, glacier area = 50%

(k) With pctlnd_pft = 2%, initial natural vegetated area 0%, glacier area 1%: natural vegetated area = 50%, glacier area = 50%

(l) With pctlnd_pft = 2%, initial natural vegetated area 2%, glacier area 1%: natural vegetated area = 50%, glacier area = 50%

[billsacks]

@ekluzek @slevisconsulting - I just pushed a somewhat-related change to this branch, where I have removed the code that hard-codes the pole points to glacier. I never liked that we had this code - e.g., this could trip someone up if they were trying to create a surface dataset for a very different climate / continental configuration - and it was leading to weirdness in the landfrac at the north pole (setting that landfrac to 1, which wasn't a problem per se, but was weird). I looked at diffs in an f09 surface dataset from removing it, and the only diff is as expected: the north pole is now classified as wetland rather than glacier; there is no change at the south pole (which is still classified as glacier).

Please let me know if you can think of any reason why this hard-coding at the poles should be retained.

[billsacks]

I have looked at the difference in LANDFRAC_MKSURFDATA for an f09 transient dataset, between 1850 and 2005. There are only 8 grid cells that have any difference in LANDFRAC_MKSURFDATA between years 1850 and 2005. Of these, only 1 has a difference exceeding 1e-4: one grid cell has a difference of 1.45e-3; this is a grid cell at the mouth of the Mediterranean Sea. This is good: there are NOT big differences in our estimated land fraction in a transient case.

@ekluzek and @slevisconsulting - at this point I feel like this is ready to be merged if you're satisfied with it.

[slevis-lmwg]

@ekluzek and @slevisconsulting - at this point I feel like this is ready to be merged if you're satisfied with it.

@billsacks I'm looking through the code and will respond soon.

[slevis-lmwg]

@billsacks thank you for your careful work in this PR! Looks good to me. I'm clicking "approve" and will let you and @ekluzek decide whether you also need his approval.

[ekluzek]

@billsacks the change you made on the poles in general makes sense to me. Although it does seem a bit weird if the north pole is now wetland rather than glacier. I'm just wondering if you found some information in git logs (or elsewhere) about why this was done way back in the day? That context can be really helpful in figuring out if we should retain a feature we now question.

[billsacks]

@ekluzek I just checked the history and this adjustment of the pole points dates all the way back to the initial version in the CLM svn repository. If I had to guess, my guess would be that it was needed at one time due to problems with mapping to the poles.

Why do you think it's weird if the north pole is now wetland? That is consistent with the raw datasets, where the north pole looks like ocean and isn't claimed by any landunit. This makes the north pole end up with the same land cover as the latitudes slightly south of the north pole, which seems more correct and consistent to me.

[slevis-lmwg]

[...] That is consistent with the raw datasets, where the north pole looks like ocean and isn't claimed by any landunit. This makes the north pole end up with the same land cover as the latitudes slightly south of the north pole, which seems more correct and consistent to me.

Am I right to also conclude that, if a user put a continent over the North Pole with evergreen trees, then the NP would be land with evergreen trees? (I.e. neither glacier nor wetland.)

[billsacks]

Am I right to also conclude that, if a user put a continent over the North Pole with evergreen trees, then the NP would be land with evergreen trees? (I.e. neither glacier nor wetland.)

With the new code, yes: the north pole should end up with whatever is specified on the raw datasets. With the old code, the north pole was forced to glacier regardless of what was specified on the raw datasets.

[ekluzek]

Ahh, OK, thanks for that clarification. I wasn't clear what this really looked like. It sounds like you used to have a completely isolated glacier point around the North pole that was disconnected to anything, the new method does make more sense. We did used to have more issues with remapping, so it makes sense it had something to do with that.

[billsacks]

It sounds like you used to have a completely isolated glacier point around the North pole that was disconnected to anything

Yes, exactly. Here's what PCT_GLACIER looked like before my recent change: notice the red band (100% glacier) in the far northern latitude band:

(In practice this probably didn't matter, because 90° North would typically be ocean, but if it were ever land for some reason, it would be classified differently from the nearby grid cells.)

[slevis-lmwg]

I talked with @ekluzek this morning and we agreed:

• He will look over this PR before it gets merged.
• Once ready, do you @billsacks have a preference for who does the merging/tagging? I think any of the three of us could do it.

[billsacks]

• Once ready, do you have a preference for who does the merging/tagging? I think any of the three of us could do it.

Thanks. No preference.

[ekluzek]

OK, I have a few questions and comments. And mostly I wanted to make sure I understood this a little better by looking through it.

The general idea is that the routehandles now normally are normalized by fractional area, and if they are NOT they have "_nonorm" at the end to say they aren't.

There's also a bug that you caught that made the North pole glacier, which is really good to catch. I'm thinking that would be good to make into an issue for visibility across the project.

[billsacks]

The general idea is that the routehandles now normally are normalized by fractional area, and if they are NOT they have "_nonorm" at the end to say they aren't.

There's also a bug that you caught that made the North pole glacier, which is really good to catch. I'm thinking that would be good to make into an issue for visibility across the project.

Yes, thanks for the summary. I just opened #1930 and marked it as resolved in the initial comment in this PR. In general, raw datasets on fractional area are mapped without normalization, whereas most parameters are mapped with normalization.

[ekluzek]

I'm going to go ahead and merge this in.