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