split explicit corrections

[mark-petersen]

This PR, transplanted from MPAS-Dev/MPAS-Model#728, tests the following issues in the Split Explicit time-stepping routine:

1. Re-split baroclinic/barotropic velocity at the beginning of every time step.
   Previously, we carried the barotropic velocity from one timestep to the next, rather than re-splitting.
2. Set the new barotropic velocity to ubar = F/H
   Currently, the new barotropic velocity is computed from the average over the subcycles.
3. Subtract barotropic average of baroclinic velocity from the new momentum velocity:
   u = ubar + u' - btrAvg( u' )
   Previously, u = ubar + u' and btrAvg( u' ) was not subtracted.
4. Set transport velocity (normalTransportVelocity) to:
   uTr = F/H + u' - btrAvg( u' ) + uBolus - btrAvg( uBolus )
   This is a simplification to remove the velocity correction term. Algebraically it is identical to before.

here:
F is the accumulated barotropic flux
H is the total column depth, including SSH, from the last time known.
u' is the baroclinic velocity from stage 1
ubar is the barotropic velocity

Previously, normalVelocity could deviate from normalTransportVelocity and become noisy. Since normalVelocity does not feed back into the pressure gradient, there was no control on this noise. The previous version also used the velocityCorrection variable, which was extremely noisy, but also misleading, because it was really just the difference between normalVelocity and normalTransportVelocity. There is no need for velocityCorrection when the velocities are built directly, as shown above.

[dengwirda]

@mark-petersen, as part of revisiting the mode-split integrators here, I'm also interested in how our approach compares with some of the current "better" algorithms:

• In the case of nonlinear thickness fluxes (eg. upwind cases for wetting+drying and Lagrangian layers), Hallberg and Adcroft suggest a locally-iterative reconciliation process: Reconciling estimates of the free surface height in Lagrangian vertical coordinate ocean models with mode-split time stepping, Ocn. Mod., 2009.
• I believe the ROMS + NEMO communities (strongly) favour the kind of forward-backward scheme of Shchepetkin and McWilliams: The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model, Ocn. Mod., 2005.

Also, just to migrate some of our slack conversation yesterday:

• For testing, I think we should do detailed RK4 vs split-explicit comparisons in a full ocn spin-up, with as many output variables active as possible. Though it will be expensive (for RK4), I think a comparison after >= 1 yr would be best, to allow enough spin-up for baroclinic signals to develop.
• Since the changes look to impact normalVelocity (and normalTransportVelocity) in general, I'm still not clear on how the differences in model behaviour are restricted to vertVelocityTop alone, since quite a few things (kinetic-energy, rel.-vorticity, divergence, etc) are computed using the same normalVelocity in mpas_ocn_diagnostics and these things drive the momentum tendencies?

[dengwirda]

@mark-petersen I did some initial testing of ideas 1-4 above but unfortunately am not seeing the best outcomes at present:

• I'm restarting at year 4 in a QU30km MPAS-O standalone spin-up (in which the 'standard' split-explicit was used).
• Running 30 days of RK4 cleans up the grid-scale noise in vertVelocityTop and maintains the spun-up kinetic-energy.
• Running 90 days of split-explicit with new 1+2 active does seem to improve the noise in vertVelocityTop, but at the expense of a precipitous reduction in kinetic-energy + circulation, so I believe this can't be quite right.
  
  In the 'standard' split-explicit code, I also noticed that the normalVelocityCorrection seems to be applied only to normalTransportVelocity and not normalVelocity?
  I'm not sure I currently understand this --- my current thinking is that normalVelocity and normalTransportVelocity should only differ by normalGMBolusVelocity; that is normalVelocity = normalTransportVelocity when GM = .false.?