More rework needed of mksurfdata_map code to set wetland percent

[billsacks]

Brief summary of bug

The code to set wetland area was reworked on the release-clm5.0 branch and then merged to master in ctsm1.0.dev093. However, I think there are outstanding issues with the current code.

General bug information

CTSM version you are using: ctsm1.0.dev093

Does this bug cause significantly incorrect results in the model's science? Yes, but probably fairly minimal

Configurations affected: Potentially all, but probably fairly small impacts in general

Details of bug

Relevant issues are #553 and #673 (also #545 for background/context). It looks to me like there are the following outstanding issues:

1. The current code doesn't implement the desired logic as stated in Ice shelf wetland fix in mksurfdata_map can lead to glacier+lake > 100% on surface datasets #673 (comment) and Ice shelf wetland fix in mksurfdata_map can lead to glacier+lake > 100% on surface datasets #673 (comment)

2. There are some possible issues with error-check code that I think need to be addressed: see More robust mksurfdata_map logic for determining where to put wetlands #553 (comment)

There may be other outstanding issues, too - #553 and #673 deserve a careful re-read before marking this as complete - but those are the two that jump out at me. Note that I think that my proposal in #553 (comment) would effectively be addressed by the logic in #673 because urban is set to 0% in the relevant block of code... but this again deserves more careful thought than what I have given it just now.

Update (2020-08-21)

See #1001 (comment) for my updated thoughts on this, which differ from what is laid out above.

[ekluzek]

I think we need to have a discussion about this issue. Should we wait for our Thursday meeting, or schedule something sooner?

[billsacks]

Thursday meeting seems fine.

[billsacks]

I started looking back at this. I can't understand exactly what I was originally thinking with my above-referenced comments in #673 , but my best guess is that I made a mistake in my thinking there. Specifically: it looks like I may have been assuming that, if the relevant code block ended up with pctlak + pctgla > 0 but < 100, and no areas of other landunits, then pctlak and pctgla would be renormalized so that those two landunits sum to 100%. But it doesn't look like that's how things work right now. It may be that the right solution is to go with my comments in #673 but then also add some code (probably in normalizencheck_landuse) so that, in areas outside the pft landmask, we first normalize all special landunits so that they sum to 100%.

I want to look back at this with fresh eyes next week.

[billsacks]

Okay, I think I have gotten my head somewhat back into this. I think my motivation here was seeing that we're setting the wetland % using fundamentally different logic for gridcells inside vs. outside the pft mask:

• Inside the pft mask, we are keeping wetlands at 0%, even if a grid cell is only partially inside this mask. Imagine a grid cell that is 50% inside and 50% outside the pft mask and has 0% special landunits (so presumably it is 1/2 land and 1/2 ocean). I think we end up making this grid cell 100% vegetated, using the vegetation distribution from the portion of the grid cell that is within the pft mask.

• Outside the pft mask, however, wetlands can be set to a fractional value. For example, imagine a grid cell that is outside the pft mask, but the glacier dataset says it's 50% glacier (so as above, it's 1/2 land – thinking of glacier as land – and 1/2 ocean). Now we end up making this grid cell 50% glacier and 50% wetland. So there's an inconsistency with what we do in the first case.

I think this inconsistency was my motivation for suggesting that we change the logic for setting pctwet so that we keep it at 0% if either lake or glacier have > 0% cover. However, I think this would involve a bit more rework of the code to ensure that lake + glacier + wetland end up summing to 100% in this case. I don't think this would be too hard, but it feels like changes like this in mksurfdata_map often ends up being harder than I initially think because of edge cases. So I'm inclined to skip this change in behavior.

There is one other issue referenced above that I think we should fix: point (1) in #553 (comment) (point (2) is only relevant if we make the other changes discussed above, I think) (edit (2020-08-21): actually, point (2) may be totally irrelevant, arising from faulty logic on my part, because with the other suggestions here, I think we shouldn't have any vegetated area outside the pft data mask... it may be that, at the time I made that comment, I was thinking that we'd fill the rest of the cell with something like bare ground, but now I'm not sure that makes sense: e.g., we don't want bare ground over the Antarctic ice shelves). I think this may be important to avoid errors or at least very wrong behavior in some edge cases where we have < 0.5% cover of glacier and/or lake. I think the best thing to do would be to move the whole block of truncation code to above the wetland adjustment. i.e., move this block:

