Convert GM-Redi to triad formulation

[mark-petersen]

Currently, the Redi along-isopycnal mixing formulation used simple spatial averaging to compute the isopycnal slope. This PR introduces triad averaging, which uses a very specific form of an expansion of the slope using three cells, the thermal expansion and saline contraction coefficients.

[mark-petersen]

Here is a quick overview. Redi mixing is a mesoscale eddy parameterization that mixes scalars across isopycnal surfaces. The formulation is:

where phi is the tracer, S is the isopycnal slope, K33 is square of S.

The current formulation in MPAS-Ocean was written 5 years ago, has never been used. The original design document is here: https://github.com/MPAS-Dev/MPAS-Documents/tree/master/ocean/gm. The current formulation produces spurious oscillations. See the test tracers on the bottom right:

This noise problem was described in detail by Beckers et al. 1998 who say:

Very disappointing is the result that even for constant slopes, uniform grid spacing, and constant diffusion coefficient ... no well-behaved scheme can be found.

and

... our theorem shows that a directional diffusion is never both consistent and min–max satisfying, however small the slopes of the rotating.

There is some help, though. Griffies et al 1998 introduces a "triad" expansion of the slope calculation, essentially his (30-34) that is total variance diminishing, and much better behaved for a nonlinear equation of state. It was used in MOM and POP. Though noise can still exist, it is much reduced. Specifically, the cross-terms produce over/undershoots beside horizontal boundaries, and the triad formulation does not help that. In practice, the vertical mixing smooths most of this out, but any Redi formulation is not locally bounds preserving.

[mark-petersen]

Here is an explanation of the triads. The horizontal flux is computed over four triads:

Looking at a single triad, the density is expanded in temperature and salinity, and spatial derivatives are taken at each side edge and top boundary:

Griffies et al 1998 defines the slope of the triad here:

and linear expansion of density here:

where the superscripts are the triad corners, and subscripts are the legs of the triad.

[mark-petersen]

@vanroekel or anyone else, the flags are a bit jumbled for GM and Redi. I've updated to move Redi into a separate list, since it is really a separate term on the tracer equations:

CURRENT

&mesoscale_eddy_parameterization
    config_use_standardGM = .false.
    config_gm_lat_variable_c2 = .false.
    config_gm_kappa_lat_depth_variable = .false.
    config_gm_min_stratification_ratio = 0.1
    config_gm_min_phase_speed = 0.1
    config_use_Redi_surface_layer_tapering = .false.
    config_Redi_surface_layer_tapering_extent = 0.0
    config_use_Redi_bottom_layer_tapering = .false.
    config_Redi_bottom_layer_tapering_depth = 0.0
    config_GM_Bolus_kappa_function = 'constant'
    config_standardGM_tracer_kappa = 600.0
    config_GM_Bolus_kappa_min = 0.0
    config_GM_Bolus_kappa_max = 600.0
    config_GM_Bolus_cell_size_min = 20e3
    config_GM_Bolus_cell_size_max = 30e3
    config_Redi_kappa = 0.0
    config_gravWaveSpeed_trunc = 0.3
    config_max_relative_slope = 0.01
/

REVISED, FOR THIS PR

&Redi_isopycnal_mixing
    config_use_Redi = .false.
    config_Redi_kappa_function = 'constant'
    config_Redi_kappa = 0.0
    config_Redi_surface_taper_layers = 3
    config_Redi_bottom_taper_height = 500.0
/
&GM_eddy_parameterization
    config_use_GM = .false.
    config_GM_kappa_function = 'constant'
    config_GM_tracer_kappa = 600.0
    config_GM_kappa_min = 0.0
    config_GM_kappa_max = 600.0
    config_GM_lat_variable_c2 = .false.
    config_GM_kappa_lat_depth_variable = .false.
    config_GM_min_stratification_ratio = 0.1
    config_GM_min_phase_speed = 0.1
    config_GM_cell_size_min = 20e3
    config_GM_cell_size_max = 30e3
    config_gravWaveSpeed_trunc = 0.3
    config_max_relative_slope = 0.01
/

where the new flag config_Redi_kappa_function may be 'constant', 'data', 'match_GM'.

Let me know if anyone has comments on any flag names on this list. They are easy to change now.

[mark-petersen]

Currently testing with analytic functions in a periodic Cartesian domain. See symbolic algebra in python of test functions and solutions for each Redi term here:
https://gist.github.com/mark-petersen/8358a67bb0dd760a77a94dd50ff451c4

