Change mapconsf mappings to mapconsd?

[billsacks]

Emerging from a discussion of ESCOMP/mizuRoute#426, I talked with @mvertens a bit about the use of mapconsf vs. mapconsd. I found an email from her from Nov, 2020 - where the context was some changes she was making to the lnd <-> glc mapping - where she said:

Actually, everything should be mapconsd (long story - I spent a week with Bob and Denise trying to sort this out). For CESM it is only round of level changes at most - but I am trying to add this back in incrementally as we move forwards and test.

However, it looks like many mappings still use mapconsf. I believe this is part of the problem behind ESCOMP/mizuRoute#426: the main lnd -> rof mapping is treated as mapconsf, but the lnd -> rof mapping of lfrac uses mapconsd, and this discrepancy means that the user-provided lnd -> rof mapping file is ignored for the mapping of lfrac.

This morning, @mvertens told me that, although she, @oehmke and @DeniseWorthen came to agree that we should be using mapconsd (at least for global grids), at the time they wanted to avoid answer changes (and/or didn't have time to do the full validation this would require), so they did not make this change. @mvertens and I felt like it would be worth revisiting this to answer the following questions:

1. Do we still agree that mapconsd should be used for all mappings? (I would like to understand this myself.)
2. Does this hold for regional grids as well as global grids?
3. For global grids, is it true that we just expect roundoff-level changes from this? Is that true for both CESM and UFS?
4. Is soon (say, within the next few months) a good time to make this change?

Things are pretty busy now both for @mvertens and the ESMF team, but maybe we could meet in a few weeks to discuss these things – or just have some discussion here.

[ekluzek]

@billsacks thanks for starting the discussion. I was going to create an issue for the mizuRoute problem, but this should cover it for us if we make this change. Those are the questions I want to know about as well, and right now I certainly do not.

Looking at the code the difference is in how the normalization is done. Here is med_map_mod.F90:

    else if (mapindex == mapconsf .or. mapindex == mapnstod_consf) then
       if (mastertask) then
          write(logunit,'(A)') trim(subname)//' creating RH '//trim(mapname)//' for '//trim(string)
       end if
       call ESMF_FieldRegridStore(fldsrc, flddst, routehandle=routehandles(mapconsf), &
            srcMaskValues=(/srcMaskValue/), &
            dstMaskValues=(/dstMaskValue/), &
            regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
            normType=ESMF_NORMTYPE_FRACAREA, &
            srcTermProcessing=srcTermProcessing_Value, &
            ignoreDegenerate=.true., &
            dstStatusField=dststatusfield, &
            unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
            rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       ldstprint = .true.
    else if (mapindex == mapconsf_aofrac) then
.
.
.
    else if (mapindex == mapconsd .or. mapindex == mapnstod_consd) then
       if (mastertask) then
          write(logunit,'(A)') trim(subname)//' creating RH '//trim(mapname)//' for '//trim(string)
       end if
       call ESMF_FieldRegridStore(fldsrc, flddst, routehandle=routehandles(mapconsd), &
            srcMaskValues=(/srcMaskValue/), &
            dstMaskValues=(/dstMaskValue/), &
            regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
            normType=ESMF_NORMTYPE_DSTAREA, &
            srcTermProcessing=srcTermProcessing_Value, &
            ignoreDegenerate=.true., &
            dstStatusField=dststatusfield, &
            unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
            rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       ldstprint = .true.

The difference is moving from normType=ESMF_NORMTYPE_FRACAREA to ESMF_NORMTYPE_DSTAREA. I would think this change would be greater than roundoff.

I think I'll need a temporary change for mizuRoute sooner than this change will happen so I'll make a separate issue for that.

[billsacks]

@mvertens just clarified: "in terms of roundoff - I was only using that term lightly. What I really meant is that answers should just have very little differences."

[billsacks]

I have spent some time thinking about this issue this week and talking with @mvertens and @oehmke about it... this is in the context of #665 but we did some more general thinking to get our heads back into this issue.

We didn't reach a point of complete confidence, and we felt this is worth one more mental revisit before actually making this change. But we have come to understand why we think consd and consf mapping should generally give the same answers in CMEPS (within roundoff) given how things are typically set up there. The key point is that masks come into play not only in the consf mapping but also in the normalization in med_map_field_normalized. As long as a given mapped field has a specified normalization field (mapnorm in the addmap_from call), and as long as that normalization field is 0 wherever the mesh's mask is 0, I think consd and consf mappings end up being equivalent. (See details at the bottom of this issue for a worked example.)

I think these conditions hold for the lnd2rof fields, where their mapnorm is set to landfrac in lines like this:

CMEPS/mediator/esmFldsExchange_cesm_mod.F90

Line 3323 in 806cce1

call addmap_from(complnd, 'Flrl_rofsur'//trim(suffix), comprof, mapconsf, map_fracname_lnd2rof, 'unset')

For a component where masks and fractions are always 1, such as atm, consd and consf mappings are the same. We didn't give a lot of thought to regional grids, but my understanding is that things should work right for regional grids, too:

• For a component that already has masks / fractions, such as lnd, I think a regional grid is equivalent to a global grid with masking as far as this mapping is concerned - i.e., I think that the existing handling of masks / fractions in the normalization should handle regional grids in the same way that it currently handles masked global grids.
• For a component that does not have masks / fractions, such as atm, I think this is still handled correctly as long as mapnorm is set to 'one'. I have spot-checked some atm fields in esmFldsExchange_cesm and they do have mapnorm set to 'one'. (I haven't extensively checked everything.) My understanding is that this maps a field of all 1s from the source to destination and uses this as the normalization field. I think this is like saying that the frac is 1 wherever the grid exists, but is 0 outside the defined grid.

If there are any fields that don't have mapnorm set in one of these ways, then consd and consf mappings could differ... if there are any cases like that, we'd need to think more carefully about whether we want consf or consd, or whether these fields should have mapnorm set.

In the discussion with @mvertens and @oehmke , we weren't sure why the earlier discussion had led to the conclusion that we should be using consd everywhere. It seems that either consd or consf would work fine.

• Arguments I see for consd:
  • I think we need consd mapping for mapping fractions (in med_fraction_mod) (I think that using consf mapping gives wrong results when mapping fractions) - so if part of the goal is to use a single mapping method everywhere, this may be required
  • It's simpler to understand what's happening
  • It doesn't rely on the mask being defined on the mesh at the point of doing the FieldRegridStore
  • It might have a minor performance advantage
• Arguments I see for consf:
  • It's what's currently used in most places, so there is less need for validation and less chance to introduce subtle bugs in changing things
  • It provides some redundancy in case mapnorm isn't set right... though I'm actually not sure if this is a pro or a con: it may be better to force mapnorm to be set appropriately rather than doing something that makes answers generally look right when in fact the implementation is subtly wrong

We should check the following in the code when making these changes:

• Look at each addmap call to verify that mapnorm is set correctly, since this would become more critical if we move to consd.
• Check custom mappings (e.g., lnd2glc, some lnd2rof, etc.) to make sure there aren't assumptions there that require either consd or consf mapping.

Validation is likely to be painful here. I still feel this is worth doing because - as seen with #665 and #409 - it is causing significant technical debt to have these different mapping types in different places. But this need for validation means that this will probably require dedicating a couple of weeks to this task. We should test any changes to ensure that they only lead to roundoff-level differences. In a configuration with active components, there is likely to be a propagation of differences, so we'd probably want to check things after the first coupling interval. Alternatively, we could use an X compset, as long as it exercises any custom mappings (which could be some of the most critical things to check).

Scenario demonstrating equivalence of consd and consf mapping given normalization

Scenario

Single destination (ROF) grid cell with 4 source (LND) grid cells, each with 25% overlap:

(1) S1: Mask=1, Frac=1, Value=10

(2) S2: Mask=1, Frac=0.5, Value=12

(3) S3: Mask=0, Frac=0

(4) S4: Mask=0, Frac=0

Note that lfrac = 0.375.

With dstarea / consd mapping

W1 = 0.25, W2 = 0.25.

frac_dest = 0.25*1 + 0.25*0.5 = 0.375

val_dest = (0.25*10 + 0.25*12*0.5)/frac_dest = 10.67

After the merge, this becomes 10.67*lfrac = 4.0

With fracarea / consf mapping

W1 = 0.5, W2 = 0.5.

frac_dest = 0.5*1 + 0.5*0.5 = 0.75

val_dest = (0.5*10 + 0.5*12*0.5)/frac_dest = 10.67

After the merge, this becomes 10.67*lfrac = 4.0

The key observation

Initially, I thought that fracarea mapping together with the normalization would be double-accounting for the source fraction. But the key point is that the impact of fracarea mapping cancels out with itself, since it appears both in the mapping of fractions (to create frac_dest) and in the mapping of the values (to create val_dest). So, given this normalization, consd and consf mapping give equivalent results.

[ekluzek]

@billsacks I'm just wondering about what you think the timing of this should be? I gather this is a good thing to do longer term, but post-release. Doing it pre-release would require a bunch of time on validation that probably no one could dedicate the time too (either in CESM or NorESM because of CMIP7 priorities). And it would be another potential answer change that potentially causes problems for CMIP7 or release priorities. So I'm assuming this wouldn't really be worked on until late Fall this year or into 2027, but let me know if you're thinking differently than that.