Add general coordinate ability to MOC analysis member

[mark-petersen]

Currently, the MOC analysis member is written for z-level coordinates, i.e. it sums horizontally by k levels. When this MOC is computed for any tilted coordinate (sigma-z, for example) it sums different depths together in zonal bins, which is incorrect.

I am numbering these comments, in case people want to reference them below:

1. I think that the only alteration required for general vertical coordinates is to interpolate these variables:

normalVelocity
layerThickness
vertVelocityTop
vertGMBolusVelocityTop

to a reference z-grid just before horizontal summations are taken, here:

E3SM/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F

Lines 506 to 507 in 3a87fa1

	* normalVelocity(k,iEdge)*dvEdge(iEdge) &
	* 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2))

E3SM/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F

Lines 524 to 525 in 3a87fa1

	sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + (vertVelocityTop(k, iCell) * &
	areaCell(iCell) * regionCellMasks(currentRegion, iCell))

E3SM/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F

Lines 605 to 606 in 417b03f

	* normalVelocity(k,iEdge)*dvEdge(iEdge) &
	* 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2))

E3SM/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F

Lines 623 to 624 in 3a87fa1

	sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + (vertGMBolusVelocityTop(k, iCell) * &
	areaCell(iCell) * regionCellMasks(currentRegion, iCell))

2. The reference z-grid could be refZMid and refBottomDepth, and then everything would work in MPAS-Analysis. It could also be a higher-resolution z-grid for this specific purpose. I propose to use the first for simplicity.

3. For background, you can compute a streamfunction psi(y,z) by integrating v=-d psi/dz vertically or w = d psi/dy horizontally, where y here is latitude. In the Fortran code we integrate v=-d psi/dz vertically along the southern transect in line 506 to get the boundary condition. Then we integrate w = d psi/dy in line 524 horizontally by latitude bins (this produces smoother statistics than integrating many transects vertically). Then we repeat that process for the GM Bolus velocity in lines 605 and 624. I know I've written that in latex a few times, can't find it right now...

4. I would add the general-coordinate MOC as a separate subroutine, chosen by a flag, so our current standard z-level runs don't spend time on unneeded interpolation, and so I do not alter z-level code that is trusted.

5. The current subroutine computes for every edge and cell and multiplies by regionCellMasks. I would change it to check the mask first, and only interpolate and sum if the mask = 1.

6. I am seeking advice on the interpolation method. The only one currently in MPAS-Ocean is the linear ocn_init_interpolation_linear_vert. We could add @dengwirda's https://github.com/dengwirda/PPR here, or perhaps @cbegeman already has one working on a general vertical coordinate branch, and we can merge that one early to use here as well.

7. Related, we need to decide what order of interpolation is sufficient.

8. When interpolating horizontal and vertical velocity, is it important to exactly conserve the transport per edge-column and cell-column, respectively? Perhaps exact conservation is not important for a basin-wide diagnostic (as opposed to numerical methods for advection, say). I'm not sure this question even makes sense for interpolating vertical velocity vertically.

[mark-petersen]

Any feedback is appreciated. Once we agree on the plan, the alterations should be straightforward. Pinging @xylar @maltrud, @vanroekel, @dengwirda @cbegeman @alicebarthel, @milenaveneziani

[cbegeman]

@mark-petersen Thanks for your efforts on this!

6. We have been using @dengwirda's PPR library for our vertical Lagrangian-remapping application. You could checkout my branch https://github.com/cbegeman/E3SM/tree/ocn/add-vertical-remap-feature (Add vertical remap feature #6) or I can quickly create another that just has the PPR library added.
7. With this library it is easy to set what order of interpolation you want, so you might just test a few options (2nd, 3rd, 5th order).
8. PPR is a conservative interpolation so you might not have to worry too much about this. Maybe @dengwirda has thoughts on this.

[xylar]

@mark-petersen, I spent a lot of time working on the overturning streamfunction (OSF) for ISOMIP+, which has a vertical coordinate that is very different from z-level (because of the ice-shelf draft and Haney number constraints).

I found I didn't get correct results when I tried to interpolate vertical velocity in the vertical from the z* to the z-level reference grid. I believe the problem is that assuming that the vertical velocity varies linearly within a cell was not consistent with the actual fluxes in and out of the cell. I think that computing an accurate vertical transport at z-level interfaces would likely require bisecting cells with z-level interfaces and summing the know transports for all other faces of the bisected cell to determine the transport through the z-level interface. This is likely too computationally intensive to be feasible.

Instead, I found that I needed to do conservative interpolation to get the total horizontal transport across faces (a rectangle defined by cell edges and reference-grid layer thickness) on the z-level grid. I then integrate these vertically to get psi at reference z-level interfaces and then bin them to get the OSF.

Here is the python code I use:
https://github.com/MPAS-Dev/compass/blob/master/compass/ocean/tests/isomip_plus/streamfunction.py#L106-L153

I also helped the ROMS group use the same technique to get their OSF for ISOMIP+.

The meshes used for ISOMIP+ are from planar_hex so they are aligned in a way that maybe leads to a cleaner OSF than we can expect for spherical JIGSAW meshes. But I think the technique will still work.