[mark-petersen]

@maltrud @milenaveneziani @vanroekel I'm putting these flags in now. Let me know if you have any further suggestions.

&Redi_isopycnal_mixing
    config_use_Redi = .false.
    config_Redi_kappa = 600.0
    config_Redi_closure = 'constant', 'data'
    config_Redi_surface_taper_layers = 10
    config_Redi_bottom_taper_height = 200.0
    config_Redi_resolution_taper = 'ramp'
    config_Redi_resolution_ramp_min = 20e3
    config_Redi_resolution_ramp_max = 30e3
/
&GM_eddy_parameterization
    config_use_GM = .false.
    config_GM_kappa = 600.0
    config_GM_closure = 'N2_dependant'  'constant',
    config_GM_min_stratification_ratio = 0.1
    config_GM_constant_gravWaveSpeed = 0.3
    config_GM_Redi_max_slope = 0.01
    config_GM_resolution_taper = 'ramp'
    config_GM_resolution_ramp_min = 20e3
    config_GM_resolution_ramp_max = 30e3

[milenaveneziani]

If we don't think we will ever have the kappas, res_min, and res_max different between GM and Redi, wouldn't it be clearer to have one section called 'eddy_parameterization' and then have the GM_Redi common config options, like we talked about in the hallway?

[mark-petersen]

I was on the fence. I like your idea of a separate namelist section that applies to all of them better. I'll post in a sec...

[mark-petersen]

Yes, this is better. Thanks for the suggestion:

&Redi_isopycnal_mixing
    config_use_Redi = .false.
    config_Redi_kappa = 600.0
    config_Redi_closure = 'constant'
    config_Redi_surface_taper_layers = 10
    config_Redi_bottom_taper_height = 200.0
/
&GM_eddy_parameterization
    config_use_GM = .false.
    config_GM_kappa = 600.0
    config_GM_closure = 'N2_dependant'
    config_GM_min_stratification_ratio = 0.1
    config_GM_constant_gravWaveSpeed = 0.3
/
&all_eddy_parameterizations
    config_eddy_max_isopycnal_slope = 0.01
    config_eddying_resolution_taper = 'ramp'
    config_eddying_resolution_ramp_min = 20e3
    config_eddying_resolution_ramp_max = 30e3
/

[maltrud]

@mark-petersen the main problem i have with this is "all_eddy_parameterizations". when we add in, say, a submesoscale eddy parameterization this will be confusing. i suggest "GM_and_Redi_parameterizations" or something like that.

[vanroekel]

I like @maltrud's idea. I would suggest "mesoscale_eddy_parameterizations" with the hope we one day move beyond GM/Redi.

[mark-petersen]

OK, currently at

&mesoscale_eddy_parameterizations
    config_MEP_isopycnal_slope_limit = 0.01
    config_eddying_resolution_taper = 'ramp'
    config_eddying_resolution_ramp_min = 20e3
    config_eddying_resolution_ramp_max = 30e3
/

Note altered name config_MEP_isopycnal_slope_limit. I'm still open to feedback.

[maltrud]

i think i'd prefer config_mesoscale_eddy_isopycnal_slope_limit

[mark-petersen]

Agreed, thanks.

[mark-petersen]

@vanroekel, please approve, based on our conversation and the above testing. Also see overview summary at:
E3SM-Project/E3SM#3391 (comment)

[vanroekel]

@mark-petersen I did a quick visual pass and most looks good. I think we should remove the new debug tracers and some of the temporary changes for debugging prior to merging though. I tried to comment where they are.

[vanroekel]

@mark-petersen do we really want to merge these debug tracer changes into ocean/develop?

[mark-petersen]

I would like to merge in the reset of debug tracers now, because we will be conducting tests in the next few weeks where this makes it very easy to control with flags. After we are satisfied, we can remove this section when we turn on Redi mixing by default, with the tuned coefficients. Also, with this merge the debug reset will be on the trunk once, so it is easy to find later.

[vanroekel]

Good point. I agree this would be good for testing.

[mark-petersen]

@vanroekel Please approve when you are ready, and I will merge and update the E3SM PR.

[vanroekel]

approved based on visual inspection and @mark-petersen's extensive testing.

[vanroekel]

Good point. I agree this would be good for testing.

[mark-petersen]

Merged locally, tested with nightly regression suite one last time, both debug and optimized.