Change renormalization of landunit areas in mksurfdata_map to be consistent with raw datasets

[billsacks]

In recent discussions (spinning off from #1445 , including @Face2sea @fang-bowen @keerzhang1 @olyson and @lawrencepj1 ), it has become clear that there are problems with the current renormalization of PCT_NATVEG and PCT_CROP in mksurfdata_map to ensure that all landunits add up to 100%. This refers mainly to the code in normalizencheck_landuse and this initial adjustment:

CTSM/tools/mksurfdata_map/src/mksurfdat.F90

Lines 751 to 763 in 7a85024

	suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
	if (suma > 250._r4) then
	write (6,*) subname, ' error: sum of pctlak, pctwet,', &
	'pcturb and pctgla is greater than 250%'
	write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', &
	n,pctlak(n),pctwet(n),pcturb(n),pctgla(n)
	call abort()
	else if (suma > 100._r4) then
	pctlak(n) = pctlak(n) * 100._r8/suma
	pctwet(n) = pctwet(n) * 100._r8/suma
	pcturb(n) = pcturb(n) * 100._r8/suma
	pctgla(n) = pctgla(n) * 100._r8/suma
	end if

My sense is that these adjustments made sense at one point (e.g., maybe for CLM4.0), but have not been updated to remain correct with other recent changes – one big one being the standard placement of crop on its own landunit.

Specific problems with what's currently done are:

• Urban should not preferentially replace bare ground: all PFTs should keep their original proportions. (The replacement of bare ground with urban may have made sense originally (see Implementing dynamic urban land cover in CESM #1445 (comment) ), but @lawrencepj1 says that it no longer makes sense with recent raw datasets, where urban is split out from bare ground.) Urban should be treated the same as other special landunits now.
• Both crop and natural veg are decreased to accommodate too-large special landunit areas. This is a problem because, in the creation of the raw datasets, only natural veg is increased beyond what's specified initially (to accomplish PCT_CROP + PCT_NATVEG = 100%), so it is inconsistent that we then decrease crop to accommodate extra special landunit areas

My proposal, based on discussions with @lawrencepj1 and @olyson, is:

• PCT_NATPFT and PCT_CFT should never be adjusted beyond what's on the raw datasets. In particular, this means that urban should not adjust PCT_NATPFT.
• The areas of all landunits other than natural veg should be treated symmetrically; this should include the crop landunit, which should now be treated similarly to special landunits in places where we assure sums don't exceed 100% and in places where we adjust all other landunits proportionally. In particular, do not decrease crop to accommodate too-large special landunit areas. The justification is that it feels equally important to capture the correct areas of all of these non-natveg landunits (including crop!).
• PCT_NATVEG should be the residual after all other landunits. So, effectively, too-large areas of non-natveg landunits will first decrease the area of natural veg; only after that will the non-natveg landunits themselves decrease below what is specified. It's arguable if this makes the most sense for a present-day surface dataset (it would arguably make more sense to decrease /everything/ proportionally – though if the input datasets were consistent, then this would be a non-issue, so it's hard to say what is really right), but this feels most correct for transient cases, particularly when the natural veg area is an unknown that is just the residual of other areas.
  • Note, though, that we'll still need to read PCT_NATVEG initially for at least one purpose: giving the weighting of each source point in the remapping of PCT_NATPFT: e.g., a source cell that is 100% natural veg should have twice the weight of a source cell thath is 50% natural veg in remapping the PFT fractions from source to model resolution. However, for determining landunit areas, PCT_NATVEG will be set to 100 minus the sum of the other landunits.
• For the reason given above (the need for PCT_NATVEG in order to provide source weightings), the raw dataset creation should be changed to not rescale PCT_NATVEG such that it plus PCT_CROP equal 100%. I'll open a separate issue for this.

To accomplish this:

• Look at the sum of all landunits other than natural veg (but including crop!): If this sum is > 100%, scale down all landunits other than natural veg (but including crop!) proportionally. This will replace the suma rescaling code that is done just before calling normalizencheck_landuse.
• Set PCT_NATVEG to 100 minus the sum of the other landunits. This will replace the current rescaling code in normalizencheck_landuse.

Note that this scaling down process differs from what's done in dynLandunitAreaMod: in the latter, landunits are decreased in a priority order. There is some argument for keeping the two consistent, although there is also some argument for allowing them to remain inconsistent. The argument for the latter is: The corrections in mksurfdata_map relate to reconciling raw datasets that disagree but should (in principle) be consistent; the corrections in dynLandunitAreaMod are more about handling dynamic areas from other components (e.g., dynamic glacier). It may be that the best solution is to just save all of the normalization to 100% to runtime (see my comment near the end of #1445 (comment) ), but it feels like that introduces a lot of potential for new issues, for uncertain benefit, so we probably won't go there for now.

When making these changes, we should double-check that there isn't any code that assumes that PCT_NATVEG + PCT_CROP adds up to 100% (or forces that to be the case). I didn't see anything like this in my initial, quick read-through, but we should double-check that. One specific thing to look for is: will code work correctly if the input PCT_NATVEG = PCT_CROP = 0%? (Or would such grid cells incorrectly be considered ocean, for example?)

• Edit (2021-12-20): I think this last (bold) paragraph is mainly or entirely relevant in the context of Remake PCT_PFT raw datasets to NOT force landunit area to 100%, and some other fixes needed to PCT_PFT raw datasets #1556, not for this issue, so I have copied it to there.

[lawrencepj1]

Thanks @billsacks and @olyson for this post and the meeting. Yes this is the approach I would support as well except I don't think there is a problem leaving the raw as it is. Having the PCT_NATVEG being the residual of the other land units ensures that the PCT_NAT_PFT x PCT_NATVEG (generated in the model at runtime) should be correct regardless of the other land units.

[olyson]

For the first task (Look at the sum of all landunits other than natural veg (but including crop!): If this sum is > 100%, scale down all landunits other than natural veg (but including crop!) proportionally.), does this bit of code accomplish what we want? Here I've included the crop landunit in the rescaling.

      ! Make sure sum of all land cover types except natural vegetation
      ! does not exceed 100. If it does, subtract excess from these land cover
      ! types proportionally.

      do n = 1,ns_o
         suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
         if (suma > 250._r4) then
            write (6,*) subname, ' error: sum of pctlak, pctwet,', &
                 'pcturb, pctgla, and pctcrop is greater than 250%'
            write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctcrop= ', &
                 n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctcft(n)%get_pct_l2g()
            call abort()
         else if (suma > 100._r4) then
            pctlak(n) = pctlak(n) * 100._r8/suma
            pctwet(n) = pctwet(n) * 100._r8/suma
            pcturb(n) = pcturb(n) * 100._r8/suma
            pctgla(n) = pctgla(n) * 100._r8/suma
            call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma)
         end if
      end do

[billsacks]

Yes, that looks right - thanks!

[billsacks]

I don't think there is a problem leaving the raw as it is. Having the PCT_NATVEG being the residual of the other land units ensures that the PCT_NAT_PFT x PCT_NATVEG (generated in the model at runtime) should be correct regardless of the other land units.

@lawrencepj1 sorry - I didn't see this comment until just now. (I think you edited your original comment, but edits don't generate a notification, so I only noticed this addition when I was rereading through this whole issue.)

To first order, you're right. However, I think the forcing of PCT_NATVEG to 100% on the raw datasets causes problems with the remapping from source to destination. I explain this more clearly in #1556 - particularly in the example I give there. Can you look through that and see if it makes sense to you? Let's move further discussion of this point to #1556 .

[billsacks]

As discussed in the ctsm-software meeting today and in comments in #1587 , this is now on the 5.2 mksurfdata branch, so I'm closing this issue.

fyi @olyson @slevisconsulting