Nearest neighbor destination to source interpolation for runoff

[dazlich]

The port of mpas-ocean into the CESM framework is not getting all of the runoff (rofl) from DROF. I believe this is due to the mpas-ocean grid not being global - only defined where there is active ocean. With offline testing I have found that a nearest neighbor destination to source interpolation will conservatively interpolate the runoff to the ocean grid. Currently, I do not see this option in med_internalstate_mod.F90.
Right now, as I see it I need to:
Add a nearest neighbor d to s option in med_internalstate_mod.F90 based on ESMF_REGRIDMETHOD_NEAREST_DTOS (I propose mapndtos).
Add a block to med_map_mod.F90 to call ESMF_FieldRegridStore for this interpolation.
define rof2ocn_liq_rmap somewhere so that it is not 'unset' in esmFldsExchange_cesm_mod.F90, but is set to mapndtos.
I've likely oversimplified and am missing things as well.

[jedwards4b]

@billsacks @mvertens any thoughts on this?

[mvertens]

@dazlich - thanks for looking into this and proposing a fix. Several points:

1. Normally the runoff->ocean mapping uses a smoothed custom mapping file since ocean instabilities due to salinity issues were found when a simple conservative mapping was used. I agree that we should add the new mapping you suggested. However, I am wondering if you don't want to use a smoothed mapping for your case as well. At that point nothing would have to be changed in CMEPS for your case.
2. If you want to rof2ocn_liq_rmap not set to 'unset' in esmFldsExchange_cesm_mod.F90 - we will have to introduce an Earthworks mode for this - or some other type of logical since for CESM cases we always want it to be 'unset'. The value of 'unset' implies that you will be reading in a mapping file. I would suggest that for now you use your wght_file_near.nc rather than modifying CMEPS. To do that out of the box you need to update the ccs_config/maps_noupc.xml to point to your mapping file(s) for the rof->ocn mapping.
3. @jedwards4b, @billsacks - can we try to set up a time to talk about a longer term solution to this to permit different mapping types for rof->ocn?

[dazlich]

@mvertens - thanks for your comments!

1. I agree the new mapping should be smoothed.
2. Ok, I did not understand the use of 'unset' for the mappings. I will try reading my mapping file as you suggest - perhaps that is all I need to do to solve my problem. In that case, I will need to regenerate my wght file so that it smooths as well.

[dazlich]

@mvertens @jedwards4b @billsacks
I followed Mariana's recommendation above in point to to edit ccs_config/maps_noupc.xml to use my mapping file and it works. I guess I was making the problem too complicated. Now I have to incorporate smoothing into the remapping, but that is an offline problem, not a cmeps one.

Thanks for the help!

[billsacks]

@mvertens , @jedwards4b , @gustavo-marques and I met to discuss this on Friday. I should have written up some notes from that meeting at the time, because now that I'm returning to it after the weekend, my memory has gotten very cloudy on what exactly we decided... here is my best recollection but others, please correct me if you remember differently... some of this is probably just restating what @dazlich has already come to realize:

Currently, our standard procedure in CESM is to provide offline-generated runoff mapping files; these do a nearest neighbor mapping but (we think) with an area correction to ensure conservation; then, in some cases, there is an extra horizontal smoothing applied to this nearest neighbor mapping. If you don't provide a runoff mapping file, then the current (untested, in a scientific sense!) behavior is to fall back to a conservative mapping. This is implemented here:

CMEPS/mediator/esmFldsExchange_cesm_mod.F90

Line 2137 in 6d4c1ff

call addmap_from(comprof, 'Forr_rofl', compocn, mapconsd, 'none', 'unset')

and here:

CMEPS/mediator/esmFldsExchange_cesm_mod.F90

Line 2161 in 6d4c1ff

call addmap_from(comprof, 'Forr_rofi', compocn, mapconsd, 'none', 'unset')

However, this probably isn't an appropriate option for this mapping, because – as @dazlich mentioned – there can be runoff in rof gridcells that don't overlap with the ocean grid.

Alternative options that are already implemented would be to use mapnstod (nearest source to destination), mapnstod_consd (nearest source to destination followed by conservative dst) or mapnstod_consf (nearest source to destination followed by conservative frac). (None of these include the smoothing option that we often use in CESM.) At first we thought that mapnstod_consd (or maybe mapnstod_consf) might be an appropriate option, but from looking at how it's implemented, our best sense is that this effectively does a conservative mapping in areas of overlap between the rof and ocn grids and a nearest neighbor (nstod) mapping for ocn cells that are outside the rof grid, and our sense is that this won't actually conserve globally.

I believe that two points would need to be addressed to have global conservation:

(1) As @dazlich pointed out in the initial issue comment, mapnstod doesn't seem like the correct nearest neighbor approach. (We didn't discuss this in our meeting Friday, but I'm coming to see this as I come back to it more carefully today.) Instead, I agree with @dazlich that we would need to use mapndtos (which is not yet implemented in CMEPS) so that each source point goes to exactly one destination point (but a given destination point can receive input from multiple source points).

See http://earthsystemmodeling.org/docs/nightly/develop/ESMF_refdoc/node9.html#opt:regridmethod for details.

(2) From the discussion Friday, we think that, for nearest neighbor to conserve, you would need to do area normalizations before and after the mapping. In effect, you need to transform the per-area flux fields on the runoff grid (specified in mm/s or, equivalently, kg / m^2 / s) into volumetric quantities (specified as kg / s), do the nearest neighbor mapping, then re-transform the mapped fields into per-area fluxes on the ocean grid. I should say: we're pretty sure that this transformation is needed for conservation; what we're less sure of is whether it is already done somewhere (in either CMEPS or ESMF).

I was asked to give some guidance to @gustavo-marques as to what simple changes he could make to test if any of the existing options might conserve, but as I have reflected on this and re-read @dazlich 's initial comment, I doubt that there are any currently-available mapping options in CMEPS that will conserve, so I don't think it's worth @gustavo-marques doing any quick tests right now.

Based on this, my recommendation is as follows; do others agree with this?

• In the short-term we disable runtime generation of runoff mappings, since they currently will not conserve. Instead, trying to make these 'unset' would generate an error.
• In the longer-term we could consider implementing nearest neighbor runtime mapping generation, taking the above points (1) and (2) into account. Whether this is worthwhile depends on whether there are configurations where people will want nearest neighbor without smoothing.

[dazlich]

@billsacks - regarding your point 2, I should have mentioned that after I generated the weights I multiplied them by the src area and divided by the dst area, accomplishing the transformation to volumetric quantities for the remapping.

[billsacks]

Thanks for the clarification, @dazlich - sounds great.

I also asked Bob Oehmke:

(1) Do you agree with this intuition that a transformation is needed to conserve fluxes in a nearest-neighbor remapping?

(2) Are there built-in methods in ESMF to achieve a conservative nearest-neighbor mapping like this, or do we need to do the pre and post normalizations in user code? (I didn't see any in my first pass of the documentation, but may be missing some way to do this.)

He replied:

1. Yes, I do agree that a transformation is needed to conserve fluxes using the nearest-neighbor.

2. We don’t currently have a conservative nearest-neighbor option in ESMF. Thinking about it briefly, I don’t think that it would be very arduous to add one, but I think that it would make sense for you to try it out first for this specific situation and see how it goes. If it works, then we can see about adding it to ESMF in a more general form to handle the various grid options, etc. If there’s a reason that it would be easier to have it in ESMF at first, then let me know and I can take a closer look at adding it.