CTSM/tools/mksurfdata_map/src/mksurfdat.F90

Lines 737 to 747 in 6c99904

	! Truncate all percentage fields on output grid. This is needed to
	! insure that wt is zero (not a very small number such as
	! 1e-16) where it really should be zero
	do k = 1,nlevsoi
	pctsand(n,k) = float(nint(pctsand(n,k)))
	pctclay(n,k) = float(nint(pctclay(n,k)))
	end do
	pctlak(n) = float(nint(pctlak(n)))
	pctwet(n) = float(nint(pctwet(n)))
	pctgla(n) = float(nint(pctgla(n)))

to above this block:

CTSM/tools/mksurfdata_map/src/mksurfdat.F90

Lines 714 to 735 in 6c99904

	! Assume wetland, glacier and/or lake when dataset landmask implies ocean
	! (assume medium soil color (15) and loamy texture).
	! Also set pftdata_mask here
	if (pctlnd_pft(n) < 1.e-6_r8) then
	pftdata_mask(n) = 0
	soicol(n) = 15
	if (pctgla(n) < 1.e-6_r8) then
	pctwet(n) = 100._r8 - pctlak(n)
	pctgla(n) = 0._r8
	else
	pctwet(n) = max(100._r8 - pctgla(n) - pctlak(n), 0.0_r8)
	end if
	pcturb(n) = 0._r8
	call pctnatpft(n)%set_pct_l2g(0._r8)
	call pctcft(n)%set_pct_l2g(0._r8)
	pctsand(n,:) = 43._r8
	pctclay(n,:) = 18._r8
	organic(n,:) = 0._r8
	else
	pftdata_mask(n) = 1
	end if

This looks like a safe thing to do, but I'd love an extra set of eyes on it.

[wwieder]

I agree this last suggestion seems safe, but also have 0% experience with mksurfdata_map.

[billsacks]

Revisiting this today, I stand by my thinking that we should NOT address the main changes referenced here. In fact, if anything, I think that the faulty logic is in what we do with wetland % in grid cells with pftdata_mask > 1e-6: I'm starting to think that our handling of wetlands with respect to lake and glacier is correct (in places with pftdata_mask = 0, fill the remainder of the grid cell with wetlands where it wasn't filled with glacier or lake), and that the real problem is that, anywhere where we have pftdata_mask > 1e-6, we fill the whole grid cell with vegetation.

I'm picturing shifting a grid cell in space: Let's ignore glaciers and lakes for now. Imagine a grid cell near the coast that is 100% within pftdata_mask, so it is 100% natural veg. Now imagine gradually moving that grid cell, so that it starts to encompass a greater and greater fraction of ocean. I believe that, with the current logic, the grid cell will stay 100% natural veg until it is (almost) entirely over the ocean, at which point it will snap to 100% wetland. It seems like it would be better if this transition were gradual.

So again, if we were to change anything here, I think we should change the behavior of wetlands with respect to vegetation, not the behavior of wetlands with respect to lakes and glaciers as I suggested elsewhere in this issue. So, something like this (note: this code assumes that, after mkpft, pctnatveg + pctcrop = 100 everywhere, as documented at the top of mkpft, so we just need to consider pctlnd_pft in the following, not pctnatveg or pctcrop):

pctwet(n) = max(100._r8 - pctlnd_pft(n) - pctgla(n) - pctlak(n), 0.0_r8)

But I don't think this needs to be changed now.

I do still think we want to move the truncation code to a bit earlier, as mentioned above. In addition to the edge case where pctlak or pctgla is < 0.5%, this could also be relevant in a case like this: pctlak = 5.4%, pctgla = 94.4%. With the current code, we'd first get pctwet = 0.2%. But then, after the rounding, we'd end up with pctlak = 5%, pctgla = 94%, pctwet = 0%. Then later, in normalizencheck_landuse, I think we'd end up with 1% natural veg, which could trigger the 'vegetation found outside the pft mask' abort. (Caveat: I may be missing something here, since there are so many corrections done at different points. It's also possible that there are other problems with doing the rounding first, but I'm not seeing any off-hand.)

To keep things clean, I'm going to go ahead and close this issue as a wontfix, and open a new issue for the one piece I think we should change.