Add vertical remap feature

[cbegeman]

This PR adds the option for vertical Lagrangian remapping (VLR). The idea behind VLR is to allow the vertical interfaces between layers to evolve in a Lagrangian fashion, bypassing the Eulerian vertical advection scheme and removing the vertical CFL limit. The layer interfaces are then remapped to a target grid every interval of timesteps. Griffies et al. (2020) provides an overview of VLR.

This development includes the following features:

• A new vertical remapping module mpas_ocn_vertical_remap.F
• Adds PPR library as a submodule https://github.com/dengwirda/PPR, used to perform conservative remapping
• Config options for remapping: the order of remapping, the remapping interval, the remapping limiter
• A diagnostic output file containing remapping-related variables
• Computation of the Eulerian vertical velocity, which is nonzero at timesteps at which regridding and remapping occur
• Vertical tracer advection budget terms must be computed after this Eulerian vertical velocity is computed. Thus, for the remapping case this computation is performed (as usual) in ocn_tracer_advection_mono_tend, where the vertical terms are zero, and repeated in diagnostics after remapping, where these terms may be nonzero.
• As noted below, roughly 20% of runtime is spent remapping when remapping occurs every timestep. However, in most realistic cases remapping at intervals of 3-5 timesteps are acceptable, reducing performance decreases to ~5%.
• VLR is currently only supported for split-explicit time integrator
• Regridding is hard-coded as z-star. Support for alternate grids is planned.

[cbegeman]

For the baroclinic channel test case, vertical remap solution is very similar to the flux-form vertical advection scheme after 12h, with 3-5% increase in run time on 4 cores:

• max difference in sa field (constant): 1e-13
• max difference in pt field: 1e-2
• max difference in uNormal field: 1e-2
• max difference in ssh field: 1e-4 (right panel below)
• left panel shows the diasurface velocity which should be 0, w=-div(hu)+dh_lagrangian computed in diagnostics
  

[xylar]

@cbegeman, would it be helpful if I update E3SM-Ocean-Discussion/master to match E3SM/master so you don't see a bunch of extraneous commits here?

[xylar]

I changed the base (target) of the PR to alternate and changed it back to master, and it looks fine now.

[cbegeman]

@xylar Great, thanks!

[cbegeman]

Tracer budgets were tested with the isomip_plus/5km/z-star/Ocean0 case. Shown below are temperature budget terms after 1 day integration. Note that only the black terms corresponding to activeTracerVerticalAdvectionTendency and activeTracerVerticalAdvectionTopFlux are computed differently in the VLR case.

[cbegeman]

For both baroclinic_channel/4km/rpe_test and isomip_plus/2km/z-star/Ocean0/performance cases, 90% of the remap time was spent in the rmap_1d ppr call. Both cases spend 20% of total simulation time remapping, when remapping occurs every timestep.

This suggests future performance work should focus on ppr routines themselves and/or exploring batch inputs to rmap_1d (currently there is one call to rmap_1d per cell/edge).

[xylar]

@cbegeman, in terms of testing fraction of the runtime use in remapping, the pefromance test case from ISOMIP+ might actually be a little too short. It might be worth running more time steps (say 10x the default) and seeing if the remapping is still a tiny fraction of the runtime. I could imagine that might account for some of the difference between the RPE test and ISOMIP+ -- a smaller fraction of the total time is spent just in initialization and such when you run many more time steps.

[cbegeman]

@cbegeman, in terms of testing fraction of the runtime use in remapping, the pefromance test case from ISOMIP+ might actually be a little too short. It might be worth running more time steps (say 10x the default) and seeing if the remapping is still a tiny fraction of the runtime. I could imagine that might account for some of the difference between the RPE test and ISOMIP+ -- a smaller fraction of the total time is spent just in initialization and such when you run many more time steps.

Thanks for pointing that out. I ran it for 10 days and remapping now takes about 20% of the run time, consistent with the other test case.

[cbegeman]

Following discussion with @dengwirda and @xylar we have decided not to support VLR for RK4.

VLR for RK4, implemented here, has demonstrated instabilities after a few timesteps for both baroclinic_channel and internal_wave test cases.