src/parameterizations/lateral/MOM_MEKE.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		#include "do_concurrent_compat.h"
6
7		!> Implements the Mesoscale Eddy Kinetic Energy framework
8		!! with topographic beta effect included in computing beta in Rhines scale
9
10		module MOM_MEKE
11
12		use iso_fortran_env, only : real32
13
14		use MOM_coms, only : PE_here
15		use MOM_database_comms, only : dbclient_type, dbcomms_CS_type
16		use MOM_debugging, only : hchksum, uvchksum
17		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
18		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
19		use MOM_diag_mediator, only : diag_ctrl, time_type
20		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
21		use MOM_domains, only : pass_vector, pass_var
22		use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, MOM_mesg, is_root_pe
23		use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
24		use MOM_grid, only : ocean_grid_type
25		use MOM_hor_index, only : hor_index_type
26		use MOM_interface_heights, only : find_eta
27		use MOM_interpolate, only : init_external_field, time_interp_external
28		use MOM_interpolate, only : time_interp_external_init
29		use MOM_interpolate, only : external_field
30		use MOM_io, only : vardesc, var_desc, slasher
31		use MOM_isopycnal_slopes, only : calc_isoneutral_slopes
32		use MOM_restart, only : MOM_restart_CS, register_restart_field, query_initialized
33		use MOM_string_functions, only : lowercase
34		use MOM_time_manager, only : time_type_to_real
35		use MOM_unit_scaling, only : unit_scale_type
36		use MOM_variables, only : vertvisc_type, thermo_var_ptrs
37		use MOM_verticalGrid, only : verticalGrid_type
38		use MOM_MEKE_types, only : MEKE_type
39
40		implicit none ; private
41
42		#include <MOM_memory.h>
43
44		public step_forward_MEKE, MEKE_init, MEKE_alloc_register_restart, MEKE_end
45
46		! Constants for this module
47		integer, parameter :: NUM_FEATURES = 4 !< How many features used to predict EKE
48		integer, parameter :: MKE_IDX = 1 !< Index of mean kinetic energy in the feature array
49		integer, parameter :: SLOPE_Z_IDX = 2 !< Index of vertically averaged isopycnal slope in the feature array
50		integer, parameter :: RV_IDX = 3 !< Index of surface relative vorticity in the feature array
51		integer, parameter :: RD_DX_Z_IDX = 4 !< Index of the radius of deformation over the grid size in the feature array
52
53		integer, parameter :: EKE_PROG = 1 !< Use prognostic equation to calculate EKE
54		integer, parameter :: EKE_FILE = 2 !< Read in EKE from a file
55		integer, parameter :: EKE_DBCLIENT = 3 !< Infer EKE using a neural network
56
57		!> Control structure that contains MEKE parameters and diagnostics handles
58		type, public :: MEKE_CS ; private
59		logical :: initialized = .false. !< True if this control structure has been initialized.
60		! Parameters
61		real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim]
62		real :: MEKE_bhFrCoeff!< Efficiency of conversion of ME into MEKE by the biharmonic dissipation [nondim]
63		real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim]
64		real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim]
65		real :: MEKE_damping !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1].
66		real :: MEKE_Cd_scale !< The ratio of the bottom eddy velocity to the column mean
67		!! eddy velocity, i.e. sqrt(2*MEKE), [nondim]. This should be less than 1
68		!! to account for the surface intensification of MEKE.
69		real :: MEKE_Cb !< Coefficient in the \f$\gamma_{bot}\f$ expression [nondim]
70		real :: MEKE_min_gamma!< Minimum value of gamma_b^2 allowed [nondim]
71		real :: MEKE_Ct !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim]
72		logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag.
73		logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC
74		!! framework (Marshall et al., 2012)
75		real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of the
76		!! GEOMETRIC thickness diffusion [nondim].
77		logical :: MEKE_equilibrium_alt !< If true, use an alternative calculation for the
78		!! equilibrium value of MEKE.
79		logical :: MEKE_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value,
80		!! which is calculated at each time step.
81		logical :: GM_src_alt !< If true, use the GM energy conversion form S^2N^2kappa rather
82		!! than the streamfunction for the MEKE GM source term.
83		real :: MEKE_min_depth_tot !< The minimum total thickness over which to distribute MEKE energy
84		!! sources from GM energy conversion [H ~> m or kg m-2]. When the total
85		!! thickness is less than this, the sources are scaled away.
86		logical :: Rd_as_max_scale !< If true the length scale can not exceed the
87		!! first baroclinic deformation radius.
88		logical :: use_old_lscale !< Use the old formula for mixing length scale.
89		logical :: use_min_lscale !< Use simple minimum for mixing length scale.
90		logical :: MEKE_positive !< If true, it guarantees that MEKE will always be >= 0.
91		real :: lscale_maxval !< The ceiling on the MEKE mixing length scale when use_min_lscale is true [L ~> m].
92		real :: cdrag !< The bottom drag coefficient for MEKE, times rescaling factors [H L-1 ~> nondim or kg m-3]
93		real :: MEKE_BGsrc !< Background energy source for MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
94		real :: MEKE_dtScale !< Scale factor to accelerate time-stepping [nondim]
95		real :: MEKE_KhCoeff !< Scaling factor to convert MEKE into Kh [nondim]
96		real :: MEKE_Uscale !< MEKE velocity scale for bottom drag [L T-1 ~> m s-1]
97		real :: MEKE_KH !< Background lateral diffusion of MEKE [L2 T-1 ~> m2 s-1]
98		real :: MEKE_K4 !< Background bi-harmonic diffusivity (of MEKE) [L4 T-1 ~> m4 s-1]
99		real :: KhMEKE_Fac !< A factor relating MEKE%Kh to the diffusivity used for
100		!! MEKE itself [nondim].
101		real :: viscosity_coeff_Ku !< The scaling coefficient in the expression for
102		!! viscosity used to parameterize lateral harmonic momentum mixing
103		!! by unresolved eddies represented by MEKE [nondim].
104		real :: viscosity_coeff_Au !< The scaling coefficient in the expression for
105		!! viscosity used to parameterize lateral biharmonic momentum mixing
106		!! by unresolved eddies represented by MEKE [nondim].
107		real :: Lfixed !< Fixed mixing length scale [L ~> m].
108		real :: aDeform !< Weighting towards deformation scale of mixing length [nondim]
109		real :: aRhines !< Weighting towards Rhines scale of mixing length [nondim]
110		real :: aFrict !< Weighting towards frictional arrest scale of mixing length [nondim]
111		real :: aEady !< Weighting towards Eady scale of mixing length [nondim]
112		real :: aGrid !< Weighting towards grid scale of mixing length [nondim]
113		real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim]
114		real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered
115		!! when computing beta in Rhines scale [nondim]
116		real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its
117		!! equilibrium value [T-1 ~> s-1].
118		logical :: MEKE_advection_bug !< If true, recover a bug in the calculation of the barotropic
119		!! transport for the advection of MEKE, wherein only the transports in the
120		!! deepest layer are used.
121		logical :: fixed_total_depth !< If true, use the nominal bathymetric depth as the estimate of
122		!! the time-varying ocean depth. Otherwise base the depth on the total
123		!! ocean mass per unit area.
124		real :: rho_fixed_total_depth !< A density used to translate the nominal bathymetric depth into an
125		!! estimate of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH
126		!! is true [R ~> kg m-3]
127		logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled.
128		logical :: initialize !< If True, invokes a steady state solver to calculate MEKE.
129		logical :: debug !< If true, write out checksums of data for debugging
130		integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),
131		!! read in from a file, or inferred via a neural network
132		logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure.
133		type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output
134		!>@{ Diagnostic handles
135		integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1
136		integer :: id_src_adv = -1, id_src_mom_K4 = -1, id_src_btm_drag = -1
137		integer :: id_src_GM = -1, id_src_mom_lp = -1, id_src_mom_bh = -1
138		integer :: id_Ub = -1, id_Ut = -1
139		integer :: id_GM_src = -1, id_mom_src = -1, id_mom_src_bh = -1, id_GME_snk = -1, id_decay = -1
140		integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1
141		integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1
142		integer :: id_Lrhines = -1, id_Leady = -1
143		integer :: id_MEKE_equilibrium = -1
144		!>@}
145		type(external_field) :: eke_handle !< Handle for reading in EKE from a file
146		! Infrastructure
147		integer :: id_clock_pass !< Clock for group pass calls
148		type(group_pass_type) :: pass_MEKE !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff
149		type(group_pass_type) :: pass_Kh !< Group halo pass handle for MEKE%Kh, MEKE%Ku, and/or MEKE%Au
150
151		! MEKE via Machine Learning
152		type(dbclient_type), pointer :: client => NULL() !< Pointer to the database client
153
154		logical :: online_analysis !< If true, post the EKE used in MOM6 at every timestep
155		character(len=5) :: model_key = 'mleke' !< Key where the ML-model is stored
156		character(len=7) :: key_suffix !< Suffix appended to every key sent to Redis
157		real :: eke_max !< The maximum value of EKE considered physically reasonable [L2 T-2 ~> m2 s-2]
158
159		! Clock ids
160		integer :: id_client_init !< Clock id to time initialization of the client
161		integer :: id_put_tensor !< Clock id to time put_tensor routine
162		integer :: id_run_model !< Clock id to time running of the ML model
163		integer :: id_unpack_tensor !< Clock id to time retrieval of EKE prediction
164
165		! Diagnostic ids
166		integer :: id_mke = -1 !< Diagnostic id for surface mean kinetic energy
167		integer :: id_slope_z = -1 !< Diagnostic id for vertically averaged horizontal slope magnitude
168		integer :: id_slope_x = -1 !< Diagnostic id for isopycnal slope in the x-direction
169		integer :: id_slope_y = -1 !< Diagnostic id for isopycnal slope in the y-direction
170		integer :: id_rv = -1 !< Diagnostic id for surface relative vorticity
171
172		! Isoneutral blocking parameters
173		integer :: niblock !< The i block size used in calc_isoneutral_slopes [nondim].
174		integer :: njblock !< The j block size used in calc_isoneutral_slopes [nondim].
175		integer :: nkblock !< The k block size used in calc_isoneutral_slopes [nondim].
176
177		end type MEKE_CS
178
179		contains
180
181		!> Integrates forward-in-time the MEKE eddy energy equation.
182		!! See \ref section_MEKE_equations.
183	24	subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, hv, u, v, tv, Time)
184		type(MEKE_type), intent(inout) :: MEKE !< MEKE data.
185		type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
186		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
187		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
188		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
189		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
190		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
191		type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type.
192		real, intent(in) :: dt !< Model(baroclinic) time-step [T ~> s].
193		type(MEKE_CS), intent(inout) :: CS !< MEKE control structure.
194		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: hu !< Accumulated zonal mass flux [H L2 ~> m3 or kg].
195		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: hv !< Accumulated meridional mass flux [H L2 ~> m3 or kg]
196		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
197		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
198		type(thermo_var_ptrs), intent(in) :: tv !< Type containing thermodynamic variables
199		type(time_type), intent(in) :: Time !< The time used for interpolating EKE
200
201		! Local variables
202		real, dimension(SZI_(G),SZJ_(G)) :: &
203	48	data_eke, & ! EKE from file [L2 T-2 ~> m2 s-2]
204	48	mass, & ! The total mass of the water column [R Z ~> kg m-2].
205	48	I_mass, & ! The inverse of mass [R-1 Z-1 ~> m2 kg-1].
206	48	depth_tot, & ! The depth of the water column [H ~> m or kg m-2].
207	48	src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3).
208	48	MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1].
209	48	src_adv, & ! The MEKE source/tendency from the horizontal advection of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
210	48	src_mom_K4, & ! The MEKE source/tendency from the bihamornic of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
211	48	src_btm_drag, & ! The MEKE source/tendency from the bottom drag acting on MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
212	48	src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3).
213	48	src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
214	48	src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
215	48	damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
216	48	MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2].
217	48	drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1]
218	48	drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
219	48	del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2].
220	48	del4MEKE, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
221	48	LmixScale, & ! Eddy mixing length [L ~> m].
222	48	barotrFac2, & ! Ratio of EKE_barotropic / EKE [nondim]
223	48	bottomFac2, & ! Ratio of EKE_bottom / EKE [nondim]
224	48	tmp, & ! Temporary variable for computation of diagnostic velocities [L T-1 ~> m s-1]
225	48	equilibrium_value, & ! The equilibrium value of MEKE to be calculated at
226		! each time step [L2 T-2 ~> m2 s-2]
227	48	damp_rate, & ! The MEKE damping rate [T-1 ~> s-1]
228	48	damping ! The net damping of a field after sdt_damp [nondim]
229
230		real, dimension(SZIB_(G),SZJ_(G)) :: &
231	48	MEKE_uflux, & ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].
232		! In one place, MEKE_uflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
233	48	Kh_u, & ! The zonal diffusivity that is actually used [L2 T-1 ~> m2 s-1].
234	24	baroHu, & ! Depth integrated accumulated zonal mass flux [R Z L2 ~> kg].
235	48	drag_vel_u ! A piston velocity associated with bottom drag at u-points [H T-1 ~> m s-1 or kg m-2 s-1]
236		real, dimension(SZI_(G),SZJB_(G)) :: &
237	48	MEKE_vflux, & ! The meridional advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].
238		! In one place, MEKE_vflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
239	48	Kh_v, & ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1].
240	48	baroHv, & ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg].
241	48	drag_vel_v ! A piston velocity associated with bottom drag at v-points [H T-1 ~> m s-1 or kg m-2 s-1]
242		real :: bh_coeff ! Biharmonic part of efficiency conversion in total MEKE [nondim]
243		real :: Kh_here ! The local horizontal viscosity [L2 T-1 ~> m2 s-1]
244		real :: Inv_Kh_max ! The inverse of the local horizontal viscosity [T L-2 ~> s m-2]
245		real :: K4_here ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1]
246		real :: Inv_K4_max ! The inverse of the local horizontal biharmonic viscosity [T L-4 ~> s m-4]
247		real :: cdrag2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
248		real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]
249		real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].
250		real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate)
251		real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split).
252		real :: damp_step ! Size of damping timestep relative to sdt [nondim]
253		logical :: use_drag_rate ! Flag to indicate drag_rate is finite
254		logical :: any_damping_diags_s1 ! True if any damped diagnostics are enabled in first stage
255		logical :: any_damping_diags ! True if any damped diagnostics are enabled
256		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
257	96	real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features
258		! needed for the machine learning inference, with different
259		! units for the various subarrays [various]
260
261	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
262	24	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
263
264	24	if (.not.CS%initialized) call MOM_error(FATAL, &
265	0	"MOM_MEKE: Module must be initialized before it is used.")
266
267	24	if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then
268	24	use_drag_rate = .true.
269		else
270	0	use_drag_rate = .false.
271		endif
272
273		! Only integrate the MEKE equations if MEKE is required.
274	24	if (.not. allocated(MEKE%MEKE)) then
275		! call MOM_error(FATAL, "MOM_MEKE: MEKE%MEKE is not associated!")
276	0	return
277		endif
278
279		! Map the per-call inputs and the local scratch. The MEKE% and CS% members were mapped
280		! and synced at init; host-written MEKE% members are refreshed below.
281		!$omp target enter data map(to: SN_u, SN_v)
282		!$omp target enter data map(alloc: mass, I_mass, depth_tot, src, drag_rate, drag_rate_visc)
283		!$omp target enter data map(alloc: drag_vel_u, drag_vel_v, bottomFac2, barotrFac2, LmixScale)
284		!$omp target enter data map(alloc: MEKE_uflux, MEKE_vflux, Kh_u, Kh_v, baroHu, baroHv)
285		!$omp target enter data map(alloc: damp_rate, damping, MEKE_current, damp_rate_s1)
286		!$omp target enter data map(alloc: del2MEKE, del4MEKE, MEKE_decay, equilibrium_value)
287		!$omp target enter data map(alloc: src_adv, src_GM, src_mom_K4, src_btm_drag, src_mom_lp, src_mom_bh)
288
289		! These MEKE members are written on the host by other modules (thickness_diffuse,
290		! hor_visc frictional-work accumulation, VarMix), so refresh the device copies.
291		if (allocated(MEKE%GM_src)) then
292		!$omp target update to(MEKE%GM_src)
293		endif
294		if (allocated(MEKE%mom_src)) then
295		!$omp target update to(MEKE%mom_src)
296		endif
297		if (allocated(MEKE%mom_src_bh)) then
298		!$omp target update to(MEKE%mom_src_bh)
299		endif
300		if (allocated(MEKE%GME_snk)) then
301		!$omp target update to(MEKE%GME_snk)
302		endif
303		if (allocated(MEKE%Rd_dx_h)) then
304		!$omp target update to(MEKE%Rd_dx_h)
305		endif
306		if (allocated(MEKE%Kh_diff)) then
307		!$omp target update to(MEKE%Kh_diff)
308		endif
309
310	48	select case(CS%eke_src)
311		case(EKE_PROG)
312	24	if (CS%debug) then
313	0	if (allocated(MEKE%mom_src)) &
314	0	call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2US%L_to_Z*2)
315	0	if (allocated(MEKE%mom_src_bh)) &
316	0	call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, unscale=US%RZ3_T3_to_W_m2US%L_to_Z*2)
317	0	if (allocated(MEKE%GME_snk)) &
318	0	call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2US%L_to_Z*2)
319	0	if (allocated(MEKE%GM_src)) &
320	0	call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, unscale=US%RZ3_T3_to_W_m2US%L_to_Z*2)
321	0	if (allocated(MEKE%MEKE)) &
322	0	call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, unscale=US%L_T_to_m_s**2)
323		call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, unscale=US%s_to_T, &
324	0	scalar_pair=.true.)
325		!$omp target update from(hu, hv)
326		call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=0, symmetric=.true., &
327	0	unscale=GV%H_to_mUS%L_to_m*2)
328		endif
329
330	24	sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping
331	24	mass_neglect = GV%H_to_RZ * GV%H_subroundoff
332	24	cdrag2 = CS%cdrag**2
333
334		! With a depth-dependent (and possibly strong) damping, it seems
335		! advisable to use Strang splitting between the damping and diffusion.
336	24	damp_step = 1.
337	24	if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5
338	24	sdt_damp = sdt * damp_step
339
340		! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]
341	24	if (CS%MEKE_advection_factor>0.) then
342	0	do concurrent (j=js:je, I=is-1:ie) DO_LOCALITY(local(k))
343	0	baroHu(I,j) = 0.
344	0	do k=1,nz
345	0	baroHu(I,j) = baroHu(I,j) + hu(I,j,k) * GV%H_to_RZ
346		enddo
347		enddo
348	0	do concurrent (J=js-1:je, i=is:ie) DO_LOCALITY(local(k))
349	0	baroHv(i,J) = 0.
350	0	do k=1,nz
351	0	baroHv(i,J) = baroHv(i,J) + hv(i,J,k) * GV%H_to_RZ
352		enddo
353		enddo
354	0	if (CS%MEKE_advection_bug) then
355		! This obviously incorrect code reproduces a bug in the original implementation of
356		! the MEKE advection.
357	0	do concurrent (j=js:je, I=is-1:ie)
358	0	baroHu(I,j) = hu(I,j,nz) * GV%H_to_RZ
359		enddo
360	0	do concurrent (J=js-1:je, i=is:ie)
361	0	baroHv(i,J) = hv(i,J,nz) * GV%H_to_RZ
362		enddo
363		endif
364		endif
365
366		! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
367	24	if (CS%visc_drag .and. allocated(visc%Kv_bbl_u) .and. allocated(visc%Kv_bbl_v)) then
368	24	do concurrent (j=js:je, I=is-1:ie)
369	174240	drag_vel_u(I,j) = 0.0
370	174240	if ((G%mask2dCu(I,j) > 0.0) .and. (visc%bbl_thick_u(I,j) > 0.0)) &
371	177168	drag_vel_u(I,j) = visc%Kv_bbl_u(I,j) / visc%bbl_thick_u(I,j)
372		enddo
373	24	do concurrent (J=js-1:je, i=is:ie)
374	175680	drag_vel_v(i,J) = 0.0
375	175680	if ((G%mask2dCv(i,J) > 0.0) .and. (visc%bbl_thick_v(i,J) > 0.0)) &
376	178584	drag_vel_v(i,J) = visc%Kv_bbl_v(i,J) / visc%bbl_thick_v(i,J)
377		enddo
378
379	24	do concurrent (j=js:je, i=is:ie)
380		drag_rate_visc(i,j) = (0.25G%IareaT(i,j) &
381		(((G%areaCu(I-1,j)*drag_vel_u(I-1,j)) + &
382		(G%areaCu(I,j)*drag_vel_u(I,j))) + &
383		((G%areaCv(i,J-1)*drag_vel_v(i,J-1)) + &
384	175704	(G%areaCv(i,J)*drag_vel_v(i,J))) ) )
385		enddo
386		else
387	0	do concurrent (j=js:je, i=is:ie)
388	0	drag_rate_visc(i,j) = 0.
389		enddo
390		endif
391
392	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1) DO_LOCALITY(local(k))
393	181536	mass(i,j) = 0.0
394	13796736	do k=1,nz
395	13796736	mass(i,j) = mass(i,j) + G%mask2dT(i,j) * (GV%H_to_RZ * h(i,j,k)) ! [R Z ~> kg m-2]
396		enddo
397	181536	I_mass(i,j) = 0.0
398	184488	if (mass(i,j) > 0.0) I_mass(i,j) = 1.0 / mass(i,j) ! [R-1 Z-1 ~> m2 kg-1]
399		enddo
400
401	24	if (CS%fixed_total_depth) then
402	24	if (GV%Boussinesq) then
403	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
404	184488	depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * GV%Z_to_H
405		enddo
406		else
407	0	do concurrent (j=js-1:je+1, i=is-1:ie+1)
408	0	depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * CS%rho_fixed_total_depth * GV%RZ_to_H
409		enddo
410		endif
411		else
412	0	do concurrent (j=js-1:je+1, i=is-1:ie+1)
413	0	depth_tot(i,j) = mass(i,j) * GV%RZ_to_H
414		enddo
415		endif
416
417	24	if (CS%initialize) then
418		! MEKE_equilibrium runs on the host, once: sync its device-computed inputs down and
419		! its result (MEKE%MEKE) back up.
420		!$omp target update from(drag_rate_visc, I_mass, depth_tot)
421	1	call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)
422	1	CS%initialize = .false.
423		!$omp target update to(MEKE%MEKE)
424		endif
425
426		! TODO: THIS RUNS ON THE HOST, BE WARY! AAAAAA
427		! MEKE_lengthScales_0d takes (1+x)0.8 / (1+x)0.25 and the intrinsic pow differs
428		! CPU vs GPU in the last bit, so this stays on the host (with the depth_tot /
429		! bottomFac2 / barotrFac2 / LmixScale syncs below) until a bit-reproducible pow
430		! lands. Only bottomFac2 (the **0.8) actually diverged on device.
431		!$omp target update from(depth_tot)
432	24	call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
433		!$omp target update to(bottomFac2, barotrFac2, LmixScale)
434	24	if (CS%debug) then
435		!$omp target update from(mass, drag_rate_visc, bottomFac2, barotrFac2, LmixScale)
436		if (CS%visc_drag) then
437		!$omp target update from(drag_vel_u, drag_vel_v)
438		endif
439	0	if (CS%visc_drag) &
440		call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, &
441	0	unscale=GV%H_to_mks*US%s_to_T, scalar_pair=.true.)
442	0	call hchksum(mass, 'MEKE mass',G%HI,haloshift=1, unscale=US%RZ_to_kg_m2)
443	0	call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', G%HI, unscale=GV%H_to_mks*US%s_to_T)
444	0	call hchksum(bottomFac2, 'MEKE bottomFac2', G%HI)
445	0	call hchksum(barotrFac2, 'MEKE barotrFac2', G%HI)
446	0	call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m)
447		endif
448
449	24	if (allocated(MEKE%Le)) then
450	0	do concurrent (j=js:je, i=is:ie)
451	0	MEKE%Le(i,j) = LmixScale(i,j)
452		enddo
453		endif
454
455		! Aggregate sources of MEKE (background, frictional and GM)
456	24	do concurrent (j=js:je, i=is:ie)
457	175704	src(i,j) = CS%MEKE_BGsrc
458		enddo
459
460		! Initialize diagnostics
461	24	if (CS%id_src_adv > 0) then
462	0	do concurrent (j=js:je, i=is:ie)
463	0	src_adv(i,j) = 0.
464		enddo
465		endif
466	24	if (CS%id_src_GM > 0) then
467	0	do concurrent (j=js:je, i=is:ie)
468	0	src_GM(i,j) = 0.
469		enddo
470		endif
471	24	if (CS%id_src_mom_lp > 0) then
472	0	do concurrent (j=js:je, i=is:ie)
473	0	src_mom_lp(i,j) = 0.
474		enddo
475		endif
476	24	if (CS%id_src_mom_bh > 0) then
477	0	do concurrent (j=js:je, i=is:ie)
478	0	src_mom_bh(i,j) = 0.
479		enddo
480		endif
481	24	if (CS%id_src_mom_K4 > 0) then
482	0	do concurrent (j=js:je, i=is:ie)
483	0	src_mom_K4(i,j) = 0.
484		enddo
485		endif
486	24	if (CS%id_src_btm_drag > 0) then
487	0	do concurrent (j=js:je, i=is:ie)
488	0	src_btm_drag(i,j) = 0.
489		enddo
490		endif
491
492		! Identify any damped diagnostics in first stage of Strang splitting
493		any_damping_diags_s1 = any([ &
494		CS%id_src_GM > 0, &
495		CS%id_src_mom_lp > 0, &
496		CS%id_src_mom_bh > 0, &
497		CS%id_src_btm_drag > 0 &
498	120	])
499
500		! Identify any damped diagnostics
501		any_damping_diags = any([ &
502		any_damping_diags_s1, &
503		CS%id_src_adv > 0, &
504		CS%id_src_mom_K4 > 0 &
505	96	])
506
507	24	if (CS%MEKE_FrCoeff > 0.) then
508	0	do concurrent (j=js:je, i=is:ie)
509	0	src(i,j) = src(i,j) - CS%MEKE_FrCoeff * I_mass(i,j) * MEKE%mom_src(i,j)
510		enddo
511		endif
512
513	24	if (allocated(MEKE%mom_src_bh)) then
514	0	if (CS%MEKE_bhFrCoeff > 0. .and. CS%MEKE_FrCoeff > 0.) then
515	0	bh_coeff = CS%MEKE_bhFrCoeff - CS%MEKE_FrCoeff
516		else
517	0	bh_coeff = CS%MEKE_bhFrCoeff
518		endif
519
520	0	do concurrent (j=js:je, i=is:ie)
521	0	src(i,j) = src(i,j) - bh_coeff * I_mass(i,j) * MEKE%mom_src_bh(i,j)
522		enddo
523
524	0	if (CS%id_src_mom_lp > 0) then
525	0	do concurrent (j=js:je, i=is:ie)
526		src_mom_lp(i,j) = -CS%MEKE_FrCoeff * I_mass(i,j) &
527	0	* (MEKE%mom_src(i,j) - MEKE%mom_src_bh(i,j))
528		enddo
529		endif
530
531	0	if (CS%id_src_mom_bh > 0) then
532	0	do concurrent (j=js:je, i=is:ie)
533	0	src_mom_bh(i,j) = -CS%MEKE_bhFrCoeff * I_mass(i,j) * MEKE%mom_src_bh(i,j)
534		enddo
535		endif
536		endif
537
538	24	if (allocated(MEKE%GME_snk)) then
539	0	do concurrent (j=js:je, i=is:ie)
540	0	src(i,j) = src(i,j) - CS%MEKE_GMECoeffI_mass(i,j)MEKE%GME_snk(i,j)
541		enddo
542		endif
543
544	24	if (allocated(MEKE%GM_src)) then
545	24	if (CS%GM_src_alt) then
546	0	do concurrent (j=js:je, i=is:ie)
547		src(i,j) = src(i,j) - CS%MEKE_GMcoeff*MEKE%GM_src(i,j) / &
548	0	(GV%H_to_RZ * MAX(CS%MEKE_min_depth_tot, depth_tot(i,j)))
549		enddo
550		else
551	24	do concurrent (j=js:je, i=is:ie)
552	175704	src(i,j) = src(i,j) - CS%MEKE_GMcoeffI_mass(i,j)MEKE%GM_src(i,j)
553		enddo
554
555	24	do concurrent (j=js:je, i=is:ie)
556	175704	src_GM(i,j) = -CS%MEKE_GMcoeffI_mass(i,j)MEKE%GM_src(i,j)
557		enddo
558		endif
559		endif
560
561	24	if (CS%MEKE_equilibrium_restoring) then
562		! MEKE_equilibrium_restoring runs on the host: sync its device-computed input down
563		! and its result up before the device loop consumes it.
564		!$omp target update from(depth_tot)
565		call MEKE_equilibrium_restoring(CS, G, GV, US, SN_u, SN_v, depth_tot, &
566	0	equilibrium_value)
567		!$omp target update to(equilibrium_value)
568	0	do concurrent (j=js:je, i=is:ie)
569	0	src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - equilibrium_value(i,j))
570		enddo
571		endif
572
573	24	if (CS%debug) then
574		!$omp target update from(src)
575	0	call hchksum(src, "MEKE src", G%HI, haloshift=0, unscale=US%L_to_m*2US%s_to_T**3)
576		endif
577
578		! Increase EKE by a full time-steps worth of source
579	24	do concurrent (j=js:je, i=is:ie)
580	172800	MEKE_current(i,j) = MEKE%MEKE(i,j)
581	175704	MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdtsrc(i,j))G%mask2dT(i,j)
582		enddo
583
584	24	if (use_drag_rate) then
585		! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
586	24	do concurrent (j=js:je, i=is:ie)
587		drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
588	175704	cdrag2 * ( max(0.0, 2.0bottomFac2(i,j)MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
589		enddo
590		else
591	0	do concurrent (j=js:je, i=is:ie)
592	0	drag_rate(i,j) = 0.
593		enddo
594		endif
595
596		! First stage of Strang splitting
597
598	24	do concurrent (j=js:je, i=is:ie)
599	172800	damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
600
601	175704	if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0.
602		! notice that the above line ensures a damping only if MEKE is positive,
603		! while leaving MEKE unchanged if it is negative
604		enddo
605
606		! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the existing
607		! bit-reproducible solution.
608	24	do concurrent (j=js:je, i=is:ie)
609	175704	MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))
610		enddo
611
612	24	if (any_damping_diags_s1) then
613	0	do concurrent (j=js:je, i=is:ie)
614	0	damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))
615		enddo
616
617	0	if (CS%id_decay > 0) then
618	0	do concurrent (j=js:je, i=is:ie)
619	0	MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j)
620		enddo
621		endif
622
623	0	if (CS%id_src_GM > 0) then
624	0	do concurrent (j=js:je, i=is:ie)
625	0	src_GM(i,j) = src_GM(i,j) * damping(i,j)
626		enddo
627		endif
628
629	0	if (CS%id_src_mom_lp > 0) then
630	0	do concurrent (j=js:je, i=is:ie)
631	0	src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)
632		enddo
633		endif
634
635	0	if (CS%id_src_mom_bh > 0) then
636	0	do concurrent (j=js:je, i=is:ie)
637	0	src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)
638		enddo
639		endif
640
641	0	if (CS%id_src_btm_drag > 0) then
642	0	do concurrent (j=js:je, i=is:ie)
643		src_btm_drag(i,j) = -MEKE_current(i,j) * ( &
644		damp_step * (damp_rate(i,j) * damping(i,j)) &
645	0	)
646		enddo
647
648		! Store the effective damping rate if sdt is split
649	0	if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) then
650	0	do concurrent (j=js:je, i=is:ie)
651	0	damp_rate_s1(i,j) = damp_rate(i,j) * damping(i,j)
652		enddo
653		endif
654		endif
655		endif
656
657	24	if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then
658		! Update MEKE in the halos for lateral or bi-harmonic diffusion
659	24	call cpu_clock_begin(CS%id_clock_pass)
660	24	call do_group_pass(CS%pass_MEKE, G%Domain, omp_offload=.true.)
661	24	call cpu_clock_end(CS%id_clock_pass)
662		endif
663
664	24	if (CS%MEKE_K4 >= 0.0) then
665		! Calculate Laplacian of MEKE using MEKE_uflux and MEKE_vflux as temporary work space.
666	0	do concurrent (j=js-1:je+1, I=is-2:ie+1)
667		! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
668		MEKE_uflux(I,j) = (G%dy_Cu(I,j)G%IdxCu_OBCmask(I,j)) &
669	0	(MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
670		! This would have units of [R Z L2 T-2 ~> kg s-2]
671		! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)G%IdxCu(I,j)) &
672		! ((2.0mass(i,j)mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
673		! (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
674		enddo
675	0	do concurrent (J=js-2:je+1, i=is-1:ie+1)
676		! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
677		MEKE_vflux(i,J) = (G%dx_Cv(i,J)G%IdyCv_OBCmask(i,J)) &
678	0	(MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
679		! This would have units of [R Z L2 T-2 ~> kg s-2]
680		! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)G%IdyCv(i,J)) &
681		! ((2.0mass(i,j)mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
682		! (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
683		enddo
684
685	0	do concurrent (j=js-1:je+1, i=is-1:ie+1) ! del2MEKE has units [T-2 ~> s-2].
686		del2MEKE(i,j) = G%IareaT(i,j) * &
687	0	((MEKE_uflux(I,j) - MEKE_uflux(I-1,j)) + (MEKE_vflux(i,J) - MEKE_vflux(i,J-1)))
688		enddo
689
690		! Bi-harmonic diffusion of MEKE
691	0	do concurrent (j=js:je, I=is-1:ie) DO_LOCALITY(local(K4_here, Inv_K4_max))
692	0	K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
693		! Limit Kh to avoid CFL violations.
694		Inv_K4_max = 64.0 * sdt * ((G%dy_Cu(I,j)G%IdxCu(I,j)) &
695	0	max(G%IareaT(i,j), G%IareaT(i+1,j)))**2
696	0	if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max
697
698		! Here the units of MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
699		MEKE_uflux(I,j) = ((K4_here * (G%dy_Cu(I,j)G%IdxCu(I,j))) &
700		((2.0mass(i,j)mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
701	0	(del2MEKE(i+1,j) - del2MEKE(i,j))
702		enddo
703	0	do concurrent (J=js-1:je, i=is:ie) DO_LOCALITY(local(K4_here, Inv_K4_max))
704	0	K4_here = CS%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
705	0	Inv_K4_max = 64.0 * sdt * ((G%dx_Cv(i,J)G%IdyCv(i,J)) max(G%IareaT(i,j), G%IareaT(i,j+1)))**2
706	0	if (K4_here*Inv_K4_max > 0.3) K4_here = 0.3 / Inv_K4_max
707
708		! Here the units of MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
709		MEKE_vflux(i,J) = ((K4_here * (G%dx_Cv(i,J)G%IdyCv(i,J))) &
710		((2.0mass(i,j)mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
711	0	(del2MEKE(i,j+1) - del2MEKE(i,j))
712		enddo
713		! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2 ~> m2 s-2].
714	0	do concurrent (j=js:je, i=is:ie)
715		del4MEKE(i,j) = (sdt(G%IareaT(i,j)I_mass(i,j))) * &
716		((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
717	0	(MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
718		src_mom_K4(i,j) = (G%IareaT(i,j)I_mass(i,j)) &
719		((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
720	0	(MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
721		enddo
722		endif !
723
724	24	if (CS%kh_flux_enabled) then
725		! Lateral diffusion of MEKE
726	24	Kh_here = max(0., CS%MEKE_Kh)
727	24	do concurrent (j=js:je, I=is-1:ie) DO_LOCALITY(local(Inv_Kh_max) local_init(Kh_here))
728		! Limit Kh to avoid CFL violations.
729	174240	if (allocated(MEKE%Kh)) &
730		Kh_here = max(0., CS%MEKE_Kh) + &
731	174240	CS%KhMEKE_Fac0.5(MEKE%Kh(i,j)+MEKE%Kh(i+1,j))
732	174240	if (allocated(MEKE%Kh_diff)) &
733		Kh_here = max(0.,CS%MEKE_Kh) + &
734	0	CS%KhMEKE_Fac0.5(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i+1,j))
735		Inv_Kh_max = 2.0sdt ((G%dy_Cu(I,j)G%IdxCu(I,j)) &
736	174240	max(G%IareaT(i,j),G%IareaT(i+1,j)))
737	174240	if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
738	174240	Kh_u(I,j) = Kh_here
739
740		! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
741		MEKE_uflux(I,j) = ((Kh_here * (G%dy_Cu(I,j)G%IdxCu(I,j))) &
742		((2.0mass(i,j)mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
743	177168	(MEKE%MEKE(i,j) - MEKE%MEKE(i+1,j))
744		enddo
745	24	do concurrent (J=js-1:je, i=is:ie) DO_LOCALITY(local(Inv_Kh_max) local_init(Kh_here))
746	175680	if (allocated(MEKE%Kh)) &
747	175680	Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac * 0.5*(MEKE%Kh(i,j)+MEKE%Kh(i,j+1))
748	175680	if (allocated(MEKE%Kh_diff)) &
749	0	Kh_here = max(0.,CS%MEKE_Kh) + CS%KhMEKE_Fac * 0.5*(MEKE%Kh_diff(i,j)+MEKE%Kh_diff(i,j+1))
750	175680	Inv_Kh_max = 2.0sdt ((G%dx_Cv(i,J)G%IdyCv(i,J)) max(G%IareaT(i,j),G%IareaT(i,j+1)))
751	175680	if (Kh_here*Inv_Kh_max > 0.25) Kh_here = 0.25 / Inv_Kh_max
752	175680	Kh_v(i,J) = Kh_here
753
754		! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
755		MEKE_vflux(i,J) = ((Kh_here * (G%dx_Cv(i,J)G%IdyCv(i,J))) &
756		((2.0mass(i,j)mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
757	178584	(MEKE%MEKE(i,j) - MEKE%MEKE(i,j+1))
758		enddo
759	24	if (CS%MEKE_advection_factor>0.) then
760	0	advFac = CS%MEKE_advection_factor / sdt ! [T-1 ~> s-1]
761	0	do concurrent (j=js:je, I=is-1:ie)
762		! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
763	0	if (baroHu(I,j)>0.) then
764	0	MEKE_uflux(I,j) = MEKE_uflux(I,j) + baroHu(I,j)MEKE%MEKE(i,j)advFac
765	0	elseif (baroHu(I,j)<0.) then
766	0	MEKE_uflux(I,j) = MEKE_uflux(I,j) + baroHu(I,j)MEKE%MEKE(i+1,j)advFac
767		endif
768		enddo
769	0	do concurrent (J=js-1:je, i=is:ie)
770		! Here the units of the quantities added to MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
771	0	if (baroHv(i,J)>0.) then
772	0	MEKE_vflux(i,J) = MEKE_vflux(i,J) + baroHv(i,J)MEKE%MEKE(i,j)advFac
773	0	elseif (baroHv(i,J)<0.) then
774	0	MEKE_vflux(i,J) = MEKE_vflux(i,J) + baroHv(i,J)MEKE%MEKE(i,j+1)advFac
775		endif
776		enddo
777		endif
778
779	24	do concurrent (j=js:je, i=is:ie)
780		MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt(G%IareaT(i,j)I_mass(i,j))) * &
781		((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
782	175704	(MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
783		enddo
784
785	24	if (CS%id_src_adv > 0) then
786	0	do concurrent (j=js:je, i=is:ie)
787		src_adv(i,j) = (G%IareaT(i,j)I_mass(i,j)) &
788		((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + &
789	0	(MEKE_vflux(i,J-1) - MEKE_vflux(i,J)))
790		enddo
791		endif
792		endif ! MEKE_KH>0
793
794		! Add on bi-harmonic tendency
795	24	if (CS%MEKE_K4 >= 0.0) then
796	0	do concurrent (j=js:je, i=is:ie)
797	0	MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + del4MEKE(i,j)
798		enddo
799		endif
800
801		! Second stage of Strang splitting
802	24	if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.0) then
803		! Recalculate the drag rate, since MEKE has changed.
804	0	if (use_drag_rate) then
805	0	do concurrent (j=js:je, i=is:ie)
806		drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
807	0	cdrag2 * ( max(0.0, 2.0bottomFac2(i,j)MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
808		enddo
809		endif
810
811	0	do concurrent (j=js:je, i=is:ie)
812	0	damp_rate(i,j) = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
813
814	0	if (MEKE%MEKE(i,j) < 0.) damp_rate(i,j) = 0.
815		! notice that the above line ensures a damping only if MEKE is positive,
816		! while leaving MEKE unchanged if it is negative
817		enddo
818
819		! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the
820		! existing bit-reproducible solution.
821	0	do concurrent (j=js:je, i=is:ie)
822	0	MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))
823		enddo
824
825	0	if (any_damping_diags) then
826	0	do concurrent (j=js:je, i=is:ie)
827	0	damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))
828		enddo
829
830	0	if (CS%id_decay > 0) then
831	0	do concurrent (j=js:je, i=is:ie)
832	0	MEKE_decay(i,j) = damp_rate(i,j) * G%mask2dT(i,j)
833		enddo
834		endif
835
836	0	if (CS%id_src_GM > 0) then
837	0	do concurrent (j=js:je, i=is:ie)
838	0	src_GM(i,j) = src_GM(i,j) * damping(i,j)
839		enddo
840		endif
841
842	0	if (CS%id_src_mom_lp > 0) then
843	0	do concurrent (j=js:je, i=is:ie)
844	0	src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)
845		enddo
846		endif
847
848	0	if (CS%id_src_mom_bh > 0) then
849	0	do concurrent (j=js:je, i=is:ie)
850	0	src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)
851		enddo
852		endif
853
854	0	if (CS%id_src_adv > 0) then
855	0	do concurrent (j=js:je, i=is:ie)
856	0	src_adv(i,j) = src_adv(i,j) * damping(i,j)
857		enddo
858		endif
859
860	0	if (CS%id_src_mom_K4 > 0) then
861	0	do concurrent (j=js:je, i=is:ie)
862	0	src_mom_K4(i,j) = src_mom_K4(i,j) * damping(i,j)
863		enddo
864		endif
865
866	0	if (CS%id_src_btm_drag > 0) then
867	0	do concurrent (j=js:je, i=is:ie)
868		src_btm_drag(i,j) = -MEKE_current(i,j) * (damp_step &
869		* ((damp_rate(i,j) + damp_rate_s1(i,j)) * damping(i,j)) &
870	0	)
871		enddo
872		endif
873		endif
874		endif ! MEKE_KH>=0
875
876	24	if (CS%debug) then
877		!$omp target update from(MEKE%MEKE)
878	0	call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, unscale=US%L_T_to_m_s**2)
879		endif
880
881		case(EKE_FILE)
882	0	call time_interp_external(CS%eke_handle, Time, data_eke, scale=US%m_s_to_L_T**2)
883	0	do j=js,je ; do i=is,ie
884	0	MEKE%MEKE(i,j) = data_eke(i,j) * G%mask2dT(i,j)
885		enddo ; enddo
886		!$omp target update to(MEKE%MEKE)
887	0	call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
888		!$omp target update to(bottomFac2, barotrFac2, LmixScale)
889		case(EKE_DBCLIENT)
890	0	call pass_vector(u, v, G%Domain)
891	0	call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale)
892		!$omp target update to(bottomFac2, barotrFac2, LmixScale)
893	0	call ML_MEKE_calculate_features(G, GV, US, CS, MEKE%Rd_dx_h, u, v, tv, h, dt, features_array)
894	0	call predict_MEKE(G, US, CS, SIZE(h), Time, features_array, MEKE%MEKE)
895		!$omp target update to(MEKE%MEKE)
896		case default
897	24	call MOM_error(FATAL,"Invalid method specified for calculating EKE")
898		end select
899
900	24	if (CS%MEKE_positive) then
901	0	do concurrent (j=js:je, i=is:ie)
902	0	MEKE%MEKE(i,j) = MAX(0., MEKE%MEKE(i,j))
903		enddo
904		endif
905
906	24	call cpu_clock_begin(CS%id_clock_pass)
907	24	call do_group_pass(CS%pass_MEKE, G%Domain, omp_offload=.true.)
908	24	call cpu_clock_end(CS%id_clock_pass)
909
910		! Calculate diffusivity for main model to use
911	24	if (CS%MEKE_KhCoeff>0.) then
912	24	if (.not.CS%MEKE_GEOMETRIC) then
913	24	if (CS%use_old_lscale) then
914	0	if (CS%Rd_as_max_scale) then
915	0	do concurrent (j=js:je, i=is:ie)
916		MEKE%Kh(i,j) = (CS%MEKE_KhCoeff * &
917		sqrt(2.max(0.,barotrFac2(i,j)MEKE%MEKE(i,j))G%areaT(i,j)) ) &
918	0	min(MEKE%Rd_dx_h(i,j), 1.0)
919		enddo
920		else
921	0	do concurrent (j=js:je, i=is:ie)
922		MEKE%Kh(i,j) = CS%MEKE_KhCoeff * &
923	0	sqrt(2.max(0., barotrFac2(i,j)MEKE%MEKE(i,j))*G%areaT(i,j))
924		enddo
925		endif
926		else
927	24	do concurrent (j=js:je, i=is:ie)
928		MEKE%Kh(i,j) = CS%MEKE_KhCoeff * &
929	175704	sqrt(2.max(0., barotrFac2(i,j)MEKE%MEKE(i,j))) * LmixScale(i,j)
930		enddo
931		endif
932		endif
933		endif
934
935		! Calculate viscosity for the main model to use
936	24	if (CS%viscosity_coeff_Ku /=0.) then
937	0	do concurrent (j=js:je, i=is:ie)
938	0	MEKE%Ku(i,j) = CS%viscosity_coeff_Ku * sqrt(2.max(0.,MEKE%MEKE(i,j))) LmixScale(i,j)
939		enddo
940		endif
941
942	24	if (CS%viscosity_coeff_Au /=0.) then
943	0	do concurrent (j=js:je, i=is:ie)
944	0	MEKE%Au(i,j) = CS%viscosity_coeff_Au * sqrt(2.max(0.,MEKE%MEKE(i,j))) LmixScale(i,j)**3
945		enddo
946		endif
947
948		if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
949	24	.or. allocated(MEKE%Le)) then
950	24	call cpu_clock_begin(CS%id_clock_pass)
951	24	call do_group_pass(CS%pass_Kh, G%Domain, omp_offload=.true.)
952	24	call cpu_clock_end(CS%id_clock_pass)
953		endif
954
955		! Everything below runs on the host: sync the prognostic fields and any posted
956		! scratch down from the device. MEKE%MEKE and MEKE%Kh in particular are read on the
957		! host by thickness_diffuse and the restart machinery.
958		!$omp target update from(MEKE%MEKE)
959		if (allocated(MEKE%Kh)) then
960		!$omp target update from(MEKE%Kh)
961		endif
962		if (allocated(MEKE%Ku)) then
963		!$omp target update from(MEKE%Ku)
964		endif
965		if (allocated(MEKE%Au)) then
966		!$omp target update from(MEKE%Au)
967		endif
968		if (allocated(MEKE%Le)) then
969		!$omp target update from(MEKE%Le)
970		endif
971		if (CS%id_Ub>0 .or. CS%id_gamma_b>0) then
972		!$omp target update from(bottomFac2)
973		endif
974		if (CS%id_Ut>0 .or. CS%id_gamma_t>0) then
975		!$omp target update from(barotrFac2)
976		endif
977		if (CS%id_Le>0) then
978		!$omp target update from(LmixScale)
979		endif
980		if (CS%id_KhMEKE_u>0) then
981		!$omp target update from(Kh_u)
982		endif
983		if (CS%id_KhMEKE_v>0) then
984		!$omp target update from(Kh_v)
985		endif
986		if (CS%id_src>0) then
987		!$omp target update from(src)
988		endif
989		if (CS%id_src_adv>0) then
990		!$omp target update from(src_adv)
991		endif
992		if (CS%id_src_mom_K4>0) then
993		!$omp target update from(src_mom_K4)
994		endif
995		if (CS%id_src_btm_drag>0) then
996		!$omp target update from(src_btm_drag)
997		endif
998		if (CS%id_src_GM>0) then
999		!$omp target update from(src_GM)
1000		endif
1001		if (CS%id_src_mom_lp>0) then
1002		!$omp target update from(src_mom_lp)
1003		endif
1004		if (CS%id_src_mom_bh>0) then
1005		!$omp target update from(src_mom_bh)
1006		endif
1007		if (CS%id_decay>0) then
1008		!$omp target update from(MEKE_decay)
1009		endif
1010
1011		! Offer fields for averaging.
1012	96	if (any([CS%id_Ue, CS%id_Ub, CS%id_Ut] > 0)) &
1013	0	tmp(:,:) = 0.
1014	24	if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag)
1015	24	if (CS%id_Ue>0) then
1016	0	do j=js,je ; do i=is,ie
1017	0	tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j)))
1018		enddo ; enddo
1019	0	call post_data(CS%id_Ue, tmp, CS%diag)
1020		endif
1021	24	if (CS%id_Ub>0) then
1022	0	do j=js,je ; do i=is,ie
1023	0	tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * bottomFac2(i,j)))
1024		enddo ; enddo
1025	0	call post_data(CS%id_Ub, tmp, CS%diag)
1026		endif
1027	24	if (CS%id_Ut>0) then
1028	0	do j=js,je ; do i=is,ie
1029	0	tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * barotrFac2(i,j)))
1030		enddo ; enddo
1031	0	call post_data(CS%id_Ut, tmp, CS%diag)
1032		endif
1033	24	if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
1034	24	if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
1035	24	if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)
1036	24	if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag)
1037	24	if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag)
1038	24	if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag)
1039	24	if (CS%id_src_adv>0) call post_data(CS%id_src_adv, src_adv, CS%diag)
1040	24	if (CS%id_src_mom_K4>0) call post_data(CS%id_src_mom_K4, src_mom_K4, CS%diag)
1041	24	if (CS%id_src_btm_drag>0) call post_data(CS%id_src_btm_drag, src_btm_drag, CS%diag)
1042	24	if (CS%id_src_GM>0) call post_data(CS%id_src_GM, src_GM, CS%diag)
1043	24	if (CS%id_src_mom_lp>0) call post_data(CS%id_src_mom_lp, src_mom_lp, CS%diag)
1044	24	if (CS%id_src_mom_bh>0) call post_data(CS%id_src_mom_bh, src_mom_bh, CS%diag)
1045	24	if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag)
1046	24	if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag)
1047	24	if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag)
1048	24	if (CS%id_mom_src_bh>0) call post_data(CS%id_mom_src_bh, MEKE%mom_src_bh, CS%diag)
1049	24	if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag)
1050	24	if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag)
1051	24	if (CS%id_gamma_b>0) then
1052	0	do j=js,je ; do i=is,ie
1053	0	bottomFac2(i,j) = sqrt(bottomFac2(i,j))
1054		enddo ; enddo
1055	0	call post_data(CS%id_gamma_b, bottomFac2, CS%diag)
1056		endif
1057	24	if (CS%id_gamma_t>0) then
1058	0	do j=js,je ; do i=is,ie
1059	0	barotrFac2(i,j) = sqrt(barotrFac2(i,j))
1060		enddo ; enddo
1061	0	call post_data(CS%id_gamma_t, barotrFac2, CS%diag)
1062		endif
1063
1064		!$omp target exit data map(release: SN_u, SN_v)
1065		!$omp target exit data map(release: mass, I_mass, depth_tot, src, drag_rate, drag_rate_visc)
1066		!$omp target exit data map(release: drag_vel_u, drag_vel_v, bottomFac2, barotrFac2, LmixScale)
1067		!$omp target exit data map(release: MEKE_uflux, MEKE_vflux, Kh_u, Kh_v, baroHu, baroHv)
1068		!$omp target exit data map(release: damp_rate, damping, MEKE_current, damp_rate_s1)
1069		!$omp target exit data map(release: del2MEKE, del4MEKE, MEKE_decay, equilibrium_value)
1070		!$omp target exit data map(release: src_adv, src_GM, src_mom_K4, src_btm_drag, src_mom_lp, src_mom_bh)
1071
1072		end subroutine step_forward_MEKE
1073
1074		!> Calculates the equilibrium solution where the source depends only on MEKE diffusivity
1075		!! and there is no lateral diffusion of MEKE.
1076		!! Results is in MEKE%MEKE.
1077	1	subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)
1078		type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
1079		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1080		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1081		type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1082		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
1083		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1084		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1085		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: drag_rate_visc !< Mean flow velocity contribution
1086		!! to the MEKE drag rate [H T-1 ~> m s-1 or kg m-2 s-1]
1087		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_mass !< Inverse of column mass [R-1 Z-1 ~> m2 kg-1].
1088		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1089
1090		! Local variables
1091		real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
1092		real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1093		real :: bottomFac2, barotrFac2 ! Vertical structure factors [nondim]
1094		real :: LmixScale, LRhines, LEady ! Various mixing length scales [L ~> m]
1095		real :: KhCoeff ! A copy of MEKE_KhCoeff from the control structure [nondim]
1096		real :: Kh ! A lateral diffusivity [L2 T-1 ~> m2 s-1]
1097		real :: Ubg2 ! Background (tidal?) velocity squared [L2 T-2 ~> m2 s-2]
1098		real :: cd2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
1099		real :: drag_rate ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
1100		real :: src ! The sum of MEKE sources [L2 T-3 ~> W kg-1]
1101		real :: ldamping ! The MEKE damping rate [T-1 ~> s-1].
1102		real :: EKE, EKEmin, EKEmax, EKEerr ! [L2 T-2 ~> m2 s-2]
1103		real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1]
1104		real :: FatH ! Coriolis parameter at h points, used to compute topographic beta [T-1 ~> s-1]
1105		real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
1106		real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]
1107		integer :: i, j, is, ie, js, je, n1, n2
1108		real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2].
1109		logical :: useSecant, debugIteration
1110
1111	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1112
1113	1	debugIteration = .false.
1114	1	KhCoeff = CS%MEKE_KhCoeff
1115	1	Ubg2 = CS%MEKE_Uscale**2
1116	1	cd2 = CS%cdrag**2
1117	1	tolerance = 1.0e-12US%m_s_to_L_T*2
1118	1	h_neglect = GV%H_subroundoff
1119
1120		!$OMP do
1121	7261	do j=js,je ; do i=is,ie
1122		! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
1123		! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
1124	7200	SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1))
1125
1126	7260	if (CS%MEKE_equilibrium_alt) then
1127	0	MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * depth_tot(i,j))**2 / cd2
1128		else
1129		FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
1130	7200	(G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points
1131
1132		! Since zero-bathymetry cells are masked, this avoids calculations on land
1133	7200	if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
1134	7200	beta_topo_x = 0. ; beta_topo_y = 0.
1135		else
1136		!### Consider different combinations of these estimates of topographic beta.
1137		beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1138		(depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) &
1139		/ max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &
1140		+ (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) &
1141	0	/ max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )
1142		beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1143		(depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) &
1144		/ max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &
1145		(depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) &
1146	0	/ max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )
1147		endif
1148		beta = sqrt(((G%dF_dx(i,j) + beta_topo_x)**2) + &
1149	7200	((G%dF_dy(i,j) + beta_topo_y)**2) )
1150
1151	7200	if (KhCoeffSNI_mass(i,j)>0.) then
1152		! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E
1153	4694	EKEmin = 0. ! Use the trivial root as the left bracket
1154	4694	ResMin = 0. ! Need to detect direction of left residual
1155	4694	EKEmax = 0.01US%m_s_to_L_T*2 ! First guess at right bracket
1156	4694	useSecant = .false. ! Start using a bisection method
1157
1158		! First find right bracket for which resid<0
1159	4694	resid = 1.0US%m_to_L2US%T_to_s**3 ; n1 = 0
1160	9388	do while (resid>0.)
1161	4694	n1 = n1 + 1
1162	4694	EKE = EKEmax
1163		call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), &
1164		MEKE%Rd_dx_h(i,j), SN, EKE, &
1165	4694	bottomFac2, barotrFac2, LmixScale, LRhines, LEady)
1166		! TODO: Should include resolution function in Kh
1167	4694	Kh = (KhCoeff * sqrt(2.barotrFac2EKE) * LmixScale)
1168	4694	src = Kh * (SN * SN)
1169	4694	drag_rate = (GV%H_to_RZ * I_mass(i,j)) * sqrt(drag_rate_visc(i,j)*2 + cd2 ( 2.0bottomFac2EKE + Ubg2 ) )
1170	4694	ldamping = CS%MEKE_damping + drag_rate * bottomFac2
1171	4694	resid = src - ldamping * EKE
1172		! if (debugIteration) then
1173		! write(0,*) n1, 'EKE=',EKE,'resid=',resid
1174		! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin
1175		! write(0,*) 'src=',src,'ldamping=',ldamping
1176		! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2
1177		! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2
1178		! endif
1179	4694	if (resid>0.) then ! EKE is to the left of the root
1180	0	EKEmin = EKE ! so we move the left bracket here
1181	0	EKEmax = 10. * EKE ! and guess again for the right bracket
1182	0	if (resid<ResMin) useSecant = .true.
1183	0	ResMin = resid
1184	0	if (EKEmax > 2.e17US%m_s_to_L_T*2) then
1185	0	if (debugIteration) stop 'Something has gone very wrong'
1186	0	debugIteration = .true.
1187	0	resid = 1. ; n1 = 0
1188	0	EKEmin = 0. ; ResMin = 0.
1189	0	EKEmax = 0.01US%m_s_to_L_T*2
1190	0	useSecant = .false.
1191		endif
1192		endif
1193		enddo ! while(resid>0.) searching for right bracket
1194	4694	ResMax = resid
1195
1196		! Bisect the bracket
1197	4694	n2 = 0 ; EKEerr = EKEmax - EKEmin
1198	137051	do while (EKEerr > tolerance)
1199	132357	n2 = n2 + 1
1200	132357	if (useSecant) then
1201	11501	EKE = EKEmin + (EKEmax - EKEmin) * (ResMin / (ResMin - ResMax))
1202		else
1203	120856	EKE = 0.5 * (EKEmin + EKEmax)
1204		endif
1205	132357	EKEerr = min( EKE-EKEmin, EKEmax-EKE )
1206		! TODO: Should include resolution function in Kh
1207	132357	Kh = (KhCoeff * sqrt(2.barotrFac2EKE) * LmixScale)
1208	132357	src = Kh * (SN * SN)
1209	132357	drag_rate = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)*2 + cd2 ( 2.0bottomFac2EKE + Ubg2 ) )
1210	132357	ldamping = CS%MEKE_damping + drag_rate * bottomFac2
1211	132357	resid = src - ldamping * EKE
1212	132357	if (useSecant .and. resid>ResMin) useSecant = .false.
1213	132357	if (resid>0.) then ! EKE is to the left of the root
1214	20693	EKEmin = EKE ! so we move the left bracket here
1215	20693	if (resid<ResMin) useSecant = .true.
1216	20693	ResMin = resid ! Save this for the secant method
1217	111664	elseif (resid<0.) then ! EKE is to the right of the root
1218	111664	EKEmax = EKE ! so we move the right bracket here
1219	111664	ResMax = resid ! Save this for the secant method
1220		else
1221	0	exit ! resid=0 => EKE is exactly at the root
1222		endif
1223	132357	if (n2>200) stop 'Failing to converge?'
1224		enddo ! while(EKEmax-EKEmin>tolerance)
1225
1226		else
1227	2506	EKE = 0.
1228		endif
1229	7200	MEKE%MEKE(i,j) = EKE
1230		endif
1231		enddo ; enddo
1232
1233	1	end subroutine MEKE_equilibrium
1234
1235
1236		!< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into
1237		!! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value
1238	0	subroutine MEKE_equilibrium_restoring(CS, G, GV, US, SN_u, SN_v, depth_tot, &
1239	0	equilibrium_value)
1240		type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
1241		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1242		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type.
1243		type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1244		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1245		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1246		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1247		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: equilibrium_value
1248		!< Equilbrium value of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2]
1249
1250		! Local variables
1251		real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1252		integer :: i, j, is, ie, js, je ! local indices
1253		real :: cd2 ! The square of the drag coefficient [nondim]
1254
1255	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1256	0	cd2 = CS%cdrag**2
1257	0	equilibrium_value(:,:) = 0.0
1258
1259		!$OMP do
1260	0	do j=js,je ; do i=is,ie
1261		! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
1262		! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
1263	0	SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1))
1264	0	equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * depth_tot(i,j))**2 / cd2
1265		enddo ; enddo
1266
1267	0	if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, equilibrium_value, CS%diag)
1268	0	end subroutine MEKE_equilibrium_restoring
1269
1270		!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
1271		!! functions that are ratios of either bottom or barotropic eddy energy to the
1272		!! column eddy energy, respectively. See \ref section_MEKE_equations.
1273	24	subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, &
1274	24	bottomFac2, barotrFac2, LmixScale)
1275		type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1276		type(MEKE_type), intent(in) :: MEKE !< MEKE field
1277		type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
1278		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1279		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1280		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
1281		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
1282		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
1283		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].
1284		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomFac2 !< gamma_b^2 [nondim]
1285		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrFac2 !< gamma_t^2 [nondim]
1286		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
1287		! Local variables
1288	48	real, dimension(SZI_(G),SZJ_(G)) :: LRhines, LEady ! Possible mixing length scales [L ~> m]
1289		real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
1290		real :: SN ! The local Eady growth rate [T-1 ~> s-1]
1291		real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1]
1292		real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
1293		real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]
1294		integer :: i, j, is, ie, js, je
1295
1296	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1297	24	h_neglect = GV%H_subroundoff
1298
1299		! TODO: THIS RUNS ON THE HOST, BE WARY! AAAAAA (see the note at the call site in
1300		! step_forward_MEKE: the intrinsic pow in MEKE_lengthScales_0d is not bitwise CPU==GPU)
1301		!$OMP do
1302	174264	do j=js,je ; do i=is,ie
1303	172800	if (.not.CS%use_old_lscale) then
1304	172800	if (CS%aEady > 0.) then
1305	172800	SN = 0.25 * ( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)) )
1306		else
1307	0	SN = 0.
1308		endif
1309		FatH = 0.25* ( ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) ) + &
1310	172800	( G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1) ) ) ! Coriolis parameter at h points
1311
1312		! If depth_tot is zero, then a division by zero FPE will be raised. In this
1313		! case, we apply Adcroft's rule of reciprocals and set the term to zero.
1314		! Since zero-bathymetry cells are masked, this should not affect values.
1315	172800	if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
1316	172800	beta_topo_x = 0. ; beta_topo_y = 0.
1317		else
1318		!### Consider different combinations of these estimates of topographic beta.
1319		beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1320		(depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) &
1321		/ max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &
1322		+ (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) &
1323	0	/ max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )
1324		beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
1325		(depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) &
1326		/ max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &
1327		(depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) &
1328	0	/ max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )
1329		endif
1330		beta = sqrt(((G%dF_dx(i,j) + beta_topo_x)**2) + &
1331	172800	((G%dF_dy(i,j) + beta_topo_y)**2) )
1332
1333		else
1334	0	beta = 0.
1335		endif
1336		! Returns bottomFac2, barotrFac2 and LmixScale
1337		call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), &
1338		MEKE%Rd_dx_h(i,j), SN, MEKE%MEKE(i,j), &
1339		bottomFac2(i,j), barotrFac2(i,j), LmixScale(i,j), &
1340	174240	LRhines(i,j), LEady(i,j))
1341		enddo ; enddo
1342	24	if (CS%id_Lrhines>0) call post_data(CS%id_LRhines, LRhines, CS%diag)
1343	24	if (CS%id_Leady>0) call post_data(CS%id_LEady, LEady, CS%diag)
1344
1345	24	end subroutine MEKE_lengthScales
1346
1347		!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
1348		!! functions that are ratios of either bottom or barotropic eddy energy to the
1349		!! column eddy energy, respectively. See \ref section_MEKE_equations.
1350	177494	pure subroutine MEKE_lengthScales_0d(CS, US, area, beta, depth_tot, Rd_dx, SN, EKE, &
1351		bottomFac2, barotrFac2, LmixScale, Lrhines, Leady)
1352		type(MEKE_CS), intent(in) :: CS !< MEKE control structure.
1353		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1354		real, intent(in) :: area !< Grid cell area [L2 ~> m2]
1355		real, intent(in) :: beta !< Planetary beta = \f$ \nabla f\f$ [T-1 L-1 ~> s-1 m-1]
1356		real, intent(in) :: depth_tot !< The total thickness of the water column [H ~> m or kg m-2]
1357		real, intent(in) :: Rd_dx !< Resolution Ld/dx [nondim].
1358		real, intent(in) :: SN !< Eady growth rate [T-1 ~> s-1].
1359		real, intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
1360		real, intent(out) :: bottomFac2 !< gamma_b^2 [nondim]
1361		real, intent(out) :: barotrFac2 !< gamma_t^2 [nondim]
1362		real, intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
1363		real, intent(out) :: Lrhines !< Rhines length scale [L ~> m].
1364		real, intent(out) :: Leady !< Eady length scale [L ~> m].
1365		! Local variables
1366		real :: Lgrid, Ldeform, Lfrict ! Length scales [L ~> m]
1367		real :: Ue ! An eddy velocity [L T-1 ~> m s-1]
1368
1369		! Length scale for MEKE derived diffusivity
1370	177494	Lgrid = sqrt(area) ! Grid scale
1371	177494	Ldeform = Lgrid * Rd_dx ! Deformation scale
1372	177494	Lfrict = depth_tot / CS%cdrag ! Frictional arrest scale
1373		! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy
1374		! used in calculating bottom drag
1375	177494	bottomFac2 = CS%MEKE_CD_SCALE**2
1376	177494	if (LfrictCS%MEKE_Cb>0.) bottomFac2 = bottomFac2 + 1./( 1. + CS%MEKE_Cb(Ldeform/Lfrict) )**0.8
1377	177494	bottomFac2 = max(bottomFac2, CS%MEKE_min_gamma)
1378		! gamma_t^2 is the ratio of barotropic eddy energy to mean column eddy energy
1379		! used in the velocity scale for diffusivity
1380	177494	barotrFac2 = 1.
1381	177494	if (LfrictCS%MEKE_Ct>0.) barotrFac2 = 1. / ( 1. + CS%MEKE_Ct(Ldeform/Lfrict) )**0.25
1382	177494	barotrFac2 = max(barotrFac2, CS%MEKE_min_gamma)
1383	177494	if (CS%use_old_lscale) then
1384	0	if (CS%Rd_as_max_scale) then
1385	0	LmixScale = min(Ldeform, Lgrid) ! The smaller of Ld or dx
1386		else
1387	0	LmixScale = Lgrid
1388		endif
1389		else
1390	177494	Ue = sqrt( 2.0 * max( 0., barotrFac2*EKE ) ) ! Barotropic eddy flow scale
1391	177494	Lrhines = sqrt( Ue / max( beta, 1.e-30US%T_to_sUS%L_to_m ) ) ! Rhines scale
1392	177494	if (CS%aEady > 0.) then
1393	177494	Leady = Ue / max( SN, 1.e-15*US%T_to_s ) ! Bound Eady time-scale < 1e15 seconds
1394		else
1395	0	Leady = 0.
1396		endif
1397	177494	if (CS%use_min_lscale) then
1398	0	LmixScale = CS%lscale_maxval
1399	0	if (CS%aDeformLdeform > 0.) LmixScale = min(LmixScale,CS%aDeformLdeform)
1400	0	if (CS%aFrict Lfrict > 0.) LmixScale = min(LmixScale,CS%aFrict Lfrict)
1401	0	if (CS%aRhinesLrhines > 0.) LmixScale = min(LmixScale,CS%aRhinesLrhines)
1402	0	if (CS%aEady Leady > 0.) LmixScale = min(LmixScale,CS%aEady Leady)
1403	0	if (CS%aGrid Lgrid > 0.) LmixScale = min(LmixScale,CS%aGrid Lgrid)
1404	0	if (CS%Lfixed > 0.) LmixScale = min(LmixScale,CS%Lfixed)
1405		else
1406	177494	LmixScale = 0.
1407	177494	if (CS%aDeformLdeform > 0.) LmixScale = LmixScale + 1./(CS%aDeformLdeform)
1408	177494	if (CS%aFrict Lfrict > 0.) LmixScale = LmixScale + 1./(CS%aFrict Lfrict)
1409	177494	if (CS%aRhinesLrhines > 0.) LmixScale = LmixScale + 1./(CS%aRhinesLrhines)
1410	177494	if (CS%aEady Leady > 0.) LmixScale = LmixScale + 1./(CS%aEady Leady)
1411	177494	if (CS%aGrid Lgrid > 0.) LmixScale = LmixScale + 1./(CS%aGrid Lgrid)
1412	177494	if (CS%Lfixed > 0.) LmixScale = LmixScale + 1./CS%Lfixed
1413	177494	if (LmixScale > 0.) LmixScale = 1. / LmixScale
1414		endif
1415		endif
1416
1417	177494	end subroutine MEKE_lengthScales_0d
1418
1419		!> Initializes the MOM_MEKE module and reads parameters.
1420		!! Returns True if module is to be used, otherwise returns False.
1421	1	logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, MEKE, restart_CS, meke_in_dynamics)
1422		type(time_type), intent(in) :: Time !< The current model time.
1423		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1424		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure.
1425		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1426		type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1427		type(dbcomms_CS_type), intent(in) :: dbcomms_CS !< Database communications control structure
1428		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
1429		type(MEKE_CS), intent(inout) :: CS !< MEKE control structure.
1430		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
1431		type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure
1432		logical, intent( out) :: meke_in_dynamics !< If true, MEKE is stepped forward in dynamics
1433		!! otherwise in tracer dynamics
1434
1435		! Local variables
1436		real :: MEKE_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value [T ~> s]
1437		real :: cdrag ! The default bottom drag coefficient [nondim].
1438		character(len=200) :: eke_filename, eke_varname, inputdir
1439		character(len=16) :: eke_source_str
1440		integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1441		logical :: laplacian, biharmonic, coldStart
1442		! This include declares and sets the variable "version".
1443		# include "version_variable.h"
1444		character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
1445
1446	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1447	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1448
1449		! Determine whether this module will be used
1450	1	call get_param(param_file, mdl, "USE_MEKE", MEKE_init, default=.false., do_not_log=.true.)
1451	1	call log_version(param_file, mdl, version, "", all_default=.not.MEKE_init)
1452		call get_param(param_file, mdl, "USE_MEKE", MEKE_init, &
1453		"If true, turns on the MEKE scheme which calculates "// &
1454		"a sub-grid mesoscale eddy kinetic energy budget.", &
1455	1	default=.false.)
1456	1	if (.not. MEKE_init) return
1457	1	CS%initialized = .true.
1458		call get_param(param_file, mdl, "MEKE_IN_DYNAMICS", meke_in_dynamics, &
1459		"If true, step MEKE forward with the dynamics "// &
1460		"otherwise with the tracer timestep.", &
1461	1	default=.true.)
1462
1463		call get_param(param_file, mdl, "EKE_SOURCE", eke_source_str, &
1464		"Determine the where EKE comes from:\n" // &
1465		" 'prog': Calculated solving EKE equation\n"// &
1466		" 'file': Read in from a file\n" // &
1467	1	" 'dbclient': Retrieved from ML-database", default='prog')
1468
1469	1	call MOM_mesg("MEKE_init: reading parameters ", 5)
1470
1471	1	select case (lowercase(eke_source_str))
1472		case("file")
1473	0	CS%eke_src = EKE_FILE
1474	0	call time_interp_external_init
1475		call get_param(param_file, mdl, "EKE_FILE", eke_filename, &
1476		"A file in which to find the eddy kineteic energy variable.", &
1477	0	default="eke_file.nc")
1478		call get_param(param_file, mdl, "EKE_VARIABLE", eke_varname, &
1479		"The name of the eddy kinetic energy variable to read from "//&
1480		"EKE_FILE to use in MEKE.", &
1481	0	default="eke")
1482		call get_param(param_file, mdl, "INPUTDIR", inputdir, &
1483		"The directory in which all input files are found.", &
1484	0	default=".", do_not_log=.true.)
1485	0	inputdir = slasher(inputdir)
1486
1487	0	eke_filename = trim(inputdir) // trim(eke_filename)
1488	0	CS%eke_handle = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain)
1489		case("prog")
1490	1	CS%eke_src = EKE_PROG
1491		! Read all relevant parameters and write them to the model log.
1492		call get_param(param_file, mdl, "MEKE_DAMPING", CS%MEKE_damping, &
1493		"The local depth-independent MEKE dissipation rate.", &
1494	1	units="s-1", default=0.0, scale=US%T_to_s)
1495		call get_param(param_file, mdl, "MEKE_CD_SCALE", CS%MEKE_Cd_scale, &
1496		"The ratio of the bottom eddy velocity to the column mean "//&
1497		"eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1 "//&
1498		"to account for the surface intensification of MEKE.", &
1499	1	units="nondim", default=0.)
1500		call get_param(param_file, mdl, "MEKE_CB", CS%MEKE_Cb, &
1501		"A coefficient in the expression for the ratio of bottom projected "//&
1502		"eddy energy and mean column energy (see Jansen et al. 2015).",&
1503	1	units="nondim", default=25.)
1504		call get_param(param_file, mdl, "MEKE_MIN_GAMMA2", CS%MEKE_min_gamma, &
1505		"The minimum allowed value of gamma_b^2.",&
1506	1	units="nondim", default=0.0001)
1507		call get_param(param_file, mdl, "MEKE_CT", CS%MEKE_Ct, &
1508		"A coefficient in the expression for the ratio of barotropic "//&
1509		"eddy energy and mean column energy (see Jansen et al. 2015).",&
1510	1	units="nondim", default=50.)
1511		call get_param(param_file, mdl, "MEKE_GMCOEFF", CS%MEKE_GMcoeff, &
1512		"The efficiency of the conversion of potential energy "//&
1513		"into MEKE by the thickness mixing parameterization. "//&
1514		"If MEKE_GMCOEFF is negative, this conversion is not "//&
1515	1	"used or calculated.", units="nondim", default=-1.0)
1516		call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, &
1517		"If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//&
1518	1	"from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
1519		call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, &
1520		"The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
1521	1	"thickness diffusion.", units="nondim", default=0.05)
1522		call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", CS%MEKE_equilibrium_alt, &
1523		"If true, use an alternative formula for computing the (equilibrium) "//&
1524	1	"initial value of MEKE.", default=.false.)
1525		call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", CS%MEKE_equilibrium_restoring, &
1526		"If true, restore MEKE back to its equilibrium value, which is calculated at "//&
1527	1	"each time step.", default=.false.)
1528	1	if (CS%MEKE_equilibrium_restoring) then
1529		call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, &
1530		"The timescale used to nudge MEKE toward its equilibrium value.", &
1531	0	units="s", default=1e6, scale=US%s_to_T)
1532	0	CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale
1533		endif
1534
1535		call get_param(param_file, mdl, "MEKE_FRCOEFF", CS%MEKE_FrCoeff, &
1536		"The efficiency of the conversion of mean energy into "//&
1537		"MEKE. If MEKE_FRCOEFF is negative, this conversion "//&
1538	1	"is not used or calculated.", units="nondim", default=-1.0)
1539		call get_param(param_file, mdl, "MEKE_BHFRCOEFF", CS%MEKE_bhFrCoeff, &
1540		"The efficiency of the conversion of mean energy into "//&
1541		"MEKE by the biharmonic dissipation. If MEKE_bhFRCOEFF is negative, this conversion "//&
1542	1	"is not used or calculated.", units="nondim", default=-1.0)
1543		call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, &
1544		"The efficiency of the conversion of MEKE into mean energy "//&
1545		"by GME. If MEKE_GMECOEFF is negative, this conversion "//&
1546	1	"is not used or calculated.", units="nondim", default=-1.0)
1547		call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, &
1548		"A background energy source for MEKE.", &
1549	1	units="W kg-1", default=0.0, scale=US%m_to_L*2US%T_to_s**3)
1550		call get_param(param_file, mdl, "MEKE_KH", CS%MEKE_Kh, &
1551		"A background lateral diffusivity of MEKE. "//&
1552		"Use a negative value to not apply lateral diffusion to MEKE.", &
1553	1	units="m2 s-1", default=-1.0, scale=US%m_to_L*2US%T_to_s)
1554		call get_param(param_file, mdl, "MEKE_K4", CS%MEKE_K4, &
1555		"A lateral bi-harmonic diffusivity of MEKE. "//&
1556		"Use a negative value to not apply bi-harmonic diffusion to MEKE.", &
1557	1	units="m4 s-1", default=-1.0, scale=US%m_to_L*4US%T_to_s)
1558		call get_param(param_file, mdl, "MEKE_DTSCALE", CS%MEKE_dtScale, &
1559		"A scaling factor to accelerate the time evolution of MEKE.", &
1560	1	units="nondim", default=1.0)
1561		call get_param(param_file, mdl, "MEKE_POSITIVE", CS%MEKE_positive, &
1562		"If true, it guarantees that MEKE will always be >= 0.", &
1563	1	default=.false.)
1564		case("dbclient")
1565	0	CS%eke_src = EKE_DBCLIENT
1566	0	call ML_MEKE_init(diag, G, US, Time, param_file, dbcomms_CS, CS)
1567		case default
1568	1	call MOM_error(FATAL, "Invalid method selected for calculating EKE")
1569		end select
1570		! GMM, make sure all parameters used to calculated MEKE are within the above if
1571
1572		call get_param(param_file, mdl, "MEKE_KHCOEFF", CS%MEKE_KhCoeff, &
1573		"A scaling factor in the expression for eddy diffusivity "//&
1574		"which is otherwise proportional to the MEKE velocity- "//&
1575		"scale times an eddy mixing-length. This factor "//&
1576		"must be >0 for MEKE to contribute to the thickness/ "//&
1577		"and tracer diffusivity in the rest of the model.", &
1578	1	units="nondim", default=1.0)
1579		call get_param(param_file, mdl, "MEKE_USCALE", CS%MEKE_Uscale, &
1580		"The background velocity that is combined with MEKE to "//&
1581	1	"calculate the bottom drag.", units="m s-1", default=0.0, scale=US%m_s_to_L_T)
1582		call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
1583		"If true, use the GM energy conversion form S^2N^2kappa rather "//&
1584	1	"than the streamfunction for the MEKE GM source term.", default=.false.)
1585		call get_param(param_file, mdl, "MEKE_MIN_DEPTH_TOT", CS%MEKE_min_depth_tot, &
1586		"The minimum total depth over which to distribute MEKE energy sources. "//&
1587		"When the total depth is less than this, the sources are scaled away.", &
1588	1	units="m", default=1.0, scale=GV%m_to_H, do_not_log=.not.CS%GM_src_alt)
1589		call get_param(param_file, mdl, "MEKE_VISC_DRAG", CS%visc_drag, &
1590		"If true, use the vertvisc_type to calculate the bottom "//&
1591	1	"drag acting on MEKE.", default=.true.)
1592		call get_param(param_file, mdl, "MEKE_KHTH_FAC", MEKE%KhTh_fac, &
1593	1	"A factor that maps MEKE%Kh to KhTh.", units="nondim", default=0.0)
1594		call get_param(param_file, mdl, "MEKE_KHTR_FAC", MEKE%KhTr_fac, &
1595	1	"A factor that maps MEKE%Kh to KhTr.", units="nondim", default=0.0)
1596		call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", CS%KhMEKE_Fac, &
1597		"A factor that maps MEKE%Kh to Kh for MEKE itself.", &
1598	1	units="nondim", default=0.0)
1599		call get_param(param_file, mdl, "MEKE_OLD_LSCALE", CS%use_old_lscale, &
1600		"If true, use the old formula for length scale which is "//&
1601		"a function of grid spacing and deformation radius.", &
1602	1	default=.false.)
1603		call get_param(param_file, mdl, "MEKE_MIN_LSCALE", CS%use_min_lscale, &
1604		"If true, use a strict minimum of provided length scales "//&
1605		"rather than harmonic mean.", &
1606	1	default=.false.)
1607		call get_param(param_file, mdl, "MEKE_LSCALE_MAX_VAL", CS%lscale_maxval, &
1608		"The ceiling on the value of the MEKE length scale when MEKE_MIN_LSCALE=True. "//&
1609		"The default is the distance from the equator to the pole on Earth, as "//&
1610		"estimated by enlightenment era scientists, but should probably scale with RAD_EARTH.", &
1611	1	units="m", default=1.0e7, scale=US%m_to_L, do_not_log=.not.CS%use_min_lscale)
1612		call get_param(param_file, mdl, "MEKE_RD_MAX_SCALE", CS%Rd_as_max_scale, &
1613		"If true, the length scale used by MEKE is the minimum of "//&
1614		"the deformation radius or grid-spacing. Only used if "//&
1615	1	"MEKE_OLD_LSCALE=True", default=.false.)
1616		call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", CS%viscosity_coeff_Ku, &
1617		"If non-zero, is the scaling coefficient in the expression for "//&
1618		"viscosity used to parameterize harmonic lateral momentum mixing by "//&
1619		"unresolved eddies represented by MEKE. Can be negative to "//&
1620		"represent backscatter from the unresolved eddies.", &
1621	1	units="nondim", default=0.0)
1622		call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", CS%viscosity_coeff_Au, &
1623		"If non-zero, is the scaling coefficient in the expression for "//&
1624		"viscosity used to parameterize biharmonic lateral momentum mixing by "//&
1625		"unresolved eddies represented by MEKE. Can be negative to "//&
1626		"represent backscatter from the unresolved eddies.", &
1627	1	units="nondim", default=0.0)
1628		call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", CS%Lfixed, &
1629		"If positive, is a fixed length contribution to the expression "//&
1630		"for mixing length used in MEKE-derived diffusivity.", &
1631	1	units="m", default=0.0, scale=US%m_to_L)
1632		call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", CS%fixed_total_depth, &
1633		"If true, use the nominal bathymetric depth as the estimate of the "//&
1634		"time-varying ocean depth. Otherwise base the depth on the total ocean mass "//&
1635	1	"per unit area.", default=.true.)
1636		call get_param(param_file, mdl, "MEKE_TOTAL_DEPTH_RHO", CS%rho_fixed_total_depth, &
1637		"A density used to translate the nominal bathymetric depth into an estimate "//&
1638		"of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH is true.", &
1639		units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, &
1640	1	do_not_log=(GV%Boussinesq.or.(.not.CS%fixed_total_depth)))
1641
1642		call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", CS%aDeform, &
1643		"If positive, is a coefficient weighting the deformation scale "//&
1644		"in the expression for mixing length used in MEKE-derived diffusivity.", &
1645	1	units="nondim", default=0.0)
1646		call get_param(param_file, mdl, "MEKE_ALPHA_RHINES", CS%aRhines, &
1647		"If positive, is a coefficient weighting the Rhines scale "//&
1648		"in the expression for mixing length used in MEKE-derived diffusivity.", &
1649	1	units="nondim", default=0.0)
1650		call get_param(param_file, mdl, "MEKE_ALPHA_EADY", CS%aEady, &
1651		"If positive, is a coefficient weighting the Eady length scale "//&
1652		"in the expression for mixing length used in MEKE-derived diffusivity.", &
1653	1	units="nondim", default=0.0)
1654		call get_param(param_file, mdl, "MEKE_ALPHA_FRICT", CS%aFrict, &
1655		"If positive, is a coefficient weighting the frictional arrest scale "//&
1656		"in the expression for mixing length used in MEKE-derived diffusivity.", &
1657	1	units="nondim", default=0.0)
1658		call get_param(param_file, mdl, "MEKE_ALPHA_GRID", CS%aGrid, &
1659		"If positive, is a coefficient weighting the grid-spacing as a scale "//&
1660		"in the expression for mixing length used in MEKE-derived diffusivity.", &
1661	1	units="nondim", default=0.0)
1662		call get_param(param_file, mdl, "MEKE_COLD_START", coldStart, &
1663		"If true, initialize EKE to zero. Otherwise a local equilibrium solution "//&
1664	1	"is used as an initial condition for EKE.", default=.false.)
1665		call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", MEKE%backscatter_Ro_c, &
1666		"The coefficient in the Rossby number function for scaling the biharmonic "//&
1667		"frictional energy source. Setting to non-zero enables the Rossby number function.", &
1668	1	units="nondim", default=0.0)
1669		call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", MEKE%backscatter_Ro_pow, &
1670		"The power in the Rossby number function for scaling the biharmonic "//&
1671	1	"frictional energy source.", units="nondim", default=0.0)
1672		call get_param(param_file, mdl, "MEKE_ADVECTION_FACTOR", CS%MEKE_advection_factor, &
1673		"A scale factor in front of advection of eddy energy. Zero turns advection off. "//&
1674		"Using unity would be normal but other values could accommodate a mismatch "//&
1675		"between the advecting barotropic flow and the vertical structure of MEKE.", &
1676	1	units="nondim", default=0.0)
1677		call get_param(param_file, mdl, "MEKE_ADVECTION_BUG", CS%MEKE_advection_bug, &
1678		"If true, recover a bug in the calculation of the barotropic transport for "//&
1679		"the advection of MEKE. With the bug, only the transports in the deepest "//&
1680	1	"layer are used.", default=.false., do_not_log=(CS%MEKE_advection_factor<=0.))
1681		call get_param(param_file, mdl, "MEKE_TOPOGRAPHIC_BETA", CS%MEKE_topographic_beta, &
1682		"A scale factor to determine how much topographic beta is weighed in " //&
1683		"computing beta in the expression of Rhines scale. Use 1 if full "//&
1684		"topographic beta effect is considered; use 0 if it's completely ignored.", &
1685	1	units="nondim", default=0.0)
1686		call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, &
1687		"If true, the eddy scale of MEKE is used for the SQG vertical structure ",&
1688	1	default=.false.)
1689
1690		! Nonlocal module parameters
1691		call get_param(param_file, mdl, "CDRAG", cdrag, &
1692		"CDRAG is the drag coefficient relating the magnitude of the velocity "//&
1693	1	"field to the bottom stress.", units="nondim", default=0.003)
1694		call get_param(param_file, mdl, "MEKE_CDRAG", CS%cdrag, &
1695		"Drag coefficient relating the magnitude of the velocity "//&
1696	1	"field to the bottom stress in MEKE.", units="nondim", default=cdrag, scale=US%L_to_m*GV%m_to_H)
1697	1	call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
1698	1	call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)
1699
1700	1	if (CS%viscosity_coeff_Ku/=0. .and. .not. laplacian) call MOM_error(FATAL, &
1701	0	"LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.")
1702
1703	1	if (CS%viscosity_coeff_Au/=0. .and. .not. biharmonic) call MOM_error(FATAL, &
1704	0	"BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.")
1705
1706	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
1707
1708		! Identify if any lateral diffusive processes are active
1709	1	CS%kh_flux_enabled = .false.
1710	1	if ((CS%MEKE_KH >= 0.0) .or. (CS%KhMEKE_FAC > 0.0) .or. (CS%MEKE_advection_factor > 0.0)) &
1711	1	CS%kh_flux_enabled = .true.
1712
1713		! Register fields for output from this module.
1714	1	CS%diag => diag
1715		CS%id_MEKE = register_diag_field('ocean_model', 'MEKE', diag%axesT1, Time, &
1716	1	'Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2)
1717	1	if (.not. allocated(MEKE%MEKE)) CS%id_MEKE = -1
1718		CS%id_Kh = register_diag_field('ocean_model', 'MEKE_KH', diag%axesT1, Time, &
1719	1	'MEKE derived diffusivity', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1720	1	if (.not. allocated(MEKE%Kh)) CS%id_Kh = -1
1721		CS%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, Time, &
1722	1	'MEKE derived lateral viscosity', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1723	1	if (.not. allocated(MEKE%Ku)) CS%id_Ku = -1
1724		CS%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, Time, &
1725	1	'MEKE derived lateral biharmonic viscosity', 'm4 s-1', conversion=US%L_to_m*4US%s_to_T)
1726	1	if (.not. allocated(MEKE%Au)) CS%id_Au = -1
1727		CS%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, Time, &
1728	1	'MEKE derived eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
1729	1	if (.not. allocated(MEKE%MEKE)) CS%id_Ue = -1
1730		CS%id_Ub = register_diag_field('ocean_model', 'MEKE_Ub', diag%axesT1, Time, &
1731	1	'MEKE derived bottom eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
1732	1	if (.not. allocated(MEKE%MEKE)) CS%id_Ub = -1
1733		CS%id_Ut = register_diag_field('ocean_model', 'MEKE_Ut', diag%axesT1, Time, &
1734	1	'MEKE derived barotropic eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s)
1735	1	if (.not. allocated(MEKE%MEKE)) CS%id_Ut = -1
1736		CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, &
1737	1	'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1738
1739		CS%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, Time, &
1740	1	'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1741
1742		CS%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, Time, &
1743	1	'MEKE energy source from the bottom drag acting on MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1744
1745	1	if (CS%MEKE_K4 >= 0.) &
1746		CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', &
1747		diag%axesT1, Time, 'MEKE energy source from the biharmonic of MEKE', &
1748	0	'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1749
1750	1	if (CS%MEKE_GMcoeff >= 0.) &
1751		CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', &
1752		diag%axesT1, Time, 'MEKE energy source from the thickness mixing (GM scheme)', &
1753	1	'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1754
1755	1	if (CS%MEKE_FrCoeff >= 0.) &
1756		CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', &
1757		diag%axesT1, Time, 'MEKE energy source from the Laplacian of resolved flows', &
1758	0	'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1759
1760	1	if (CS%MEKE_bhFrCoeff >= 0.) &
1761		CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', &
1762		diag%axesT1, Time, 'MEKE energy source from the biharmonic of resolved flows', &
1763	0	'm2 s-3', conversion=(US%L_T_to_m_s*2)US%s_to_T)
1764
1765		CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, &
1766	1	'MEKE decay rate', 's-1', conversion=US%s_to_T)
1767		CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, &
1768		'MEKE energy available from thickness mixing', &
1769	1	'W m-2', conversion=US%RZ3_T3_to_W_m2US%L_to_Z*2)
1770	1	if (.not. allocated(MEKE%GM_src)) CS%id_GM_src = -1
1771		CS%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, Time, &
1772		'MEKE energy available from momentum', &
1773	1	'W m-2', conversion=US%RZ3_T3_to_W_m2US%L_to_Z*2)
1774	1	if (.not. allocated(MEKE%mom_src)) CS%id_mom_src = -1
1775		CS%id_mom_src_bh = register_diag_field('ocean_model', 'MEKE_mom_src_bh',diag%axesT1, Time, &
1776		'MEKE energy available from the biharmonic dissipation of momentum', &
1777	1	'W m-2', conversion=US%RZ3_T3_to_W_m2US%L_to_Z*2)
1778	1	if (.not. allocated(MEKE%mom_src_bh)) CS%id_mom_src_bh = -1
1779		CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, &
1780		'MEKE energy lost to GME backscatter', &
1781	1	'W m-2', conversion=US%RZ3_T3_to_W_m2US%L_to_Z*2)
1782	1	if (.not. allocated(MEKE%GME_snk)) CS%id_GME_snk = -1
1783		CS%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, Time, &
1784	1	'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1785		CS%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, Time, &
1786	1	'Rhines length scale used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1787		CS%id_Leady = register_diag_field('ocean_model', 'MEKE_Leady', diag%axesT1, Time, &
1788	1	'Eady length scale used in the MEKE derived eddy diffusivity', 'm', conversion=US%L_to_m)
1789		CS%id_gamma_b = register_diag_field('ocean_model', 'MEKE_gamma_b', diag%axesT1, Time, &
1790	1	'Ratio of bottom-projected eddy velocity to column-mean eddy velocity', 'nondim')
1791		CS%id_gamma_t = register_diag_field('ocean_model', 'MEKE_gamma_t', diag%axesT1, Time, &
1792	1	'Ratio of barotropic eddy velocity to column-mean eddy velocity', 'nondim')
1793
1794	1	if (CS%kh_flux_enabled) then
1795		CS%id_KhMEKE_u = register_diag_field('ocean_model', 'KHMEKE_u', diag%axesCu1, Time, &
1796	1	'Zonal diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1797		CS%id_KhMEKE_v = register_diag_field('ocean_model', 'KHMEKE_v', diag%axesCv1, Time, &
1798	1	'Meridional diffusivity of MEKE', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1799		endif
1800
1801	1	if (CS%MEKE_equilibrium_restoring) then
1802		CS%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, Time, &
1803	0	'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=US%L_T_to_m_s**2)
1804		endif
1805
1806	1	CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE)
1807
1808
1809		! Detect whether this instance of MEKE_init() is at the beginning of a run
1810		! or after a restart. If at the beginning, we will initialize MEKE to a local
1811		! equilibrium.
1812	1	CS%initialize = .not.query_initialized(MEKE%MEKE, "MEKE", restart_CS)
1813	1	if (coldStart) CS%initialize = .false.
1814	1	if (CS%initialize) call MOM_error(WARNING, &
1815	1	"MEKE_init: Initializing MEKE with a local equilibrium balance.")
1816	1	if (allocated(MEKE%Le)) then
1817	0	if (.not.query_initialized(MEKE%Le, "MEKE_Le", restart_CS)) then
1818		!$OMP parallel do default(shared)
1819	0	do j=js,je ; do i=is,ie
1820	0	MEKE%Le(i,j) = sqrt(G%areaT(i,j))
1821		enddo ; enddo
1822		endif
1823		endif
1824
1825		! Set up group passes. In the case of a restart, these fields need a halo update now.
1826	1	if (allocated(MEKE%MEKE)) then
1827	1	call create_group_pass(CS%pass_MEKE, MEKE%MEKE, G%Domain)
1828	1	if (allocated(MEKE%Kh_diff)) call create_group_pass(CS%pass_MEKE, MEKE%Kh_diff, G%Domain)
1829	1	if (.not.CS%initialize) call do_group_pass(CS%pass_MEKE, G%Domain)
1830		endif
1831	1	if (allocated(MEKE%Kh)) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain)
1832	1	if (allocated(MEKE%Ku)) call create_group_pass(CS%pass_Kh, MEKE%Ku, G%Domain)
1833	1	if (allocated(MEKE%Au)) call create_group_pass(CS%pass_Kh, MEKE%Au, G%Domain)
1834	1	if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain)
1835
1836		if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) &
1837	1	.or. allocated(MEKE%Le)) &
1838	1	call do_group_pass(CS%pass_Kh, G%Domain)
1839
1840		! This CS and the MEKE type are inline members of MOM_control_struct, whose whole address
1841		! range the solo driver puts on the device with `enter data map(alloc: MOM_CSp)` -- allocated,
1842		! never copied. Any later map(to:) of them is a present-table no-op, so device kernels would
1843		! read zeros for CS scalars. `target update` always copies; issue it here, once the
1844		! parameters are final, and BEFORE the component maps below, whose attach then fixes the
1845		! MEKE arrays' descriptors in the device copy.
1846		!$omp target update to(CS)
1847		if (allocated(MEKE%MEKE)) then
1848		!$omp target enter data map(to: MEKE%MEKE)
1849		endif
1850		if (allocated(MEKE%GM_src)) then
1851		!$omp target enter data map(to: MEKE%GM_src)
1852		endif
1853		if (allocated(MEKE%mom_src)) then
1854		!$omp target enter data map(to: MEKE%mom_src)
1855		endif
1856		if (allocated(MEKE%mom_src_bh)) then
1857		!$omp target enter data map(to: MEKE%mom_src_bh)
1858		endif
1859		if (allocated(MEKE%GME_snk)) then
1860		!$omp target enter data map(to: MEKE%GME_snk)
1861		endif
1862		if (allocated(MEKE%Kh)) then
1863		!$omp target enter data map(to: MEKE%Kh)
1864		endif
1865		if (allocated(MEKE%Kh_diff)) then
1866		!$omp target enter data map(to: MEKE%Kh_diff)
1867		endif
1868		if (allocated(MEKE%Rd_dx_h)) then
1869		!$omp target enter data map(to: MEKE%Rd_dx_h)
1870		endif
1871		if (allocated(MEKE%Ku)) then
1872		!$omp target enter data map(to: MEKE%Ku)
1873		endif
1874		if (allocated(MEKE%Au)) then
1875		!$omp target enter data map(to: MEKE%Au)
1876		endif
1877		if (allocated(MEKE%Le)) then
1878		!$omp target enter data map(to: MEKE%Le)
1879		endif
1880
1881	1	end function MEKE_init
1882
1883		!> Initializer for the variant of MEKE that uses ML to predict eddy kinetic energy
1884	0	subroutine ML_MEKE_init(diag, G, US, Time, param_file, dbcomms_CS, CS)
1885		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
1886		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1887		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1888		type(time_type), intent(in) :: Time !< The current model time.
1889		type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1890		type(dbcomms_CS_type), intent(in) :: dbcomms_CS !< Control structure for database communication
1891		type(MEKE_CS), intent(inout) :: CS !< Control structure for this module
1892
1893		character(len=200) :: inputdir, backend, model_filename
1894		integer :: db_return_code, batch_size
1895		character(len=40) :: mdl = "MOM_ML_MEKE"
1896		#ifdef __NVCOMPILER_OPENMP_GPU
1897		integer, parameter :: default_niblock = 0
1898		integer, parameter :: default_njblock = 0
1899		integer, parameter :: default_nkblock = 0
1900		#else
1901		integer, parameter :: default_niblock = 0
1902		integer, parameter :: default_njblock = 1
1903		integer, parameter :: default_nkblock = 1
1904		#endif
1905
1906		! Store pointers in control structure
1907	0	write(CS%key_suffix, '(A,I6.6)') '_', PE_here()
1908		! Put some basic information into the database
1909	0	db_return_code = 0
1910		db_return_code = CS%client%put_tensor("meta"//CS%key_suffix, &
1911	0	REAL([G%isd_global, G%idg_offset, G%jsd_global, G%jdg_offset]),[4]) + db_return_code
1912	0	db_return_code = CS%client%put_tensor("geolat"//CS%key_suffix, G%geoLatT, shape(G%geoLatT)) + db_return_code
1913	0	db_return_code = CS%client%put_tensor("geolon"//CS%key_suffix, G%geoLonT, shape(G%geoLonT)) + db_return_code
1914	0	db_return_code = CS%client%put_tensor("EKE_shape"//CS%key_suffix, shape(G%geolonT), [2]) + db_return_code
1915
1916	0	if (CS%client%SR_error_parser(db_return_code)) call MOM_error(FATAL, "Putting metadata into the database failed")
1917
1918	0	call read_param(param_file, "INPUTDIR", inputdir)
1919	0	inputdir = slasher(inputdir)
1920
1921	0	call get_param(param_file, mdl, "BATCH_SIZE", batch_size, "Batch size to use for inference", default=1)
1922		call get_param(param_file, mdl, "EKE_BACKEND", backend, &
1923	0	"The computational backend to use for EKE inference (CPU or GPU)", default="GPU")
1924		call get_param(param_file, mdl, "EKE_MODEL", model_filename, &
1925	0	"Filename of the a saved pyTorch model to use", fail_if_missing = .true.)
1926		call get_param(param_file, mdl, "EKE_MAX", CS%eke_max, &
1927		"Maximum value of EKE allowed when inferring EKE", &
1928	0	units="m2 s-2", default=2., scale=US%m_s_to_L_T**2)
1929
1930		! Set the machine learning model
1931	0	if (dbcomms_CS%colocated) then
1932	0	if (modulo(PE_here(),dbcomms_CS%colocated_stride) == 0) then
1933		db_return_code = CS%client%set_model_from_file(CS%model_key, trim(inputdir)//trim(model_filename), &
1934	0	"TORCH", backend, batch_size=batch_size)
1935		endif
1936		else
1937	0	if (is_root_pe()) then
1938		db_return_code = CS%client%set_model_from_file(CS%model_key, trim(inputdir)//trim(model_filename), &
1939	0	"TORCH", backend, batch_size=batch_size)
1940		endif
1941		endif
1942	0	if (CS%client%SR_error_parser(db_return_code)) then
1943	0	call MOM_error(FATAL, "MEKE: set_model failed")
1944		endif
1945
1946		call get_param(param_file, mdl, "ONLINE_ANALYSIS", CS%online_analysis, &
1947	0	"If true, post EKE used in MOM6 to the database for analysis", default=.true.)
1948
1949		! Set various clock ids
1950	0	CS%id_client_init = cpu_clock_id('(ML_MEKE client init)', grain=CLOCK_ROUTINE)
1951	0	CS%id_put_tensor = cpu_clock_id('(ML_MEKE put tensor)', grain=CLOCK_ROUTINE)
1952	0	CS%id_run_model = cpu_clock_id('(ML_MEKE run model)', grain=CLOCK_ROUTINE)
1953	0	CS%id_unpack_tensor = cpu_clock_id('(ML_MEKE unpack tensor )', grain=CLOCK_ROUTINE)
1954
1955		! Diagnostics for ML_MEKE
1956		CS%id_mke = register_diag_field('ocean_model', 'MEKE_MKE', diag%axesT1, Time, &
1957	0	'Surface mean (resolved) kinetic energy used in MEKE', 'm2 s-2', conversion=US%L_T_to_m_s**2)
1958		CS%id_slope_z= register_diag_field('ocean_model', 'MEKE_slope_z', diag%axesT1, Time, &
1959	0	'Vertically averaged isopyncal slope magnitude used in MEKE', 'nondim', conversion=US%Z_to_L)
1960		CS%id_slope_x= register_diag_field('ocean_model', 'MEKE_slope_x', diag%axesCui, Time, &
1961	0	'Isopycnal slope in the x-direction used in MEKE', 'nondim', conversion=US%Z_to_L)
1962		CS%id_slope_y= register_diag_field('ocean_model', 'MEKE_slope_y', diag%axesCvi, Time, &
1963	0	'Isopycnal slope in the y-direction used in MEKE', 'nondim', conversion=US%Z_to_L)
1964		CS%id_rv = register_diag_field('ocean_model', 'MEKE_RV', diag%axesT1, Time, &
1965	0	'Surface relative vorticity used in MEKE', 's-1', conversion=US%s_to_T)
1966
1967		! Isoneutral blocking parameters
1968		call get_param(param_file, mdl, "ISOPYCNAL_NIBLOCK", CS%niblock, &
1969		"The i-direction block size used to calculate isopycnal slopes. "//&
1970		"If 0, or when running with OpenMP offload, "//&
1971		"the full computational domain width is used. "//&
1972		"If USE_STANLEY_ISO is true, ISOPYCNAL_NIBLOCK cannot equal 1.", &
1973	0	default=default_niblock, layoutParam=.true.)
1974		call get_param(param_file, mdl, "ISOPYCNAL_NJBLOCK", CS%njblock, &
1975		"The j-direction block size used to calculate isopycnal slopes. "//&
1976		"If 0, defaults to 1, except when running with OpenMP offload, "//&
1977		"in which case the full computational domain height is used. " //&
1978		"If USE_STANLEY_ISO is true, ISOPYCNAL_NJBLOCK cannot equal 1.", &
1979	0	default=default_njblock, layoutParam=.true.)
1980		call get_param(param_file, mdl, "ISOPYCNAL_NKBLOCK", CS%nkblock, &
1981		"The k-direction block size used to calculate isopycnal slopes. "//&
1982		"If 0, defaults to 1, except when "//&
1983		"running with OpenMP offload, in which case the full computational "//&
1984	0	"domain depth is used.", default=default_nkblock, layoutParam=.true.)
1985
1986	0	if (CS%niblock < 0) &
1987		call MOM_error(FATAL, "ISOPYCNAL_NIBLOCK must be nonnegative; "//&
1988	0	"use 0 to select the default block size.")
1989	0	if (CS%njblock < 0) &
1990		call MOM_error(FATAL, "ISOPYCNAL_NJBLOCK must be nonnegative; "//&
1991	0	"use 0 to select the default block size.")
1992	0	if (CS%nkblock < 0) &
1993		call MOM_error(FATAL, "ISOPYCNAL_NKBLOCK must be nonnegative; "//&
1994	0	"use 0 to select the default block size.")
1995
1996
1997	0	end subroutine ML_MEKE_init
1998
1999		!> Calculate the various features used for the machine learning prediction
2000	0	subroutine ML_MEKE_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, features_array)
2001		type(ocean_grid_type), intent(inout) :: G !< Ocean grid
2002		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2003		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2004		type(MEKE_CS), intent(in) :: CS !< Control structure for MEKE
2005		real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: Rd_dx_h !< Rossby radius of deformation over
2006		!! the grid length scale [nondim]
2007		real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
2008		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
2009		type(thermo_var_ptrs), intent(in) :: tv !< Type containing thermodynamic variables
2010		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
2011		real, intent(in) :: dt !< Model(baroclinic) time-step [T ~> s].
2012		real(kind=real32), dimension(SIZE(h),num_features), intent( out) :: features_array
2013		!< The array of features needed for machine
2014		!! learning inference, with different units
2015		!! for the various subarrays [various]
2016
2017	0	real, dimension(SZI_(G),SZJ_(G)) :: mke ! Surface kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
2018	0	real, dimension(SZI_(G),SZJ_(G)) :: slope_z ! Vertically averaged isoneutral slopes [Z L-1 ~> nondim]
2019	0	real, dimension(SZIB_(G),SZJB_(G)) :: rv_z ! Surface relative vorticity [T-1 ~> s-1]
2020	0	real, dimension(SZIB_(G),SZJB_(G)) :: rv_z_t ! Surface relative vorticity interpolated to tracer points [T-1 ~> s-1]
2021
2022	0	real, dimension(SZIB_(G),SZJ_(G), SZK_(G)) :: h_u ! Thickness at u point [H ~> m or kg m-2]
2023	0	real, dimension(SZI_(G),SZJB_(G), SZK_(G)) :: h_v ! Thickness at v point [H ~> m or kg m-2]
2024	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: slope_x ! Isoneutral slope at U point [Z L-1 ~> nondim]
2025	0	real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: slope_y ! Isoneutral slope at V point [Z L-1 ~> nondim]
2026	0	real, dimension(SZIB_(G),SZJ_(G)) :: slope_x_vert_avg ! Isoneutral slope at U point [Z L-1 ~> nondim]
2027	0	real, dimension(SZI_(G),SZJB_(G)) :: slope_y_vert_avg ! Isoneutral slope at V point [Z L-1 ~> nondim]
2028	0	real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: e ! The interface heights relative to mean sea level [Z ~> m].
2029		real :: slope_t ! Slope interpolated to thickness points [Z L-1 ~> nondim]
2030		real :: u_t, v_t ! u and v interpolated to thickness points [L T-1 ~> m s-1]
2031		real :: dvdx, dudy ! Components of relative vorticity [T-1 ~> s-1]
2032		real :: a_e, a_w, a_n, a_s ! Fractional areas of neighboring cells for interpolating velocities [nondim]
2033		real :: Idenom ! A normalizing factor in calculating weighted averages of areas [L-2 ~> m-2]
2034		real :: sum_area ! A sum of adjacent cell areas [L2 ~> m2]
2035
2036		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
2037		integer :: niblock, njblock, nkblock
2038
2039	0	niblock = CS%niblock
2040	0	njblock = CS%njblock
2041	0	nkblock = CS%nkblock
2042
2043	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2044	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
2045
2046		! Calculate various features for used to infer eddy kinetic energy
2047		! Linear interpolation to estimate thickness at a velocity points
2048	0	do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
2049	0	h_u(I,j,k) = 0.5(h(i,j,k)G%mask2dT(i,j) + h(i+1,j,k)*G%mask2dT(i+1,j)) + GV%Angstrom_H
2050	0	h_v(i,J,k) = 0.5(h(i,j,k)G%mask2dT(i,j) + h(i,j+1,k)*G%mask2dT(i,j+1)) + GV%Angstrom_H
2051		enddo ; enddo ; enddo
2052
2053	0	if (niblock == 0) niblock = ie - is + 1
2054	0	if (njblock == 0) njblock = je - js + 1
2055	0	if (nkblock == 0) nkblock = nz
2056
2057		!$omp target update to(h)
2058		!$omp target enter data map(alloc: e)
2059	0	call find_eta(h, tv, G, GV, US, e, halo_size=2)
2060		! Note the hard-coded dimenisional constant in the following line.
2061		! UMW: Below is untested
2062		!$omp target enter data map(to: tv, tv%T, tv%S, slope_x, slope_y)
2063		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
2064		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
2065		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt1.e-7GV%m2_s_to_HZ_T, .false., slope_x, slope_y, &
2066	0	niblock, njblock, nkblock )
2067		!$omp target exit data map(release: tv, tv%T, tv%S)
2068		!$omp target exit data map(delete: e)
2069		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
2070		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
2071		!$omp target exit data map(from: slope_x, slope_y)
2072	0	call pass_vector(slope_x, slope_y, G%Domain)
2073	0	do j=js-1,je+1 ; do i=is-1,ie+1
2074	0	slope_x_vert_avg(I,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), GV%H_subroundoff)
2075	0	slope_y_vert_avg(i,J) = vertical_average_interface(slope_y(i,j,:), h_v(i,j,:), GV%H_subroundoff)
2076		enddo ; enddo
2077	0	slope_z(:,:) = 0.
2078
2079	0	call pass_vector(slope_x_vert_avg, slope_y_vert_avg, G%Domain)
2080	0	do j=js,je ; do i=is,ie
2081		! Calculate weights for interpolation from velocity points to h points
2082	0	sum_area = G%areaCu(I-1,j) + G%areaCu(I,j)
2083	0	if (sum_area>0.0) then
2084	0	Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
2085	0	a_w = G%areaCu(I-1,j) * Idenom
2086	0	a_e = G%areaCu(I,j) * Idenom
2087		else
2088	0	a_w = 0.0 ; a_e = 0.0
2089		endif
2090
2091	0	sum_area = G%areaCv(i,J-1) + G%areaCv(i,J)
2092	0	if (sum_area>0.0) then
2093	0	Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
2094	0	a_s = G%areaCv(i,J-1) * Idenom
2095	0	a_n = G%areaCv(i,J) * Idenom
2096		else
2097	0	a_s = 0.0 ; a_n = 0.0
2098		endif
2099
2100		! Calculate mean kinetic energy
2101	0	u_t = (a_eu(I,j,1)) + (a_wu(I-1,j,1))
2102	0	v_t = (a_nv(i,J,1)) + (a_sv(i,J-1,1))
2103	0	mke(i,j) = 0.5( (u_tu_t) + (v_t*v_t) )
2104
2105		! Calculate the magnitude of the slope
2106	0	slope_t = slope_x_vert_avg(I,j)a_e+slope_x_vert_avg(I-1,j)a_w
2107	0	slope_z(i,j) = sqrt(slope_t*slope_t)
2108	0	slope_t = slope_y_vert_avg(i,J)a_n+slope_y_vert_avg(i,J-1)a_s
2109	0	slope_z(i,j) = 0.5(slope_z(i,j) + sqrt(slope_tslope_t))*G%mask2dT(i,j)
2110		enddo ; enddo
2111	0	call pass_var(slope_z, G%Domain)
2112
2113		! Calculate relative vorticity
2114	0	do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
2115	0	dvdx = ((v(i+1,J,1)G%dyCv(i+1,J)) - (v(i,J,1)G%dyCv(i,J)))
2116	0	dudy = ((u(I,j+1,1)G%dxCu(I,j+1)) - (u(I,j,1)G%dxCu(I,j)))
2117		! Assumed no slip
2118	0	rv_z(I,J) = (2.0-G%mask2dBu(I,J)) * (dvdx - dudy) * G%IareaBu(I,J)
2119		enddo ; enddo
2120		! Interpolate RV to t-point, revisit this calculation to include metrics
2121	0	do j=js,je ; do i=is,ie
2122	0	rv_z_t(i,j) = 0.25*(rv_z(i-1,j) + rv_z(i,j) + rv_z(i-1,j-1) + rv_z(i,j-1))
2123		enddo ; enddo
2124
2125
2126		! Construct the feature array
2127	0	features_array(:,mke_idx) = pack(mke,.true.)
2128	0	features_array(:,slope_z_idx) = pack(slope_z,.true.)
2129	0	features_array(:,rd_dx_z_idx) = pack(Rd_dx_h,.true.)
2130	0	features_array(:,rv_idx) = pack(rv_z_t,.true.)
2131
2132	0	if (CS%id_rv>0) call post_data(CS%id_rv, rv_z, CS%diag)
2133	0	if (CS%id_mke>0) call post_data(CS%id_mke, mke, CS%diag)
2134	0	if (CS%id_slope_z>0) call post_data(CS%id_slope_z, slope_z, CS%diag)
2135	0	if (CS%id_slope_x>0) call post_data(CS%id_slope_x, slope_x, CS%diag)
2136	0	if (CS%id_slope_y>0) call post_data(CS%id_slope_y, slope_y, CS%diag)
2137	0	end subroutine ML_MEKE_calculate_features
2138
2139		!> Use the machine learning interface to predict EKE
2140	0	subroutine predict_MEKE(G, US, CS, npts, Time, features_array, MEKE)
2141		type(ocean_grid_type), intent(inout) :: G !< Ocean grid
2142		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2143		type(MEKE_CS), intent(in ) :: CS !< Control structure for MEKE
2144		integer, intent(in ) :: npts !< Number of T-grid cells on the local
2145		!! domain
2146		type(time_type), intent(in ) :: Time !< The current model time
2147		real(kind=real32), dimension(npts,num_features), intent(in ) :: features_array
2148		!< The array of features needed for machine
2149		!! learning inference, with different units
2150		!! for the various subarrays [various]
2151		real, dimension(SZI_(G),SZJ_(G)), intent( out) :: MEKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2]
2152
2153		! Local variables
2154		integer :: db_return_code
2155		character(len=255), dimension(1) :: model_out, model_in
2156		character(len=255) :: time_suffix
2157	0	real(kind=real32), dimension(SIZE(MEKE)) :: MEKE_vec ! A one-dimensional array of the natural log of eddy kinetic
2158		! energy in mks units [m2 s-2]
2159	0	real, dimension(size(MEKE,1),size(MEKE,2)) :: ln_MEKE ! the natural log of eddy kinetic energy
2160		! in mks units [m2 s-2]
2161	0	real, dimension(size(MEKE,1),size(MEKE,2)) :: MEKE_mks ! The eddy kinetic energy in mks units [m2 s-2]
2162		integer :: i, j, is, ie, js, je
2163
2164	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
2165		!> Use the database client to call a machine learning model to predict eddy kinetic energy
2166	0	call cpu_clock_begin(CS%id_put_tensor)
2167	0	db_return_code = CS%client%put_tensor("features"//CS%key_suffix, features_array, shape(features_array))
2168	0	call cpu_clock_end(CS%id_put_tensor)
2169
2170		! Run the ML model to predict EKE and return the result
2171	0	model_out(1) = "EKE"//CS%key_suffix
2172	0	model_in(1) = "features"//CS%key_suffix
2173	0	call cpu_clock_begin(CS%id_run_model)
2174	0	db_return_code = CS%client%run_model(CS%model_key, model_in, model_out)
2175	0	call cpu_clock_end(CS%id_run_model)
2176	0	if (CS%client%SR_error_parser(db_return_code)) then
2177	0	call MOM_error(FATAL, "MEKE: run_model failed")
2178		endif
2179	0	call cpu_clock_begin(CS%id_unpack_tensor)
2180	0	db_return_code = CS%client%unpack_tensor( model_out(1), MEKE_vec, shape(MEKE_vec) )
2181	0	call cpu_clock_end(CS%id_unpack_tensor)
2182
2183	0	ln_MEKE = reshape(MEKE_vec, shape(MEKE))
2184		! Zero out the halos. These will usually be reset by the pass_var in a few lines.
2185	0	MEKE_mks(:,:) = 0.0
2186	0	do j=js,je ; do i=is,ie
2187	0	MEKE_mks(i,j) = MIN(exp(ln_MEKE(i,j)), US%L_T_to_m_s*2CS%eke_max)
2188		enddo ; enddo
2189	0	call pass_var(MEKE_mks, G%Domain, halo=1)
2190
2191	0	if (CS%online_analysis) then
2192	0	write(time_suffix,"(F16.0)") time_type_to_real(Time)
2193		db_return_code = CS%client%put_tensor(trim("EKE_")//trim(adjustl(time_suffix))//CS%key_suffix, &
2194	0	MEKE_mks, shape(MEKE))
2195		endif
2196
2197		! Copy MEKE_mks into the argument in rescaled units.
2198		! MEKE(:,:) = 0.0 ! This would fill in the wider halos of this intent(out) array.
2199	0	do j=js-1,je+1 ; do i=is-1,ie+1
2200	0	MEKE(i,j) = US%m_s_to_L_T*2 MEKE_mks(i,j)
2201		enddo ; enddo
2202
2203	0	end subroutine predict_MEKE
2204
2205		!> Compute average of interface quantities weighted by the thickness of the surrounding
2206		!! layers [arbitrary]
2207	0	real function vertical_average_interface(h, w, h_min)
2208
2209		real, dimension(:), intent(in) :: h !< Layer Thicknesses [H ~> m or kg m-2]
2210		real, dimension(:), intent(in) :: w !< Quantity to average [arbitrary]
2211		real, intent(in) :: h_min !< The vanishingly small layer thickness [H ~> m or kg m-2]
2212
2213		real :: htot ! Twice the sum of the layer thicknesses interpolated to interior interfaces [H ~> m or kg m-2]
2214		real :: inv_htot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1]
2215		integer :: k, nk
2216
2217	0	nk = size(h)
2218	0	htot = h_min
2219	0	do k=2,nk
2220	0	htot = htot + (h(k-1)+h(k))
2221		enddo
2222	0	inv_htot = 1./htot
2223
2224	0	vertical_average_interface = 0.
2225	0	do K=2,nk
2226	0	vertical_average_interface = vertical_average_interface + (w(k)(h(k-1)+h(k)))inv_htot
2227		enddo
2228	0	end function vertical_average_interface
2229
2230		!> Allocates memory and register restart fields for the MOM_MEKE module.
2231	1	subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)
2232		! Arguments
2233		type(hor_index_type), intent(in) :: HI !< Horizontal index structure
2234		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2235		type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
2236		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
2237		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct
2238
2239		! Local variables
2240		real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_bhFrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim]
2241		real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim]
2242		logical :: Use_KH_in_MEKE
2243		logical :: useMEKE
2244		logical :: sqg_use_MEKE
2245		integer :: isd, ied, jsd, jed
2246
2247		! Determine whether this module will be used
2248	1	useMEKE = .false. ; call read_param(param_file,"USE_MEKE",useMEKE)
2249
2250		! Read these parameters to determine what should be in the restarts
2251	1	MEKE_GMcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff)
2252	1	MEKE_FrCoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff)
2253	1	MEKE_bhFrCoeff = -1. ; call read_param(param_file,"MEKE_bhFRCOEFF",MEKE_bhFrCoeff)
2254	1	MEKE_GMEcoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff)
2255	1	MEKE_KhCoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff)
2256	1	MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku)
2257	1	MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au)
2258	1	Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE)
2259	1	sqg_use_MEKE = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_MEKE)
2260
2261	1	if (.not. useMEKE) return
2262
2263		! Allocate memory
2264	1	call MOM_mesg("MEKE_alloc_register_restart: allocating and registering", 5)
2265	1	isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed
2266	8773	allocate(MEKE%MEKE(isd:ied,jsd:jed), source=0.0)
2267		call register_restart_field(MEKE%MEKE, "MEKE", .false., restart_CS, &
2268	1	longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2)
2269
2270	8773	if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0)
2271	1	if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) &
2272	0	allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0)
2273	1	if (MEKE_bhFrCoeff >= 0.) &
2274	0	allocate(MEKE%mom_src_bh(isd:ied,jsd:jed), source=0.0)
2275	1	if (MEKE_FrCoeff<0.) MEKE_FrCoeff = 0.
2276	1	if (MEKE_bhFrCoeff<0.) MEKE_bhFrCoeff = 0.
2277	1	if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0)
2278	1	if (MEKE_KhCoeff>=0.) then
2279	8773	allocate(MEKE%Kh(isd:ied,jsd:jed), source=0.0)
2280		call register_restart_field(MEKE%Kh, "MEKE_Kh", .false., restart_CS, &
2281		longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy", &
2282	1	units="m2 s-1", conversion=US%L_to_m*2US%s_to_T)
2283		endif
2284	8773	allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
2285	1	if (MEKE_viscCoeff_Ku/=0.) then
2286	0	allocate(MEKE%Ku(isd:ied,jsd:jed), source=0.0)
2287		call register_restart_field(MEKE%Ku, "MEKE_Ku", .false., restart_CS, &
2288		longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &
2289	0	units="m2 s-1", conversion=US%L_to_m*2US%s_to_T)
2290		endif
2291	1	if (sqg_use_MEKE) then
2292	0	allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0)
2293		call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, &
2294		longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &
2295	0	units="m", conversion=US%L_to_m)
2296		endif
2297	1	if (Use_Kh_in_MEKE) then
2298	0	allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0)
2299		call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, &
2300		longname="Copy of thickness diffusivity for diffusing MEKE", &
2301	0	units="m2 s-1", conversion=US%L_to_m*2US%s_to_T)
2302		endif
2303
2304	1	if (MEKE_viscCoeff_Au/=0.) then
2305	0	allocate(MEKE%Au(isd:ied,jsd:jed), source=0.0)
2306		call register_restart_field(MEKE%Au, "MEKE_Au", .false., restart_CS, &
2307		longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy", &
2308	0	units="m4 s-1", conversion=US%L_to_m*4US%s_to_T)
2309		endif
2310
2311		end subroutine MEKE_alloc_register_restart
2312
2313		!> Deallocates any variables allocated in MEKE_alloc_register_restart.
2314	1	subroutine MEKE_end(MEKE)
2315		type(MEKE_type), intent(inout) :: MEKE !< A structure with MEKE-related fields.
2316
2317		! NOTE: MEKE will always be allocated by MEKE_init, even if MEKE is disabled.
2318		! So these must all be conditional, even though MEKE%MEKE and MEKE%Rd_dx_h
2319		! are always allocated (when MEKE is enabled)
2320
2321		! Unmap before deallocating, mirroring the maps at the end of MEKE_init.
2322		if (allocated(MEKE%MEKE)) then
2323		!$omp target exit data map(delete: MEKE%MEKE)
2324		endif
2325		if (allocated(MEKE%GM_src)) then
2326		!$omp target exit data map(delete: MEKE%GM_src)
2327		endif
2328		if (allocated(MEKE%mom_src)) then
2329		!$omp target exit data map(delete: MEKE%mom_src)
2330		endif
2331		if (allocated(MEKE%mom_src_bh)) then
2332		!$omp target exit data map(delete: MEKE%mom_src_bh)
2333		endif
2334		if (allocated(MEKE%GME_snk)) then
2335		!$omp target exit data map(delete: MEKE%GME_snk)
2336		endif
2337		if (allocated(MEKE%Kh)) then
2338		!$omp target exit data map(delete: MEKE%Kh)
2339		endif
2340		if (allocated(MEKE%Kh_diff)) then
2341		!$omp target exit data map(delete: MEKE%Kh_diff)
2342		endif
2343		if (allocated(MEKE%Rd_dx_h)) then
2344		!$omp target exit data map(delete: MEKE%Rd_dx_h)
2345		endif
2346		if (allocated(MEKE%Ku)) then
2347		!$omp target exit data map(delete: MEKE%Ku)
2348		endif
2349		if (allocated(MEKE%Au)) then
2350		!$omp target exit data map(delete: MEKE%Au)
2351		endif
2352		if (allocated(MEKE%Le)) then
2353		!$omp target exit data map(delete: MEKE%Le)
2354		endif
2355
2356	1	if (allocated(MEKE%Au)) deallocate(MEKE%Au)
2357	1	if (allocated(MEKE%Kh_diff)) deallocate(MEKE%Kh_diff)
2358	1	if (allocated(MEKE%Ku)) deallocate(MEKE%Ku)
2359	1	if (allocated(MEKE%Rd_dx_h)) deallocate(MEKE%Rd_dx_h)
2360	1	if (allocated(MEKE%Kh)) deallocate(MEKE%Kh)
2361	1	if (allocated(MEKE%GME_snk)) deallocate(MEKE%GME_snk)
2362	1	if (allocated(MEKE%mom_src)) deallocate(MEKE%mom_src)
2363	1	if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh)
2364	1	if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src)
2365	1	if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE)
2366	1	if (allocated(MEKE%Le)) deallocate(MEKE%Le)
2367	1	end subroutine MEKE_end
2368
2369		!> \namespace mom_meke
2370		!!
2371		!! \section section_MEKE The Mesoscale Eddy Kinetic Energy (MEKE) framework
2372		!!
2373		!! The MEKE framework accounts for the mean potential energy removed by
2374		!! the first order closures used to parameterize mesoscale eddies.
2375		!! It requires closure at the second order, namely dissipation and transport
2376		!! of eddy energy.
2377		!!
2378		!! Monitoring the sub-grid scale eddy energy budget provides a means to predict
2379		!! a sub-grid eddy-velocity scale which can be used in the lower order closures.
2380		!!
2381		!! \subsection section_MEKE_equations MEKE equations
2382		!!
2383		!! The eddy kinetic energy equation is:
2384		!! \f[ \partial_{\tilde{t}} E =
2385		!! \overbrace{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v
2386		!! }^\text{sources}
2387		!! - \overbrace{ ( \lambda + C_d \| U_d \| \gamma_b^2 ) E
2388		!! }^\text{local dissipation}
2389		!! + \overbrace{ \nabla \cdot ( ( \kappa_E + \gamma_M \kappa_M ) \nabla E
2390		!! - \kappa_4 \nabla^3 E )
2391		!! }^\text{smoothing}
2392		!! \f]
2393		!! where \f$ E \f$ is the eddy kinetic energy (variable <code>MEKE</code>) with units of
2394		!! m<sup>2</sup>s<sup>-2</sup>,
2395		!! and \f$\tilde{t} = a t\f$ is a scaled time. The non-dimensional factor
2396		!! \f$ a\geq 1 \f$ is used to accelerate towards equilibrium.
2397		!!
2398		!! The MEKE equation is two-dimensional and obtained by depth averaging the
2399		!! the three-dimensional eddy energy equation. In the following expressions
2400		!! \f$ \left< \phi \right> = \frac{1}{H} \int^\eta_{-D} \phi \, dz \f$ maps
2401		!! three dimensional terms into the two-dimensional quantities needed.
2402		!!
2403		!! \subsubsection section_MEKE_source_terms MEKE source terms
2404		!!
2405		!! The source term \f$ \dot{E}_b \f$ is a constant background source
2406		!! of energy intended to avoid the limit \f$E\rightarrow 0\f$.
2407		!!
2408		!! The "GM" source term
2409		!! \f[ \dot{E}_\eta = - \left< \overline{w^\prime b^\prime} \right>
2410		!! = \left< \kappa_h N^2S^2 \right>
2411		!! \approx \left< \kappa_h g\prime \|\nabla_\sigma \eta\|^2 \right>\f]
2412		!! equals the mean potential energy removed by the Gent-McWilliams closure,
2413		!! and is excluded/included in the MEKE budget by the efficiency parameter
2414		!! \f$ \gamma_\eta \in [0,1] \f$.
2415		!!
2416		!! The "frictional" source term
2417		!! \f[ \dot{E}_{v} = \left< \partial_i u_j \tau_{ij} \right> \f]
2418		!! equals the mean kinetic energy removed by lateral viscous fluxes, and
2419		!! is excluded/included in the MEKE budget by the efficiency parameter
2420		!! \f$ \gamma_v \in [0,1] \f$.
2421		!!
2422		!! \subsubsection section_MEKE_dissipation_terms MEKE dissipation terms
2423		!!
2424		!! The local dissipation of \f$ E \f$ is parameterized through a linear
2425		!! damping, \f$\lambda\f$, and bottom drag, \f$ C_d \| U_d \| \gamma_b^2 \f$.
2426		!! The \f$ \gamma_b \f$ accounts for the weak projection of the column-mean
2427		!! eddy velocity to the bottom. In other words, the bottom velocity is
2428		!! estimated as \f$ \gamma_b U_e \f$.
2429		!! The bottom drag coefficient, \f$ C_d \f$ is the same as that used in the bottom
2430		!! friction in the mean model equations.
2431		!!
2432		!! The bottom drag velocity scale, \f$ U_d \f$, has contributions from the
2433		!! resolved state and \f$ E \f$:
2434		!! \f[ U_d = \sqrt{ U_b^2 + \|u\|^2_{z=-D} + \|\gamma_b U_e\|^2 } .\f]
2435		!! where the eddy velocity scale, \f$ U_e \f$, is given by:
2436		!! \f[ U_e = \sqrt{ 2 E } .\f]
2437		!! \f$ U_b \f$ is a constant background bottom velocity scale and is
2438		!! typically not used (i.e. set to zero).
2439		!!
2440		!! Following \cite jansen2015, the projection of eddy energy on to the bottom
2441		!! is given by the ratio of bottom energy to column mean energy:
2442		!! \f[
2443		!! \gamma_b^2 = \frac{E_b}{E} = \gamma_{d0}
2444		!! + \left( 1 + c_{b} \frac{L_d}{L_f} \right)^{-\frac{4}{5}}
2445		!! ,
2446		!! \f]
2447		!! \f[
2448		!! \gamma_b^2 \leftarrow \max{\left( \gamma_b^2, \gamma_{min}^2 \right)}
2449		!! .
2450		!! \f]
2451		!!
2452		!! \subsection section_MEKE_smoothing MEKE smoothing terms
2453		!!
2454		!! \f$ E \f$ is laterally diffused by a diffusivity \f$ \kappa_E + \gamma_M
2455		!! \kappa_M \f$ where \f$ \kappa_E \f$ is a constant diffusivity and the term
2456		!! \f$ \gamma_M \kappa_M \f$ is a "self diffusion" using the diffusivity
2457		!! calculated in the section \ref section_MEKE_diffusivity.
2458		!! \f$ \kappa_4 \f$ is a constant bi-harmonic diffusivity.
2459		!!
2460		!! \subsection section_MEKE_diffusivity Diffusivity derived from MEKE
2461		!!
2462		!! The predicted eddy velocity scale, \f$ U_e \f$, can be combined with a
2463		!! mixing length scale to form a diffusivity.
2464		!! The primary use of a MEKE derived diffusivity is for use in thickness
2465		!! diffusion (module mom_thickness_diffuse) and optionally in along
2466		!! isopycnal mixing of tracers (module mom_tracer_hor_diff).
2467		!! The original form used (enabled with MEKE_OLD_LSCALE=True):
2468		!!
2469		!! \f[ \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \f]
2470		!!
2471		!! where \f$ A_\Delta \f$ is the area of the grid cell.
2472		!! Following \cite jansen2015, we now use
2473		!!
2474		!! \f[ \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \f]
2475		!!
2476		!! where \f$ \gamma_\kappa \in [0,1] \f$ is a non-dimensional factor and,
2477		!! following \cite jansen2015, \f$\gamma_t^2\f$ is the ratio of barotropic
2478		!! eddy energy to column mean eddy energy given by
2479		!! \f[
2480		!! \gamma_t^2 = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}}
2481		!! ,
2482		!! \f]
2483		!! \f[
2484		!! \gamma_t^2 \leftarrow \max{\left( \gamma_t^2, \gamma_{min}^2 \right)}
2485		!! .
2486		!! \f]
2487		!!
2488		!! The length-scale is a configurable combination of multiple length scales:
2489		!!
2490		!! \f[
2491		!! l_M = \left(
2492		!! \frac{\alpha_d}{L_d}
2493		!! + \frac{\alpha_f}{L_f}
2494		!! + \frac{\alpha_R}{L_R}
2495		!! + \frac{\alpha_e}{L_e}
2496		!! + \frac{\alpha_\Delta}{L_\Delta}
2497		!! + \frac{\delta[L_c]}{L_c}
2498		!! \right)^{-1}
2499		!! \f]
2500		!!
2501		!! where
2502		!!
2503		!! \f{eqnarray*}{
2504		!! L_d & = & \sqrt{\frac{c_g^2}{f^2+2\beta c_g}} \sim \frac{ c_g }{f} \\\\
2505		!! L_R & = & \sqrt{\frac{U_e}{\beta^*}} \\\\
2506		!! L_e & = & \frac{U_e}{\|S\| N} \\\\
2507		!! L_f & = & \frac{H}{c_d} \\\\
2508		!! L_\Delta & = & \sqrt{A_\Delta} .
2509		!! \f}
2510		!!
2511		!! \f$L_c\f$ is a constant and \f$\delta[L_c]\f$ is the impulse function so that the term
2512		!! \f$\frac{\delta[L_c]}{L_c}\f$ evaluates to \f$\frac{1}{L_c}\f$ when \f$L_c\f$ is non-zero
2513		!! but is dropped if \f$L_c=0\f$.
2514		!!
2515		!! \f$\beta^*\f$ is the effective \f$\beta\f$ that combines both the planetary vorticity
2516		!! gradient (i.e. \f$\beta=\nabla f\f$) and the topographic \f$\beta\f$ effect,
2517		!! with the latter weighed by a weighting constant, \f$c_\beta\f$, that varies
2518		!! from 0 to 1, so that \f$c_\beta=0\f$ means the topographic \f$\beta\f$ effect is ignored,
2519		!! while \f$c_\beta=1\f$ means it is fully considered. The new \f$\beta^*\f$ therefore
2520		!! takes the form of
2521		!!
2522		!! \f[
2523		!! \beta^* = \sqrt{( \partial_xf - c_\beta\frac{f}{D}\partial_xD )^2 +
2524		!! ( \partial_yf - c_\beta\frac{f}{D}\partial_yD )^2}
2525		!! \f]
2526		!! where \f$D\f$ is water column depth at T points.
2527		!!
2528		!! \subsection section_MEKE_viscosity Viscosity derived from MEKE
2529		!!
2530		!! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be
2531		!! used to form a harmonic eddy viscosity,
2532		!!
2533		!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } \f]
2534		!!
2535		!! as well as a biharmonic eddy viscosity,
2536		!!
2537		!! \f[ \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 } \f]
2538		!!
2539		!! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance
2540		!!
2541		!! Note that in steady-state (or when \f$ a>>1 \f$) and there is no
2542		!! diffusion of \f$ E \f$ then
2543		!! \f[ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2544		!! \gamma_v \dot{E}_v }{ \lambda + C_d\|U_d\|\gamma_b^2 } . \f]
2545		!!
2546		!! In the linear drag limit, where
2547		!! \f$ U_e << \min(U_b, \|u\|_{z=-D}, C_d^{-1}\lambda) \f$, the equilibrium becomes
2548		!! \f$ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2549		!! \gamma_v \dot{E}_v }{ \lambda + C_d \sqrt{ U_b^2 + \|u\|^2_{z=-D} } } \f$.
2550		!!
2551		!! In the nonlinear drag limit, where \f$ U_e >> \max(U_b, \|u\|_{z=-D}, C_d^{-1}\lambda) \f$,
2552		!! the equilibrium becomes
2553		!! \f$ \overline{E} \approx \left( \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
2554		!! \gamma_v \dot{E}_v }{ \sqrt{2} C_d \gamma_b^3 } \right)^\frac{2}{3} \f$.
2555		!!
2556		!! \subsubsection section_MEKE_module_parameters MEKE module parameters
2557		!!
2558		!! \| Symbol \| Module parameter \|
2559		!! \| ------ \| --------------- \|
2560		!! \| - \| <code>USE_MEKE</code> \|
2561		!! \| \f$ a \f$ \| <code>MEKE_DTSCALE</code> \|
2562		!! \| \f$ \dot{E}_b \f$ \| <code>MEKE_BGSRC</code> \|
2563		!! \| \f$ \gamma_\eta \f$ \| <code>MEKE_GMCOEFF</code> \|
2564		!! \| \f$ \gamma_v \f$ \| <code>MEKE_FrCOEFF</code> \|
2565		!! \| \f$ \lambda \f$ \| <code>MEKE_DAMPING</code> \|
2566		!! \| \f$ U_b \f$ \| <code>MEKE_USCALE</code> \|
2567		!! \| \f$ \gamma_{d0} \f$ \| <code>MEKE_CD_SCALE</code> \|
2568		!! \| \f$ c_{b} \f$ \| <code>MEKE_CB</code> \|
2569		!! \| \f$ c_{t} \f$ \| <code>MEKE_CT</code> \|
2570		!! \| \f$ \kappa_E \f$ \| <code>MEKE_KH</code> \|
2571		!! \| \f$ \kappa_4 \f$ \| <code>MEKE_K4</code> \|
2572		!! \| \f$ \gamma_\kappa \f$ \| <code>MEKE_KHCOEFF</code> \|
2573		!! \| \f$ \gamma_M \f$ \| <code>MEKE_KHMEKE_FAC</code> \|
2574		!! \| \f$ \gamma_u \f$ \| <code>MEKE_VISCOSITY_COEFF_KU</code> \|
2575		!! \| \f$ \gamma_4 \f$ \| <code>MEKE_VISCOSITY_COEFF_AU</code> \|
2576		!! \| \f$ \gamma_{min}^2 \f$\| <code>MEKE_MIN_GAMMA2</code> \|
2577		!! \| \f$ \alpha_d \f$ \| <code>MEKE_ALPHA_DEFORM</code> \|
2578		!! \| \f$ \alpha_f \f$ \| <code>MEKE_ALPHA_FRICT</code> \|
2579		!! \| \f$ \alpha_R \f$ \| <code>MEKE_ALPHA_RHINES</code> \|
2580		!! \| \f$ \alpha_e \f$ \| <code>MEKE_ALPHA_EADY</code> \|
2581		!! \| \f$ \alpha_\Delta \f$ \| <code>MEKE_ALPHA_GRID</code> \|
2582		!! \| \f$ L_c \f$ \| <code>MEKE_FIXED_MIXING_LENGTH</code> \|
2583		!! \| \f$ c_\beta \f$ \| <code>MEKE_TOPOGRAPHIC_BETA</code> \|
2584		!! \| - \| <code>MEKE_KHTH_FAC</code> \|
2585		!! \| - \| <code>MEKE_KHTR_FAC</code> \|
2586		!!
2587		!! \| Symbol \| Model parameter \|
2588		!! \| ------ \| --------------- \|
2589		!! \| \f$ C_d \f$ \| <code>CDRAG</code> \|
2590		!!
2591		!! \subsection section_MEKE_references References
2592		!!
2593		!! Jansen, M. F., A. J. Adcroft, R. Hallberg, and I. M. Held, 2015: Parameterization of eddy fluxes based on a
2594		!! mesoscale energy budget. Ocean Modelling, 92, 28--41, http://doi.org/10.1016/j.ocemod.2015.05.007 .
2595		!!
2596		!! Marshall, D. P., and A. J. Adcroft, 2010: Parameterization of ocean eddies: Potential vorticity mixing, energetics
2597		!! and Arnold first stability theorem. Ocean Modelling, 32, 188--204, http://doi.org/10.1016/j.ocemod.2010.02.001 .
2598
2599	0	end module MOM_MEKE
2600