MOM_MEKE.F90

! This file is part of MOM6, the Modular Ocean Model version 6.

! See the LICENSE file for licensing information.

! SPDX-License-Identifier: Apache-2.0

!> Implements the Mesoscale Eddy Kinetic Energy framework

!! with topographic beta effect included in computing beta in Rhines scale

module mom_meke

use iso_fortran_env,       only : real32

use mom_coms,              only : pe_here

use mom_database_comms,    only : dbclient_type, dbcomms_cs_type

use mom_debugging,         only : hchksum, uvchksum

use mom_cpu_clock,         only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine

use mom_diag_mediator,     only : post_data, register_diag_field, safe_alloc_ptr

use mom_diag_mediator,     only : diag_ctrl, time_type

use mom_domains,           only : create_group_pass, do_group_pass, group_pass_type

use mom_domains,           only : pass_vector, pass_var

use mom_error_handler,     only : mom_error, fatal, warning, note, mom_mesg, is_root_pe

use mom_file_parser,       only : read_param, get_param, log_version, param_file_type

use mom_grid,              only : ocean_grid_type

use mom_hor_index,         only : hor_index_type

use mom_interface_heights, only : find_eta

use mom_interpolate,       only : init_external_field, time_interp_external

use mom_interpolate,       only : time_interp_external_init

use mom_interpolate,       only : external_field

use mom_io,                only : vardesc, var_desc, slasher

use mom_isopycnal_slopes,  only : calc_isoneutral_slopes

use mom_restart,           only : mom_restart_cs, register_restart_field, query_initialized

use mom_string_functions,  only : lowercase

use mom_time_manager,      only : time_type_to_real

use mom_unit_scaling,      only : unit_scale_type

use mom_variables,         only : vertvisc_type, thermo_var_ptrs

use mom_verticalgrid,      only : verticalgrid_type

use mom_meke_types,        only : meke_type

implicit none ; private

#include <MOM_memory.h>

public step_forward_meke, meke_init, meke_alloc_register_restart, meke_end

! Constants for this module

integer, parameter :: num_features = 4 !< How many features used to predict EKE

integer, parameter :: mke_idx = 1     !< Index of mean kinetic energy in the feature array

integer, parameter :: slope_z_idx = 2 !< Index of vertically averaged isopycnal slope in the feature array

integer, parameter :: rv_idx = 3      !< Index of surface relative vorticity in the feature array

integer, parameter :: rd_dx_z_idx = 4 !< Index of the radius of deformation over the grid size in the feature array

integer, parameter :: eke_prog = 1     !< Use prognostic equation to calculate EKE

integer, parameter :: eke_file = 2     !< Read in EKE from a file

integer, parameter :: eke_dbclient = 3 !< Infer EKE using a neural network

!> Control structure that contains MEKE parameters and diagnostics handles

type, public :: meke_cs ; private

  logical :: initialized = .false. !< True if this control structure has been initialized.

  ! Parameters

  real :: meke_frcoeff  !< Efficiency of conversion of ME into MEKE [nondim]

  real :: meke_bhfrcoeff!< Efficiency of conversion of ME into MEKE by the biharmonic dissipation [nondim]

  real :: meke_gmcoeff  !< Efficiency of conversion of PE into MEKE [nondim]

  real :: meke_gmecoeff !< Efficiency of conversion of MEKE into ME by GME [nondim]

  real :: meke_damping  !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1].

  real :: meke_cd_scale !< The ratio of the bottom eddy velocity to the column mean

                        !! eddy velocity, i.e. sqrt(2*MEKE), [nondim]. This should be less than 1

                        !! to account for the surface intensification of MEKE.

  real :: meke_cb       !< Coefficient in the \f$\gamma_{bot}\f$ expression [nondim]

  real :: meke_min_gamma!< Minimum value of gamma_b^2 allowed [nondim]

  real :: meke_ct       !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim]

  logical :: visc_drag  !< If true use the vertvisc_type to calculate bottom drag.

  logical :: meke_geometric !< If true, uses the GM coefficient formulation from the GEOMETRIC

                        !! framework (Marshall et al., 2012)

  real    :: meke_geometric_alpha !< The nondimensional coefficient governing the efficiency of the

                        !! GEOMETRIC thickness diffusion [nondim].

  logical :: meke_equilibrium_alt !< If true, use an alternative calculation for the

                        !! equilibrium value of MEKE.

  logical :: meke_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value,

                        !!  which is calculated at each time step.

  logical :: gm_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather

                        !! than the streamfunction for the MEKE GM source term.

  real    :: meke_min_depth_tot  !< The minimum total thickness over which to distribute MEKE energy

                        !! sources from GM energy conversion [H ~> m or kg m-2].  When the total

                        !! thickness is less than this, the sources are scaled away.

  logical :: rd_as_max_scale !< If true the length scale can not exceed the

                        !! first baroclinic deformation radius.

  logical :: use_old_lscale !< Use the old formula for mixing length scale.

  logical :: use_min_lscale !< Use simple minimum for mixing length scale.

  logical :: meke_positive  !< If true, it guarantees that MEKE will always be >= 0.

  real :: lscale_maxval !< The ceiling on the MEKE mixing length scale when use_min_lscale is true [L ~> m].

  real :: cdrag         !< The bottom drag coefficient for MEKE, times rescaling factors [H L-1 ~> nondim or kg m-3]

  real :: meke_bgsrc    !< Background energy source for MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).

  real :: meke_dtscale  !< Scale factor to accelerate time-stepping [nondim]

  real :: meke_khcoeff  !< Scaling factor to convert MEKE into Kh [nondim]

  real :: meke_uscale   !< MEKE velocity scale for bottom drag [L T-1 ~> m s-1]

  real :: meke_kh       !< Background lateral diffusion of MEKE [L2 T-1 ~> m2 s-1]

  real :: meke_k4       !< Background bi-harmonic diffusivity (of MEKE) [L4 T-1 ~> m4 s-1]

  real :: khmeke_fac    !< A factor relating MEKE%Kh to the diffusivity used for

                        !! MEKE itself [nondim].

  real :: viscosity_coeff_ku !< The scaling coefficient in the expression for

                        !! viscosity used to parameterize lateral harmonic momentum mixing

                        !! by unresolved eddies represented by MEKE [nondim].

  real :: viscosity_coeff_au !< The scaling coefficient in the expression for

                        !! viscosity used to parameterize lateral biharmonic momentum mixing

                        !! by unresolved eddies represented by MEKE [nondim].

  real :: lfixed        !< Fixed mixing length scale [L ~> m].

  real :: adeform       !< Weighting towards deformation scale of mixing length [nondim]

  real :: arhines       !< Weighting towards Rhines scale of mixing length [nondim]

  real :: africt        !< Weighting towards frictional arrest scale of mixing length [nondim]

  real :: aeady         !< Weighting towards Eady scale of mixing length [nondim]

  real :: agrid         !< Weighting towards grid scale of mixing length [nondim]

  real :: meke_advection_factor !< A scaling in front of the advection of MEKE [nondim]

  real :: meke_topographic_beta !< Weight for how much topographic beta is considered

                                !! when computing beta in Rhines scale [nondim]

  real :: meke_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its

                        !! equilibrium value [T-1 ~> s-1].

  logical :: meke_advection_bug !< If true, recover a bug in the calculation of the barotropic

                        !! transport for the advection of MEKE, wherein only the transports in the

                        !! deepest layer are used.

  logical :: fixed_total_depth  !< If true, use the nominal bathymetric depth as the estimate of

                        !! the time-varying ocean depth.  Otherwise base the depth on the total

                        !! ocean mass per unit area.

  real :: rho_fixed_total_depth !< A density used to translate the nominal bathymetric depth into an

                        !! estimate of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH

                        !! is true [R ~> kg m-3]

  logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled.

  logical :: initialize !< If True, invokes a steady state solver to calculate MEKE.

  logical :: debug      !< If true, write out checksums of data for debugging

  integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default),

                     !! read in from a file, or inferred via a neural network

  logical :: sqg_use_meke !< If True, use MEKE%Le for the SQG vertical structure.

  type(diag_ctrl), pointer :: diag => null() !< A type that regulates diagnostics output

  !>@{ Diagnostic handles

  integer :: id_meke = -1, id_ue = -1, id_kh = -1, id_src = -1

  integer :: id_src_adv = -1, id_src_mom_k4 = -1, id_src_btm_drag = -1

  integer :: id_src_gm = -1, id_src_mom_lp = -1, id_src_mom_bh = -1

  integer :: id_ub = -1, id_ut = -1

  integer :: id_gm_src = -1, id_mom_src = -1, id_mom_src_bh = -1, id_gme_snk = -1, id_decay = -1

  integer :: id_khmeke_u = -1, id_khmeke_v = -1, id_ku = -1, id_au = -1

  integer :: id_le = -1, id_gamma_b = -1, id_gamma_t = -1

  integer :: id_lrhines = -1, id_leady = -1

  integer :: id_meke_equilibrium = -1

  !>@}

  type(external_field) :: eke_handle   !< Handle for reading in EKE from a file

  ! Infrastructure

  integer :: id_clock_pass !< Clock for group pass calls

  type(group_pass_type) :: pass_meke !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff

  type(group_pass_type) :: pass_kh   !< Group halo pass handle for MEKE%Kh, MEKE%Ku, and/or MEKE%Au

  ! MEKE via Machine Learning

  type(dbclient_type), pointer :: client => null() !< Pointer to the database client

  logical :: online_analysis !< If true, post the EKE used in MOM6 at every timestep

  character(len=5) :: model_key  = 'mleke'  !< Key where the ML-model is stored

  character(len=7) :: key_suffix !< Suffix appended to every key sent to Redis

  real :: eke_max !< The maximum value of EKE considered physically reasonable [L2 T-2 ~> m2 s-2]

  ! Clock ids

  integer :: id_client_init   !< Clock id to time initialization of the client

  integer :: id_put_tensor    !< Clock id to time put_tensor routine

  integer :: id_run_model     !< Clock id to time running of the ML model

  integer :: id_unpack_tensor !< Clock id to time retrieval of EKE prediction

  ! Diagnostic ids

  integer :: id_mke     = -1 !< Diagnostic id for surface mean kinetic energy

  integer :: id_slope_z = -1 !< Diagnostic id for vertically averaged horizontal slope magnitude

  integer :: id_slope_x = -1 !< Diagnostic id for isopycnal slope in the x-direction

  integer :: id_slope_y = -1 !< Diagnostic id for isopycnal slope in the y-direction

  integer :: id_rv      = -1 !< Diagnostic id for surface relative vorticity

end type meke_cs

contains

!> Integrates forward-in-time the MEKE eddy energy equation.

!! See \ref section_MEKE_equations.

subroutine step_forward_meke(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, hv, u, v, tv, Time)

  type(meke_type),                          intent(inout) :: meke !< MEKE data.

  type(ocean_grid_type),                    intent(inout) :: g    !< Ocean grid.

  type(verticalgrid_type),                  intent(in)    :: gv   !< Ocean vertical grid structure.

  type(unit_scale_type),                    intent(in)    :: us   !< A dimensional unit scaling type

  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)   :: h    !< Layer thickness [H ~> m or kg m-2].

  real, dimension(SZIB_(G),SZJ_(G)),        intent(in)    :: sn_u !< Eady growth rate at u-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJB_(G)),        intent(in)    :: sn_v !< Eady growth rate at v-points [T-1 ~> s-1].

  type(vertvisc_type),                      intent(in)    :: visc !< The vertical viscosity type.

  real,                                     intent(in)    :: dt   !< Model(baroclinic) time-step [T ~> s].

  type(meke_cs),                            intent(inout) :: cs   !< MEKE control structure.

  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in)  :: hu   !< Accumulated zonal mass flux [H L2 ~> m3 or kg].

  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in)  :: hv   !< Accumulated meridional mass flux [H L2 ~> m3 or kg]

  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u  !< Zonal velocity [L T-1 ~> m s-1]

  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v  !< Meridional velocity [L T-1 ~> m s-1]

  type(thermo_var_ptrs),                    intent(in)    :: tv   !< Type containing thermodynamic variables

  type(time_type),                          intent(in)    :: time !< The time used for interpolating EKE

  ! Local variables

  real, dimension(SZI_(G),SZJ_(G)) :: &

    data_eke, &     ! EKE from file [L2 T-2 ~> m2 s-2]

    mass, &         ! The total mass of the water column [R Z ~> kg m-2].

    i_mass, &       ! The inverse of mass [R-1 Z-1 ~> m2 kg-1].

    depth_tot, &    ! The depth of the water column [H ~> m or kg m-2].

    src, &          ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3).

    meke_decay, &   ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1].

    src_adv, &      ! The MEKE source/tendency from the horizontal advection of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).

    src_mom_k4, &   ! The MEKE source/tendency from the bihamornic of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).

    src_btm_drag, & ! The MEKE source/tendency from the bottom drag acting on MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).

    src_gm, &       ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3).

    src_mom_lp, &   ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).

    src_mom_bh, &   ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).

    damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].

    meke_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2].

    drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1]

    drag_rate, &    ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].

    del2meke, &     ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2].

    del4meke, &     ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].

    lmixscale, &    ! Eddy mixing length [L ~> m].

    barotrfac2, &   ! Ratio of EKE_barotropic / EKE [nondim]

    bottomfac2, &   ! Ratio of EKE_bottom / EKE [nondim]

    tmp, &          ! Temporary variable for computation of diagnostic velocities [L T-1 ~> m s-1]

    equilibrium_value, & ! The equilibrium value of MEKE to be calculated at

                    ! each time step [L2 T-2 ~> m2 s-2]

    damp_rate, &    ! The MEKE damping rate [T-1 ~> s-1]

    damping         ! The net damping of a field after sdt_damp [nondim]

  real, dimension(SZIB_(G),SZJ_(G)) :: &

    meke_uflux, &   ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].

                    ! In one place, MEKE_uflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].

    kh_u, &         ! The zonal diffusivity that is actually used [L2 T-1 ~> m2 s-1].

    barohu, &       ! Depth integrated accumulated zonal mass flux [R Z L2 ~> kg].

    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]

  real, dimension(SZI_(G),SZJB_(G)) :: &

    meke_vflux, &   ! The meridional advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m2 s-3].

                    ! In one place, MEKE_vflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].

    kh_v, &         ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1].

    barohv, &       ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg].

    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]

  real :: bh_coeff  ! Biharmonic part of efficiency conversion in total MEKE [nondim]

  real :: kh_here   ! The local horizontal viscosity [L2 T-1 ~> m2 s-1]

  real :: inv_kh_max ! The inverse of the local horizontal viscosity [T L-2 ~> s m-2]

  real :: k4_here   ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1]

  real :: inv_k4_max ! The inverse of the local horizontal biharmonic viscosity [T L-4 ~> s m-4]

  real :: cdrag2    ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]

  real :: advfac    ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]

  real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].

  real :: sdt       ! dt to use locally [T ~> s] (could be scaled to accelerate)

  real :: sdt_damp  ! dt for damping [T ~> s] (sdt could be split).

  real :: damp_step ! Size of damping timestep relative to sdt [nondim]

  logical :: use_drag_rate ! Flag to indicate drag_rate is finite

  logical :: any_damping_diags_s1 ! True if any damped diagnostics are enabled in first stage

  logical :: any_damping_diags  ! True if any damped diagnostics are enabled

  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz

  real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features

                                        ! needed for the machine learning inference, with different

                                        ! units for the various subarrays [various]

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  if (.not.cs%initialized) call mom_error(fatal, &

         "MOM_MEKE: Module must be initialized before it is used.")

  if ((cs%MEKE_Cd_scale > 0.0) .or. (cs%MEKE_Cb>0.) .or. cs%visc_drag) then

    use_drag_rate = .true.

  else

    use_drag_rate = .false.

  endif

  ! Only integrate the MEKE equations if MEKE is required.

  if (.not. allocated(meke%MEKE)) then

!   call MOM_error(FATAL, "MOM_MEKE: MEKE%MEKE is not associated!")

    return

  endif

  select case(cs%eke_src)

  case(eke_prog)

    if (cs%debug) then

      if (allocated(meke%mom_src)) &

        call hchksum(meke%mom_src, 'MEKE mom_src', g%HI, unscale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

      if (allocated(meke%mom_src_bh)) &

        call hchksum(meke%mom_src_bh, 'MEKE mom_src_bh', g%HI, unscale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

      if (allocated(meke%GME_snk)) &

        call hchksum(meke%GME_snk, 'MEKE GME_snk', g%HI, unscale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

      if (allocated(meke%GM_src)) &

        call hchksum(meke%GM_src, 'MEKE GM_src', g%HI, unscale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

      if (allocated(meke%MEKE)) &

        call hchksum(meke%MEKE, 'MEKE MEKE', g%HI, unscale=us%L_T_to_m_s**2)

      call uvchksum("MEKE SN_[uv]", sn_u, sn_v, g%HI, unscale=us%s_to_T, &

                    scalar_pair=.true.)

      call uvchksum("MEKE h[uv]", hu, hv, g%HI, haloshift=0, symmetric=.true., &

                    unscale=gv%H_to_m*us%L_to_m**2)

    endif

    sdt = dt*cs%MEKE_dtScale ! Scaled dt to use for time-stepping

    mass_neglect = gv%H_to_RZ * gv%H_subroundoff

    cdrag2 = cs%cdrag**2

    ! With a depth-dependent (and possibly strong) damping, it seems

    ! advisable to use Strang splitting between the damping and diffusion.

    damp_step = 1.

    if (cs%MEKE_KH >= 0. .or. cs%MEKE_K4 >= 0.) damp_step = 0.5

    sdt_damp = sdt * damp_step

    ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]

    if (cs%MEKE_advection_factor>0.) then

      do j=js,je ; do i=is-1,ie

        barohu(i,j) = 0.

      enddo ; enddo

      do k=1,nz

        do j=js,je ; do i=is-1,ie

          barohu(i,j) = barohu(i,j) + hu(i,j,k) * gv%H_to_RZ

        enddo ; enddo

      enddo

      do j=js-1,je ; do i=is,ie

        barohv(i,j) = 0.

      enddo ; enddo

      do k=1,nz

        do j=js-1,je ; do i=is,ie

          barohv(i,j) = barohv(i,j) + hv(i,j,k) * gv%H_to_RZ

        enddo ; enddo

      enddo

      if (cs%MEKE_advection_bug) then

        ! This obviously incorrect code reproduces a bug in the original implementation of

        ! the MEKE advection.

        do j=js,je ; do i=is-1,ie

          barohu(i,j) = hu(i,j,nz) * gv%H_to_RZ

        enddo ; enddo

        do j=js-1,je ; do i=is,ie

          barohv(i,j) = hv(i,j,nz) * gv%H_to_RZ

        enddo ; enddo

      endif

    endif

    ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow

    if (cs%visc_drag .and. allocated(visc%Kv_bbl_u) .and. allocated(visc%Kv_bbl_v)) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is-1,ie

        drag_vel_u(i,j) = 0.0

        if ((g%mask2dCu(i,j) > 0.0) .and. (visc%bbl_thick_u(i,j) > 0.0)) &

          drag_vel_u(i,j) = visc%Kv_bbl_u(i,j) / visc%bbl_thick_u(i,j)

      enddo ; enddo

      !$OMP parallel do default(shared)

      do j=js-1,je ; do i=is,ie

        drag_vel_v(i,j) = 0.0

        if ((g%mask2dCv(i,j) > 0.0) .and. (visc%bbl_thick_v(i,j) > 0.0)) &

          drag_vel_v(i,j) = visc%Kv_bbl_v(i,j) / visc%bbl_thick_v(i,j)

      enddo ; enddo

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        drag_rate_visc(i,j) = (0.25*g%IareaT(i,j) * &

                (((g%areaCu(i-1,j)*drag_vel_u(i-1,j)) + &

                  (g%areaCu(i,j)*drag_vel_u(i,j))) + &

                 ((g%areaCv(i,j-1)*drag_vel_v(i,j-1)) + &

                  (g%areaCv(i,j)*drag_vel_v(i,j))) ) )

      enddo ; enddo

    else

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        drag_rate_visc(i,j) = 0.

      enddo ; enddo

    endif

    !$OMP parallel do default(shared)

    do j=js-1,je+1

      do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo

      do k=1,nz ; do i=is-1,ie+1

        mass(i,j) = mass(i,j) + g%mask2dT(i,j) * (gv%H_to_RZ * h(i,j,k)) ! [R Z ~> kg m-2]

      enddo ; enddo

      do i=is-1,ie+1

        i_mass(i,j) = 0.0

        if (mass(i,j) > 0.0) i_mass(i,j) = 1.0 / mass(i,j) ! [R-1 Z-1 ~> m2 kg-1]

      enddo

    enddo

    if (cs%fixed_total_depth) then

      if (gv%Boussinesq) then

        !$OMP parallel do default(shared)

        do j=js-1,je+1 ; do i=is-1,ie+1

          depth_tot(i,j) = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) * gv%Z_to_H

        enddo ; enddo

      else

        !$OMP parallel do default(shared)

        do j=js-1,je+1 ; do i=is-1,ie+1

          depth_tot(i,j) = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) * cs%rho_fixed_total_depth * gv%RZ_to_H

        enddo ; enddo

      endif

    else

      !$OMP parallel do default(shared)

      do j=js-1,je+1 ; do i=is-1,ie+1

        depth_tot(i,j) = mass(i,j) * gv%RZ_to_H

      enddo ; enddo

    endif

    if (cs%initialize) then

      call meke_equilibrium(cs, meke, g, gv, us, sn_u, sn_v, drag_rate_visc, i_mass, depth_tot)

      cs%initialize = .false.

    endif

    ! Calculates bottomFac2, barotrFac2 and LmixScale

    call meke_lengthscales(cs, meke, g, gv, us, sn_u, sn_v, meke%MEKE, depth_tot, bottomfac2, barotrfac2, lmixscale)

    if (cs%debug) then

      if (cs%visc_drag) &

        call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, g%HI, &

                      unscale=gv%H_to_mks*us%s_to_T, scalar_pair=.true.)

      call hchksum(mass, 'MEKE mass',g%HI,haloshift=1, unscale=us%RZ_to_kg_m2)

      call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', g%HI, unscale=gv%H_to_mks*us%s_to_T)

      call hchksum(bottomfac2, 'MEKE bottomFac2', g%HI)

      call hchksum(barotrfac2, 'MEKE barotrFac2', g%HI)

      call hchksum(lmixscale, 'MEKE LmixScale', g%HI, unscale=us%L_to_m)

    endif

    if (allocated(meke%Le)) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        meke%Le(i,j) = lmixscale(i,j)

      enddo ; enddo

    endif

    ! Aggregate sources of MEKE (background, frictional and GM)

    !$OMP parallel do default(shared)

    do j=js,je ; do i=is,ie

      src(i,j) = cs%MEKE_BGsrc

    enddo ; enddo

    ! Initialize diagnostics

    if (cs%id_src_adv > 0) src_adv(is:ie, js:je) = 0.

    if (cs%id_src_GM > 0) src_gm(is:ie, js:je) = 0.

    if (cs%id_src_mom_lp > 0) src_mom_lp(is:ie, js:je) = 0.

    if (cs%id_src_mom_bh > 0) src_mom_bh(is:ie, js:je) = 0.

    if (cs%id_src_mom_K4 > 0) src_mom_k4(is:ie, js:je) = 0.

    if (cs%id_src_btm_drag > 0) src_btm_drag(is:ie, js:je) = 0.

    ! Identify any damped diagnostics in first stage of Strang splitting

    any_damping_diags_s1 = any([ &

        cs%id_src_GM > 0, &

        cs%id_src_mom_lp > 0, &

        cs%id_src_mom_bh > 0, &

        cs%id_src_btm_drag > 0 &

    ])

    ! Identify any damped diagnostics

    any_damping_diags = any([ &

        any_damping_diags_s1, &

        cs%id_src_adv > 0, &

        cs%id_src_mom_K4 > 0 &

    ])

    if (cs%MEKE_FrCoeff > 0.) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        src(i,j) = src(i,j) - cs%MEKE_FrCoeff * i_mass(i,j) * meke%mom_src(i,j)

      enddo ; enddo

    endif

    if (allocated(meke%mom_src_bh)) then

      if (cs%MEKE_bhFrCoeff > 0. .and. cs%MEKE_FrCoeff > 0.) then

        bh_coeff = cs%MEKE_bhFrCoeff - cs%MEKE_FrCoeff

      else

        bh_coeff = cs%MEKE_bhFrCoeff

      endif

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        src(i,j) = src(i,j) - bh_coeff * i_mass(i,j) * meke%mom_src_bh(i,j)

      enddo ; enddo

      if (cs%id_src_mom_lp > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_mom_lp(i,j) = -cs%MEKE_FrCoeff * i_mass(i,j) &

              * (meke%mom_src(i,j) - meke%mom_src_bh(i,j))

        enddo ; enddo

      endif

      if (cs%id_src_mom_bh > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_mom_bh(i,j) = -cs%MEKE_bhFrCoeff * i_mass(i,j) * meke%mom_src_bh(i,j)

        enddo ; enddo

      endif

    endif

    if (allocated(meke%GME_snk)) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        src(i,j) = src(i,j) - cs%MEKE_GMECoeff*i_mass(i,j)*meke%GME_snk(i,j)

      enddo ; enddo

    endif

    if (allocated(meke%GM_src)) then

      if (cs%GM_src_alt) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src(i,j) = src(i,j) - cs%MEKE_GMcoeff*meke%GM_src(i,j) / &

                     (gv%H_to_RZ * max(cs%MEKE_min_depth_tot, depth_tot(i,j)))

        enddo ; enddo

      else

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src(i,j) = src(i,j) - cs%MEKE_GMcoeff*i_mass(i,j)*meke%GM_src(i,j)

        enddo ; enddo

        do j=js,je ; do i=is,ie

          src_gm(i,j) = -cs%MEKE_GMcoeff*i_mass(i,j)*meke%GM_src(i,j)

        enddo ; enddo

      endif

    endif

    if (cs%MEKE_equilibrium_restoring) then

      call meke_equilibrium_restoring(cs, g, gv, us, sn_u, sn_v, depth_tot, &

                                      equilibrium_value)

      do j=js,je ; do i=is,ie

        src(i,j) = src(i,j) - cs%MEKE_restoring_rate*(meke%MEKE(i,j) - equilibrium_value(i,j))

      enddo ; enddo

    endif

    if (cs%debug) then

      call hchksum(src, "MEKE src", g%HI, haloshift=0, unscale=us%L_to_m**2*us%s_to_T**3)

    endif

    ! Increase EKE by a full time-steps worth of source

    !$OMP parallel do default(shared)

    do j=js,je ; do i=is,ie

      meke_current(i,j) = meke%MEKE(i,j)

      meke%MEKE(i,j) = (meke%MEKE(i,j) + sdt*src(i,j))*g%mask2dT(i,j)

    enddo ; enddo

    if (use_drag_rate) then

      ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        drag_rate(i,j) = (gv%H_to_RZ * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &

                 cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )

      enddo ; enddo

    else

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        drag_rate(i,j) = 0.

      enddo ; enddo

    endif

    ! First stage of Strang splitting

    !$OMP parallel do default(shared)

    do j=js,je ; do i=is,ie

      damp_rate(i,j) = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)

      if (meke%MEKE(i,j) < 0.) damp_rate(i,j) = 0.

      ! notice that the above line ensures a damping only if MEKE is positive,

      ! while leaving MEKE unchanged if it is negative

    enddo ; enddo

    ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the existing

    !   bit-reproducible solution.

    !$OMP parallel do default(shared)

    do j=js,je ; do i=is,ie

      meke%MEKE(i,j) =  meke%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))

    enddo ; enddo

    if (any_damping_diags_s1) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))

      enddo ; enddo

      if (cs%id_decay > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          meke_decay(i,j) = damp_rate(i,j) * g%mask2dT(i,j)

        enddo ; enddo

      endif

      if (cs%id_src_GM > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_gm(i,j) = src_gm(i,j) * damping(i,j)

        enddo ; enddo

      endif

      if (cs%id_src_mom_lp > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)

        enddo ; enddo

      endif

      if (cs%id_src_mom_bh > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)

        enddo ; enddo

      endif

      if (cs%id_src_btm_drag > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_btm_drag(i,j) = -meke_current(i,j) * ( &

              damp_step * (damp_rate(i,j) * damping(i,j)) &

          )

        enddo ; enddo

        ! Store the effective damping rate if sdt is split

        if (cs%MEKE_KH >= 0. .or. cs%MEKE_K4 >= 0.) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            damp_rate_s1(i,j) = damp_rate(i,j) * damping(i,j)

          enddo ; enddo

        endif

      endif

    endif

    if (cs%kh_flux_enabled .or. cs%MEKE_K4 >= 0.0) then

      ! Update MEKE in the halos for lateral or bi-harmonic diffusion

      call cpu_clock_begin(cs%id_clock_pass)

      call do_group_pass(cs%pass_MEKE, g%Domain)

      call cpu_clock_end(cs%id_clock_pass)

    endif

    if (cs%MEKE_K4 >= 0.0) then

      ! Calculate Laplacian of MEKE using MEKE_uflux and MEKE_vflux as temporary work space.

      !$OMP parallel do default(shared)

      do j=js-1,je+1 ; do i=is-2,ie+1

        ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].

        meke_uflux(i,j) = (g%dy_Cu(i,j)*g%IdxCu_OBCmask(i,j)) * &

            (meke%MEKE(i+1,j) - meke%MEKE(i,j))

      ! This would have units of [R Z L2 T-2 ~> kg s-2]

      ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &

      !     ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &

      !     (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))

      enddo ; enddo

      !$OMP parallel do default(shared)

      do j=js-2,je+1 ; do i=is-1,ie+1

        ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].

        meke_vflux(i,j) = (g%dx_Cv(i,j)*g%IdyCv_OBCmask(i,j)) * &

            (meke%MEKE(i,j+1) - meke%MEKE(i,j))

      ! This would have units of [R Z L2 T-2 ~> kg s-2]

      ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &

      !     ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &

      !     (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))

      enddo ; enddo

      !$OMP parallel do default(shared)

      do j=js-1,je+1 ; do i=is-1,ie+1 ! del2MEKE has units [T-2 ~> s-2].

        del2meke(i,j) = g%IareaT(i,j) * &

            ((meke_uflux(i,j) - meke_uflux(i-1,j)) + (meke_vflux(i,j) - meke_vflux(i,j-1)))

      enddo ; enddo

      ! Bi-harmonic diffusion of MEKE

      !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)

      do j=js,je ; do i=is-1,ie

        k4_here = cs%MEKE_K4 ! [L4 T-1 ~> m4 s-1]

        ! Limit Kh to avoid CFL violations.

        inv_k4_max = 64.0 * sdt * ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &

                     max(g%IareaT(i,j), g%IareaT(i+1,j)))**2

        if (k4_here*inv_k4_max > 0.3) k4_here = 0.3 / inv_k4_max

        ! Here the units of MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].

        meke_uflux(i,j) = ((k4_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &

            ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &

            (del2meke(i+1,j) - del2meke(i,j))

      enddo ; enddo

      !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)

      do j=js-1,je ; do i=is,ie

        k4_here = cs%MEKE_K4 ! [L4 T-1 ~> m4 s-1]

        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

        if (k4_here*inv_k4_max > 0.3) k4_here = 0.3 / inv_k4_max

        ! Here the units of MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].

        meke_vflux(i,j) = ((k4_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &

            ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &

            (del2meke(i,j+1) - del2meke(i,j))

      enddo ; enddo

      ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2 ~> m2 s-2].

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        del4meke(i,j) = (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &

            ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &

             (meke_vflux(i,j-1) - meke_vflux(i,j)))

        src_mom_k4(i,j) = (g%IareaT(i,j)*i_mass(i,j))  * &

            ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &

             (meke_vflux(i,j-1) - meke_vflux(i,j)))

      enddo ; enddo

    endif !

    if (cs%kh_flux_enabled) then

      ! Lateral diffusion of MEKE

      kh_here = max(0., cs%MEKE_Kh)

      !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)

      do j=js,je ; do i=is-1,ie

        ! Limit Kh to avoid CFL violations.

        if (allocated(meke%Kh)) &

          kh_here = max(0., cs%MEKE_Kh) + &

              cs%KhMEKE_Fac*0.5*(meke%Kh(i,j)+meke%Kh(i+1,j))

        if (allocated(meke%Kh_diff)) &

          kh_here = max(0.,cs%MEKE_Kh) + &

              cs%KhMEKE_Fac*0.5*(meke%Kh_diff(i,j)+meke%Kh_diff(i+1,j))

        inv_kh_max = 2.0*sdt * ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &

                     max(g%IareaT(i,j),g%IareaT(i+1,j)))

        if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max

        kh_u(i,j) = kh_here

        ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].

        meke_uflux(i,j) = ((kh_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &

            ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &

            (meke%MEKE(i,j) - meke%MEKE(i+1,j))

      enddo ; enddo

      !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)

      do j=js-1,je ; do i=is,ie

        if (allocated(meke%Kh)) &

          kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac * 0.5*(meke%Kh(i,j)+meke%Kh(i,j+1))

        if (allocated(meke%Kh_diff)) &

          kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac * 0.5*(meke%Kh_diff(i,j)+meke%Kh_diff(i,j+1))

        inv_kh_max = 2.0*sdt * ((g%dx_Cv(i,j)*g%IdyCv(i,j)) * max(g%IareaT(i,j),g%IareaT(i,j+1)))

        if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max

        kh_v(i,j) = kh_here

        ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].

        meke_vflux(i,j) = ((kh_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &

            ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &

            (meke%MEKE(i,j) - meke%MEKE(i,j+1))

      enddo ; enddo

      if (cs%MEKE_advection_factor>0.) then

        advfac = cs%MEKE_advection_factor / sdt ! [T-1 ~> s-1]

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is-1,ie

          ! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].

          if (barohu(i,j)>0.) then

            meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i,j)*advfac

          elseif (barohu(i,j)<0.) then

            meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i+1,j)*advfac

          endif

        enddo ; enddo

        !$OMP parallel do default(shared)

        do j=js-1,je ; do i=is,ie

          ! Here the units of the quantities added to MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].

          if (barohv(i,j)>0.) then

            meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j)*advfac

          elseif (barohv(i,j)<0.) then

            meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j+1)*advfac

          endif

        enddo ; enddo

      endif

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        meke%MEKE(i,j) = meke%MEKE(i,j) + (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &

            ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &

             (meke_vflux(i,j-1) - meke_vflux(i,j)))

      enddo ; enddo

      if (cs%id_src_adv > 0) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          src_adv(i,j) = (g%IareaT(i,j)*i_mass(i,j)) * &

              ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &

               (meke_vflux(i,j-1) - meke_vflux(i,j)))

        enddo ; enddo

      endif

    endif ! MEKE_KH>0

    ! Add on bi-harmonic tendency

    if (cs%MEKE_K4 >= 0.0) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        meke%MEKE(i,j) = meke%MEKE(i,j) + del4meke(i,j)

      enddo ; enddo

    endif

    ! Second stage of Strang splitting

    if (cs%MEKE_KH >= 0.0 .or. cs%MEKE_K4 >= 0.0) then

      ! Recalculate the drag rate, since MEKE has changed.

      if (use_drag_rate) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          drag_rate(i,j) = (gv%H_to_RZ * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &

                 cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )

        enddo ; enddo

      endif

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        damp_rate(i,j) = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)

        if (meke%MEKE(i,j) < 0.) damp_rate(i,j) = 0.

        ! notice that the above line ensures a damping only if MEKE is positive,

        ! while leaving MEKE unchanged if it is negative

      enddo ; enddo

      ! NOTE: MEKE%MEKE cannot use `damping` since we must preserve the

      !   existing bit-reproducible solution.

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        meke%MEKE(i,j) =  meke%MEKE(i,j) / (1. + sdt_damp * damp_rate(i,j))

      enddo ; enddo

      if (any_damping_diags) then

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          damping(i,j) = 1. / (1. + sdt_damp * damp_rate(i,j))

        enddo ; enddo

        if (cs%id_decay > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            meke_decay(i,j) = damp_rate(i,j) * g%mask2dT(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_GM > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_gm(i,j) = src_gm(i,j) * damping(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_mom_lp > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_mom_lp(i,j) = src_mom_lp(i,j) * damping(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_mom_bh > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_mom_bh(i,j) = src_mom_bh(i,j) * damping(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_adv > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_adv(i,j) = src_adv(i,j) * damping(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_mom_K4 > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_mom_k4(i,j) = src_mom_k4(i,j) * damping(i,j)

          enddo ; enddo

        endif

        if (cs%id_src_btm_drag > 0) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            src_btm_drag(i,j) = -meke_current(i,j) * (damp_step &

                * ((damp_rate(i,j) + damp_rate_s1(i,j)) * damping(i,j)) &

            )

          enddo ; enddo

        endif

      endif

    endif ! MEKE_KH>=0

    if (cs%debug) then

      call hchksum(meke%MEKE, "MEKE post-update MEKE", g%HI, haloshift=0, unscale=us%L_T_to_m_s**2)

    endif

  case(eke_file)

    call time_interp_external(cs%eke_handle, time, data_eke, scale=us%m_s_to_L_T**2)

    do j=js,je ; do i=is,ie

      meke%MEKE(i,j) = data_eke(i,j) * g%mask2dT(i,j)

    enddo ; enddo

    call meke_lengthscales(cs, meke, g, gv, us, sn_u, sn_v, meke%MEKE, depth_tot, bottomfac2, barotrfac2, lmixscale)

  case(eke_dbclient)

    call pass_vector(u, v, g%Domain)

    call meke_lengthscales(cs, meke, g, gv, us, sn_u, sn_v, meke%MEKE, depth_tot, bottomfac2, barotrfac2, lmixscale)

    call ml_meke_calculate_features(g, gv, us, cs, meke%Rd_dx_h, u, v, tv, h, dt, features_array)

    call predict_meke(g, us, cs, SIZE(h), time, features_array, meke%MEKE)

  case default

    call mom_error(fatal,"Invalid method specified for calculating EKE")

  end select

  if (cs%MEKE_positive) then

    !$OMP parallel do default(shared)

    do j=js,je ; do i=is,ie

      meke%MEKE(i,j) = max(0., meke%MEKE(i,j))

    enddo ; enddo

  endif

  call cpu_clock_begin(cs%id_clock_pass)

  call do_group_pass(cs%pass_MEKE, g%Domain)

  call cpu_clock_end(cs%id_clock_pass)

  ! Calculate diffusivity for main model to use

  if (cs%MEKE_KhCoeff>0.) then

    if (.not.cs%MEKE_GEOMETRIC) then

      if (cs%use_old_lscale) then

        if (cs%Rd_as_max_scale) then

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            meke%Kh(i,j) = (cs%MEKE_KhCoeff * &

                       sqrt(2.*max(0.,barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j)) ) * &

                       min(meke%Rd_dx_h(i,j), 1.0)

          enddo ; enddo

        else

          !$OMP parallel do default(shared)

          do j=js,je ; do i=is,ie

            meke%Kh(i,j) = cs%MEKE_KhCoeff * &

                sqrt(2.*max(0., barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j))

          enddo ; enddo

        endif

      else

        !$OMP parallel do default(shared)

        do j=js,je ; do i=is,ie

          meke%Kh(i,j) = cs%MEKE_KhCoeff * &

              sqrt(2.*max(0., barotrfac2(i,j)*meke%MEKE(i,j))) * lmixscale(i,j)

        enddo ; enddo

      endif

    endif

  endif

  ! Calculate viscosity for the main model to use

  if (cs%viscosity_coeff_Ku /=0.) then

    do j=js,je ; do i=is,ie

      meke%Ku(i,j) = cs%viscosity_coeff_Ku * sqrt(2.*max(0.,meke%MEKE(i,j))) * lmixscale(i,j)

    enddo ; enddo

  endif

  if (cs%viscosity_coeff_Au /=0.) then

    do j=js,je ; do i=is,ie

      meke%Au(i,j) = cs%viscosity_coeff_Au * sqrt(2.*max(0.,meke%MEKE(i,j))) * lmixscale(i,j)**3

    enddo ; enddo

  endif

  if (allocated(meke%Kh) .or. allocated(meke%Ku) .or. allocated(meke%Au) &

      .or. allocated(meke%Le)) then

    call cpu_clock_begin(cs%id_clock_pass)

    call do_group_pass(cs%pass_Kh, g%Domain)

    call cpu_clock_end(cs%id_clock_pass)

  endif

  ! Offer fields for averaging.

  if (any([cs%id_Ue, cs%id_Ub, cs%id_Ut] > 0)) &

    tmp(:,:) = 0.

  if (cs%id_MEKE>0) call post_data(cs%id_MEKE, meke%MEKE, cs%diag)

  if (cs%id_Ue>0) then

    do j=js,je ; do i=is,ie

      tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j)))

    enddo ; enddo

    call post_data(cs%id_Ue, tmp, cs%diag)

  endif

  if (cs%id_Ub>0) then

    do j=js,je ; do i=is,ie

      tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j) * bottomfac2(i,j)))

    enddo ; enddo

    call post_data(cs%id_Ub, tmp, cs%diag)

  endif

  if (cs%id_Ut>0) then

    do j=js,je ; do i=is,ie

      tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j) * barotrfac2(i,j)))

    enddo ; enddo

    call post_data(cs%id_Ut, tmp, cs%diag)

  endif

  if (cs%id_Kh>0) call post_data(cs%id_Kh, meke%Kh, cs%diag)

  if (cs%id_Ku>0) call post_data(cs%id_Ku, meke%Ku, cs%diag)

  if (cs%id_Au>0) call post_data(cs%id_Au, meke%Au, cs%diag)

  if (cs%id_KhMEKE_u>0) call post_data(cs%id_KhMEKE_u, kh_u, cs%diag)

  if (cs%id_KhMEKE_v>0) call post_data(cs%id_KhMEKE_v, kh_v, cs%diag)

  if (cs%id_src>0) call post_data(cs%id_src, src, cs%diag)

  if (cs%id_src_adv>0) call post_data(cs%id_src_adv, src_adv, cs%diag)

  if (cs%id_src_mom_K4>0) call post_data(cs%id_src_mom_K4, src_mom_k4, cs%diag)

  if (cs%id_src_btm_drag>0) call post_data(cs%id_src_btm_drag, src_btm_drag, cs%diag)

  if (cs%id_src_GM>0) call post_data(cs%id_src_GM, src_gm, cs%diag)

  if (cs%id_src_mom_lp>0) call post_data(cs%id_src_mom_lp, src_mom_lp, cs%diag)

  if (cs%id_src_mom_bh>0) call post_data(cs%id_src_mom_bh, src_mom_bh, cs%diag)

  if (cs%id_decay>0) call post_data(cs%id_decay, meke_decay, cs%diag)

  if (cs%id_GM_src>0) call post_data(cs%id_GM_src, meke%GM_src, cs%diag)

  if (cs%id_mom_src>0) call post_data(cs%id_mom_src, meke%mom_src, cs%diag)

  if (cs%id_mom_src_bh>0) call post_data(cs%id_mom_src_bh, meke%mom_src_bh, cs%diag)

  if (cs%id_GME_snk>0) call post_data(cs%id_GME_snk, meke%GME_snk, cs%diag)

  if (cs%id_Le>0) call post_data(cs%id_Le, lmixscale, cs%diag)

  if (cs%id_gamma_b>0) then

    do j=js,je ; do i=is,ie

      bottomfac2(i,j) = sqrt(bottomfac2(i,j))

    enddo ; enddo

    call post_data(cs%id_gamma_b, bottomfac2, cs%diag)

  endif

  if (cs%id_gamma_t>0) then

    do j=js,je ; do i=is,ie

      barotrfac2(i,j) = sqrt(barotrfac2(i,j))

    enddo ; enddo

    call post_data(cs%id_gamma_t, barotrfac2, cs%diag)

  endif

end subroutine step_forward_meke

!> Calculates the equilibrium solution where the source depends only on MEKE diffusivity

!! and there is no lateral diffusion of MEKE.

!! Results is in MEKE%MEKE.

subroutine meke_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)

  type(ocean_grid_type),             intent(inout) :: g    !< Ocean grid.

  type(verticalgrid_type),           intent(in)    :: GV   !< Ocean vertical grid structure.

  type(unit_scale_type),             intent(in)    :: US   !< A dimensional unit scaling type

  type(meke_cs),                     intent(in)    :: CS   !< MEKE control structure.

  type(meke_type),                   intent(inout) :: MEKE !< MEKE fields

  real, dimension(SZIB_(G),SZJ_(G)), intent(in)    :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJB_(G)), intent(in)    :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: drag_rate_visc !< Mean flow velocity contribution

                                                           !! to the MEKE drag rate [H T-1 ~> m s-1 or kg m-2 s-1]

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: I_mass  !< Inverse of column mass [R-1 Z-1 ~> m2 kg-1].

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].

  ! Local variables

  real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]

  real :: SN   ! The local Eady growth rate [T-1 ~> s-1]

  real :: bottomFac2, barotrFac2    ! Vertical structure factors [nondim]

  real :: LmixScale, LRhines, LEady ! Various mixing length scales [L ~> m]

  real :: KhCoeff ! A copy of MEKE_KhCoeff from the control structure [nondim]

  real :: Kh    ! A lateral diffusivity [L2 T-1 ~> m2 s-1]

  real :: Ubg2  ! Background (tidal?) velocity squared [L2 T-2 ~> m2 s-2]

  real :: cd2   ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]

  real :: drag_rate ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].

  real :: src   ! The sum of MEKE sources [L2 T-3 ~> W kg-1]

  real :: ldamping  ! The MEKE damping rate [T-1 ~> s-1].

  real :: EKE, EKEmin, EKEmax, EKEerr ! [L2 T-2 ~> m2 s-2]

  real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1]

  real :: FatH    ! Coriolis parameter at h points, used to compute topographic beta [T-1 ~> s-1]

  real :: beta_topo_x, beta_topo_y    ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]

  real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]

  integer :: i, j, is, ie, js, je, n1, n2

  real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2].

  logical :: useSecant, debugIteration

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  debugiteration = .false.

  khcoeff = cs%MEKE_KhCoeff

  ubg2 = cs%MEKE_Uscale**2

  cd2 = cs%cdrag**2

  tolerance = 1.0e-12*us%m_s_to_L_T**2

  h_neglect = gv%H_subroundoff

!$OMP do

  do j=js,je ; do i=is,ie

    ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)

    ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v

    sn = min(sn_u(i,j), sn_u(i-1,j), sn_v(i,j), sn_v(i,j-1))

    if (cs%MEKE_equilibrium_alt) then

      meke%MEKE(i,j) = (cs%MEKE_GEOMETRIC_alpha * sn * depth_tot(i,j))**2 / cd2

    else

      fath = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &

                   (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1))) ! Coriolis parameter at h points

      ! Since zero-bathymetry cells are masked, this avoids calculations on land

      if (cs%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then

        beta_topo_x = 0. ; beta_topo_y = 0.

      else

        !### Consider different combinations of these estimates of topographic beta.

        beta_topo_x = -cs%MEKE_topographic_beta * fath * 0.5 * ( &

                      (depth_tot(i+1,j)-depth_tot(i,j)) * g%IdxCu(i,j)  &

                  / max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &

              +       (depth_tot(i,j)-depth_tot(i-1,j)) * g%IdxCu(i-1,j) &

                  / max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )

        beta_topo_y = -cs%MEKE_topographic_beta * fath * 0.5 * ( &

                      (depth_tot(i,j+1)-depth_tot(i,j)) * g%IdyCv(i,j)  &

                  / max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &

                      (depth_tot(i,j)-depth_tot(i,j-1)) * g%IdyCv(i,j-1) &

                  / max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )

      endif

      beta =  sqrt(((g%dF_dx(i,j) + beta_topo_x)**2) + &

                   ((g%dF_dy(i,j) + beta_topo_y)**2) )

      if (khcoeff*sn*i_mass(i,j)>0.) then

        ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E

        ekemin = 0.   ! Use the trivial root as the left bracket

        resmin = 0.   ! Need to detect direction of left residual

        ekemax = 0.01*us%m_s_to_L_T**2 ! First guess at right bracket

        usesecant = .false. ! Start using a bisection method

        ! First find right bracket for which resid<0

        resid = 1.0*us%m_to_L**2*us%T_to_s**3 ; n1 = 0

        do while (resid>0.)

          n1 = n1 + 1

          eke = ekemax

          call meke_lengthscales_0d(cs, us, g%areaT(i,j), beta, depth_tot(i,j), &

                                    meke%Rd_dx_h(i,j), sn, eke, &

                                    bottomfac2, barotrfac2, lmixscale, lrhines, leady)

          ! TODO: Should include resolution function in Kh

          kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)

          src = kh * (sn * sn)

          drag_rate = (gv%H_to_RZ * i_mass(i,j)) * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )

          ldamping = cs%MEKE_damping + drag_rate * bottomfac2

          resid = src - ldamping * eke

          ! if (debugIteration) then

          !   write(0,*) n1, 'EKE=',EKE,'resid=',resid

          !   write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin

          !   write(0,*) 'src=',src,'ldamping=',ldamping

          !   write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2

          !   write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2

          ! endif

          if (resid>0.) then    ! EKE is to the left of the root

            ekemin = eke        ! so we move the left bracket here

            ekemax = 10. * eke  ! and guess again for the right bracket

            if (resid<resmin) usesecant = .true.

            resmin = resid

            if (ekemax > 2.e17*us%m_s_to_L_T**2) then

              if (debugiteration) stop 'Something has gone very wrong'

              debugiteration = .true.

              resid = 1. ; n1 = 0

              ekemin = 0. ; resmin = 0.

              ekemax = 0.01*us%m_s_to_L_T**2

              usesecant = .false.

            endif

          endif

        enddo ! while(resid>0.) searching for right bracket

        resmax = resid

        ! Bisect the bracket

        n2 = 0 ; ekeerr = ekemax - ekemin

        do while (ekeerr > tolerance)

          n2 = n2 + 1

          if (usesecant) then

            eke = ekemin + (ekemax - ekemin) * (resmin / (resmin - resmax))

          else

            eke = 0.5 * (ekemin + ekemax)

          endif

          ekeerr = min( eke-ekemin, ekemax-eke )

          ! TODO: Should include resolution function in Kh

          kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)

          src = kh * (sn * sn)

          drag_rate = (gv%H_to_RZ * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )

          ldamping = cs%MEKE_damping + drag_rate * bottomfac2

          resid = src - ldamping * eke

          if (usesecant .and. resid>resmin) usesecant = .false.

          if (resid>0.) then              ! EKE is to the left of the root

            ekemin = eke                  ! so we move the left bracket here

            if (resid<resmin) usesecant = .true.

            resmin = resid                ! Save this for the secant method

          elseif (resid<0.) then          ! EKE is to the right of the root

            ekemax = eke                  ! so we move the right bracket here

            resmax = resid                ! Save this for the secant method

          else

            exit                          ! resid=0 => EKE is exactly at the root

          endif

          if (n2>200) stop 'Failing to converge?'

        enddo ! while(EKEmax-EKEmin>tolerance)

      else

        eke = 0.

      endif

      meke%MEKE(i,j) = eke

    endif

  enddo ; enddo

end subroutine meke_equilibrium

!< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into

!! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value

subroutine meke_equilibrium_restoring(CS, G, GV, US, SN_u, SN_v, depth_tot, &

                                      equilibrium_value)

  type(ocean_grid_type),             intent(inout) :: G    !< Ocean grid.

  type(verticalgrid_type),           intent(in)    :: GV   !< Ocean vertical grid structure.

  type(unit_scale_type),             intent(in)    :: US   !< A dimensional unit scaling type.

  type(meke_cs),                     intent(in)    :: CS   !< MEKE control structure.

  real, dimension(SZIB_(G),SZJ_(G)), intent(in)    :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJB_(G)), intent(in)    :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].

  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: equilibrium_value

      !< Equilbrium value of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2]

  ! Local variables

  real :: SN                      ! The local Eady growth rate [T-1 ~> s-1]

  integer :: i, j, is, ie, js, je ! local indices

  real :: cd2                     ! The square of the drag coefficient [nondim]

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  cd2 = cs%cdrag**2

  equilibrium_value(:,:) = 0.0

!$OMP do

  do j=js,je ; do i=is,ie

    ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)

    ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v

    sn = min(sn_u(i,j), sn_u(i-1,j), sn_v(i,j), sn_v(i,j-1))

    equilibrium_value(i,j) = (cs%MEKE_GEOMETRIC_alpha * sn * depth_tot(i,j))**2 / cd2

  enddo ; enddo

  if (cs%id_MEKE_equilibrium>0) call post_data(cs%id_MEKE_equilibrium, equilibrium_value, cs%diag)

end subroutine meke_equilibrium_restoring

!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$

!! functions that are ratios of either bottom or barotropic eddy energy to the

!! column eddy energy, respectively.  See \ref section_MEKE_equations.

subroutine meke_lengthscales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, &

                             bottomFac2, barotrFac2, LmixScale)

  type(meke_cs),                     intent(in)    :: CS   !< MEKE control structure.

  type(meke_type),                   intent(in)    :: MEKE !< MEKE field

  type(ocean_grid_type),             intent(inout) :: G    !< Ocean grid.

  type(verticalgrid_type),           intent(in)    :: GV   !< Ocean vertical grid structure.

  type(unit_scale_type),             intent(in)    :: US   !< A dimensional unit scaling type

  real, dimension(SZIB_(G),SZJ_(G)), intent(in)    :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJB_(G)), intent(in)    :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: EKE  !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].

  real, dimension(SZI_(G),SZJ_(G)),  intent(in)    :: depth_tot !< The thickness of the water column [H ~> m or kg m-2].

  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: bottomFac2 !< gamma_b^2 [nondim]

  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: barotrFac2 !< gamma_t^2 [nondim]

  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: LmixScale !< Eddy mixing length [L ~> m].

  ! Local variables

  real, dimension(SZI_(G),SZJ_(G)) :: LRhines, LEady  ! Possible mixing length scales [L ~> m]

  real :: beta ! Combined topographic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]

  real :: SN   ! The local Eady growth rate [T-1 ~> s-1]

  real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1]

  real :: beta_topo_x, beta_topo_y  ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]

  real :: h_neglect ! A negligible thickness [H ~> m or kg m-2]

  integer :: i, j, is, ie, js, je

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  h_neglect = gv%H_subroundoff

!$OMP do

  do j=js,je ; do i=is,ie

    if (.not.cs%use_old_lscale) then

      if (cs%aEady > 0.) then

        sn = 0.25 * ( (sn_u(i,j) + sn_u(i-1,j)) + (sn_v(i,j) + sn_v(i,j-1)) )

      else

        sn = 0.

      endif

      fath = 0.25* ( ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1) ) + &

                     ( g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1) ) )  ! Coriolis parameter at h points

      ! If depth_tot is zero, then a division by zero FPE will be raised.  In this

      ! case, we apply Adcroft's rule of reciprocals and set the term to zero.

      ! Since zero-bathymetry cells are masked, this should not affect values.

      if (cs%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then

        beta_topo_x = 0. ; beta_topo_y = 0.

      else

        !### Consider different combinations of these estimates of topographic beta.

        beta_topo_x = -cs%MEKE_topographic_beta * fath * 0.5 * ( &

                      (depth_tot(i+1,j)-depth_tot(i,j)) * g%IdxCu(i,j)  &

                 / max(depth_tot(i+1,j), depth_tot(i,j), h_neglect) &

              +       (depth_tot(i,j)-depth_tot(i-1,j)) * g%IdxCu(i-1,j) &

                 / max(depth_tot(i,j), depth_tot(i-1,j), h_neglect) )

        beta_topo_y = -cs%MEKE_topographic_beta * fath * 0.5 * ( &

                      (depth_tot(i,j+1)-depth_tot(i,j)) * g%IdyCv(i,j)  &

                 / max(depth_tot(i,j+1), depth_tot(i,j), h_neglect) + &

                      (depth_tot(i,j)-depth_tot(i,j-1)) * g%IdyCv(i,j-1) &

                 / max(depth_tot(i,j), depth_tot(i,j-1), h_neglect) )

      endif

      beta =  sqrt(((g%dF_dx(i,j) + beta_topo_x)**2) + &

                   ((g%dF_dy(i,j) + beta_topo_y)**2) )

    else

      beta = 0.

    endif

    ! Returns bottomFac2, barotrFac2 and LmixScale

    call meke_lengthscales_0d(cs, us, g%areaT(i,j), beta, depth_tot(i,j),  &

                              meke%Rd_dx_h(i,j), sn, meke%MEKE(i,j), &

                              bottomfac2(i,j), barotrfac2(i,j), lmixscale(i,j), &

                              lrhines(i,j), leady(i,j))

  enddo ; enddo

  if (cs%id_Lrhines>0) call post_data(cs%id_LRhines, lrhines, cs%diag)

  if (cs%id_Leady>0) call post_data(cs%id_LEady, leady, cs%diag)

end subroutine meke_lengthscales

!> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$

!! functions that are ratios of either bottom or barotropic eddy energy to the

!! column eddy energy, respectively.  See \ref section_MEKE_equations.

subroutine meke_lengthscales_0d(CS, US, area, beta, depth_tot, Rd_dx, SN, EKE, &

                                bottomFac2, barotrFac2, LmixScale, Lrhines, Leady)

  type(meke_cs), intent(in)    :: CS         !< MEKE control structure.

  type(unit_scale_type), intent(in) :: US    !< A dimensional unit scaling type

  real,          intent(in)    :: area       !< Grid cell area [L2 ~> m2]

  real,          intent(in)    :: beta       !< Planetary beta = \f$ \nabla f\f$  [T-1 L-1 ~> s-1 m-1]

  real,          intent(in)    :: depth_tot  !< The total thickness of the water column [H ~> m or kg m-2]

  real,          intent(in)    :: Rd_dx      !< Resolution Ld/dx [nondim].

  real,          intent(in)    :: SN         !< Eady growth rate [T-1 ~> s-1].

  real,          intent(in)    :: EKE        !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].

  real,          intent(out)   :: bottomFac2 !< gamma_b^2 [nondim]

  real,          intent(out)   :: barotrFac2 !< gamma_t^2 [nondim]

  real,          intent(out)   :: LmixScale  !< Eddy mixing length [L ~> m].

  real,          intent(out)   :: Lrhines    !< Rhines length scale [L ~> m].

  real,          intent(out)   :: Leady      !< Eady length scale [L ~> m].

  ! Local variables

  real :: Lgrid, Ldeform, Lfrict ! Length scales [L ~> m]

  real :: Ue  ! An eddy velocity [L T-1 ~> m s-1]

  ! Length scale for MEKE derived diffusivity

  lgrid = sqrt(area)               ! Grid scale

  ldeform = lgrid * rd_dx          ! Deformation scale

  lfrict = depth_tot / cs%cdrag    ! Frictional arrest scale

  ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy

  ! used in calculating bottom drag

  bottomfac2 = cs%MEKE_CD_SCALE**2

  if (lfrict*cs%MEKE_Cb>0.) bottomfac2 = bottomfac2 + 1./( 1. + cs%MEKE_Cb*(ldeform/lfrict) )**0.8

  bottomfac2 = max(bottomfac2, cs%MEKE_min_gamma)

  ! gamma_t^2 is the ratio of barotropic eddy energy to mean column eddy energy

  ! used in the velocity scale for diffusivity

  barotrfac2 = 1.

  if (lfrict*cs%MEKE_Ct>0.) barotrfac2 = 1. / ( 1. + cs%MEKE_Ct*(ldeform/lfrict) )**0.25

  barotrfac2 = max(barotrfac2, cs%MEKE_min_gamma)

  if (cs%use_old_lscale) then

    if (cs%Rd_as_max_scale) then

      lmixscale = min(ldeform, lgrid) ! The smaller of Ld or dx

    else

      lmixscale = lgrid

    endif

  else

    ue = sqrt( 2.0 * max( 0., barotrfac2*eke ) ) ! Barotropic eddy flow scale

    lrhines = sqrt( ue / max( beta, 1.e-30*us%T_to_s*us%L_to_m ) )       ! Rhines scale

    if (cs%aEady > 0.) then

      leady = ue / max( sn, 1.e-15*us%T_to_s ) ! Bound Eady time-scale < 1e15 seconds

    else

      leady = 0.

    endif

    if (cs%use_min_lscale) then

      lmixscale = cs%lscale_maxval

      if (cs%aDeform*ldeform > 0.) lmixscale = min(lmixscale,cs%aDeform*ldeform)

      if (cs%aFrict *lfrict  > 0.) lmixscale = min(lmixscale,cs%aFrict *lfrict)

      if (cs%aRhines*lrhines > 0.) lmixscale = min(lmixscale,cs%aRhines*lrhines)

      if (cs%aEady  *leady   > 0.) lmixscale = min(lmixscale,cs%aEady  *leady)

      if (cs%aGrid  *lgrid   > 0.) lmixscale = min(lmixscale,cs%aGrid  *lgrid)

      if (cs%Lfixed          > 0.) lmixscale = min(lmixscale,cs%Lfixed)

    else

      lmixscale = 0.

      if (cs%aDeform*ldeform > 0.) lmixscale = lmixscale + 1./(cs%aDeform*ldeform)

      if (cs%aFrict *lfrict  > 0.) lmixscale = lmixscale + 1./(cs%aFrict *lfrict)

      if (cs%aRhines*lrhines > 0.) lmixscale = lmixscale + 1./(cs%aRhines*lrhines)

      if (cs%aEady  *leady   > 0.) lmixscale = lmixscale + 1./(cs%aEady  *leady)

      if (cs%aGrid  *lgrid   > 0.) lmixscale = lmixscale + 1./(cs%aGrid  *lgrid)

      if (cs%Lfixed          > 0.) lmixscale = lmixscale + 1./cs%Lfixed

      if (lmixscale > 0.) lmixscale = 1. / lmixscale

    endif

  endif

end subroutine meke_lengthscales_0d

!> Initializes the MOM_MEKE module and reads parameters.

!! Returns True if module is to be used, otherwise returns False.

logical function meke_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, MEKE, restart_CS, meke_in_dynamics)

  type(time_type),         intent(in)    :: time       !< The current model time.

  type(ocean_grid_type),   intent(inout) :: g          !< The ocean's grid structure.

  type(verticalgrid_type), intent(in)    :: gv         !< Ocean vertical grid structure.

  type(unit_scale_type),   intent(in)    :: us         !< A dimensional unit scaling type

  type(param_file_type),   intent(in)    :: param_file !< Parameter file parser structure.

  type(dbcomms_cs_type),   intent(in)    :: dbcomms_cs !< Database communications control structure

  type(diag_ctrl), target, intent(inout) :: diag       !< Diagnostics structure.

  type(meke_cs),           intent(inout) :: cs         !< MEKE control structure.

  type(meke_type),         intent(inout) :: meke       !< MEKE fields

  type(mom_restart_cs),    intent(in)    :: restart_cs !< MOM restart control structure

  logical,                 intent(  out) :: meke_in_dynamics !< If true, MEKE is stepped forward in dynamics

                                                             !! otherwise in tracer dynamics

  ! Local variables

  real :: meke_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value [T ~> s]

  real :: cdrag            ! The default bottom drag coefficient [nondim].

  character(len=200) :: eke_filename, eke_varname, inputdir

  character(len=16) :: eke_source_str

  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed

  logical :: laplacian, biharmonic, coldstart

  ! This include declares and sets the variable "version".

# include "version_variable.h"

  character(len=40)  :: mdl = "MOM_MEKE" ! This module's name.

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  ! Determine whether this module will be used

  call get_param(param_file, mdl, "USE_MEKE", meke_init, default=.false., do_not_log=.true.)

  call log_version(param_file, mdl, version, "", all_default=.not.meke_init)

  call get_param(param_file, mdl, "USE_MEKE", meke_init, &

                 "If true, turns on the MEKE scheme which calculates "// &

                 "a sub-grid mesoscale eddy kinetic energy budget.", &

                 default=.false.)

  if (.not. meke_init) return

  cs%initialized = .true.

  call get_param(param_file, mdl, "MEKE_IN_DYNAMICS", meke_in_dynamics, &

                 "If true, step MEKE forward with the dynamics "// &

                 "otherwise with the tracer timestep.", &

                 default=.true.)

  call get_param(param_file, mdl, "EKE_SOURCE", eke_source_str, &

                 "Determine the where EKE comes from:\n" // &

                 "  'prog': Calculated solving EKE equation\n"// &

                 "  'file': Read in from a file\n"            // &

                 "  'dbclient': Retrieved from ML-database", default='prog')

  call mom_mesg("MEKE_init: reading parameters ", 5)

  select case (lowercase(eke_source_str))

  case("file")

    cs%eke_src = eke_file

    call time_interp_external_init

    call get_param(param_file, mdl, "EKE_FILE", eke_filename, &

                 "A file in which to find the eddy kineteic energy variable.", &

                 default="eke_file.nc")

    call get_param(param_file, mdl, "EKE_VARIABLE", eke_varname, &

                 "The name of the eddy kinetic energy variable to read from "//&

                 "EKE_FILE to use in MEKE.", &

                 default="eke")

    call get_param(param_file, mdl, "INPUTDIR", inputdir, &

                 "The directory in which all input files are found.", &

                 default=".", do_not_log=.true.)

    inputdir = slasher(inputdir)

    eke_filename = trim(inputdir) // trim(eke_filename)

    cs%eke_handle = init_external_field(eke_filename, eke_varname, domain=g%Domain%mpp_domain)

  case("prog")

    cs%eke_src = eke_prog

    ! Read all relevant parameters and write them to the model log.

    call get_param(param_file, mdl, "MEKE_DAMPING", cs%MEKE_damping, &

                   "The local depth-independent MEKE dissipation rate.", &

                   units="s-1", default=0.0, scale=us%T_to_s)

    call get_param(param_file, mdl, "MEKE_CD_SCALE", cs%MEKE_Cd_scale, &

                   "The ratio of the bottom eddy velocity to the column mean "//&

                   "eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1 "//&

                   "to account for the surface intensification of MEKE.", &

                   units="nondim", default=0.)

    call get_param(param_file, mdl, "MEKE_CB", cs%MEKE_Cb, &

                   "A coefficient in the expression for the ratio of bottom projected "//&

                   "eddy energy and mean column energy (see Jansen et al. 2015).",&

                   units="nondim", default=25.)

    call get_param(param_file, mdl, "MEKE_MIN_GAMMA2", cs%MEKE_min_gamma, &

                   "The minimum allowed value of gamma_b^2.",&

                   units="nondim", default=0.0001)

    call get_param(param_file, mdl, "MEKE_CT", cs%MEKE_Ct, &

                   "A coefficient in the expression for the ratio of barotropic "//&

                   "eddy energy and mean column energy (see Jansen et al. 2015).",&

                   units="nondim", default=50.)

    call get_param(param_file, mdl, "MEKE_GMCOEFF", cs%MEKE_GMcoeff, &

                   "The efficiency of the conversion of potential energy "//&

                   "into MEKE by the thickness mixing parameterization. "//&

                   "If MEKE_GMCOEFF is negative, this conversion is not "//&

                   "used or calculated.", units="nondim", default=-1.0)

    call get_param(param_file, mdl, "MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &

                   "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//&

                   "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)

    call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &

                   "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&

                   "thickness diffusion.", units="nondim", default=0.05)

    call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", cs%MEKE_equilibrium_alt, &

                   "If true, use an alternative formula for computing the (equilibrium) "//&

                   "initial value of MEKE.", default=.false.)

    call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", cs%MEKE_equilibrium_restoring, &

                   "If true, restore MEKE back to its equilibrium value, which is calculated at "//&

                   "each time step.", default=.false.)

    if (cs%MEKE_equilibrium_restoring) then

      call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", meke_restoring_timescale, &

                     "The timescale used to nudge MEKE toward its equilibrium value.", &

                     units="s", default=1e6, scale=us%s_to_T)

      cs%MEKE_restoring_rate = 1.0 / meke_restoring_timescale

    endif

    call get_param(param_file, mdl, "MEKE_FRCOEFF", cs%MEKE_FrCoeff, &

                   "The efficiency of the conversion of mean energy into "//&

                   "MEKE.  If MEKE_FRCOEFF is negative, this conversion "//&

                   "is not used or calculated.", units="nondim", default=-1.0)

    call get_param(param_file, mdl, "MEKE_BHFRCOEFF", cs%MEKE_bhFrCoeff, &

                 "The efficiency of the conversion of mean energy into "//&

                 "MEKE by the biharmonic dissipation.  If MEKE_bhFRCOEFF is negative, this conversion "//&

                 "is not used or calculated.", units="nondim", default=-1.0)

    call get_param(param_file, mdl, "MEKE_GMECOEFF", cs%MEKE_GMECoeff, &

                   "The efficiency of the conversion of MEKE into mean energy "//&

                   "by GME.  If MEKE_GMECOEFF is negative, this conversion "//&

                   "is not used or calculated.", units="nondim", default=-1.0)

    call get_param(param_file, mdl, "MEKE_BGSRC", cs%MEKE_BGsrc, &

                   "A background energy source for MEKE.", &

                   units="W kg-1", default=0.0, scale=us%m_to_L**2*us%T_to_s**3)

    call get_param(param_file, mdl, "MEKE_KH", cs%MEKE_Kh, &

                   "A background lateral diffusivity of MEKE. "//&

                   "Use a negative value to not apply lateral diffusion to MEKE.", &

                   units="m2 s-1", default=-1.0, scale=us%m_to_L**2*us%T_to_s)

    call get_param(param_file, mdl, "MEKE_K4", cs%MEKE_K4, &

                   "A lateral bi-harmonic diffusivity of MEKE. "//&

                   "Use a negative value to not apply bi-harmonic diffusion to MEKE.", &

                   units="m4 s-1", default=-1.0, scale=us%m_to_L**4*us%T_to_s)

    call get_param(param_file, mdl, "MEKE_DTSCALE", cs%MEKE_dtScale, &

                   "A scaling factor to accelerate the time evolution of MEKE.", &

                   units="nondim", default=1.0)

    call get_param(param_file, mdl, "MEKE_POSITIVE", cs%MEKE_positive, &

                   "If true, it guarantees that MEKE will always be >= 0.", &

                   default=.false.)

  case("dbclient")

    cs%eke_src = eke_dbclient

    call ml_meke_init(diag, g, us, time, param_file, dbcomms_cs, cs)

  case default

    call mom_error(fatal, "Invalid method selected for calculating EKE")

  end select

  ! GMM, make sure all parameters used to calculated MEKE are within the above if

  call get_param(param_file, mdl, "MEKE_KHCOEFF", cs%MEKE_KhCoeff, &

                 "A scaling factor in the expression for eddy diffusivity "//&

                 "which is otherwise proportional to the MEKE velocity- "//&

                 "scale times an eddy mixing-length. This factor "//&

                 "must be >0 for MEKE to contribute to the thickness/ "//&

                 "and tracer diffusivity in the rest of the model.", &

                 units="nondim", default=1.0)

  call get_param(param_file, mdl, "MEKE_USCALE", cs%MEKE_Uscale, &

                 "The background velocity that is combined with MEKE to "//&

                 "calculate the bottom drag.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)

  call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", cs%GM_src_alt, &

                 "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&

                 "than the streamfunction for the MEKE GM source term.", default=.false.)

  call get_param(param_file, mdl, "MEKE_MIN_DEPTH_TOT", cs%MEKE_min_depth_tot, &

                 "The minimum total depth over which to distribute MEKE energy sources.  "//&

                 "When the total depth is less than this, the sources are scaled away.", &

                 units="m", default=1.0, scale=gv%m_to_H, do_not_log=.not.cs%GM_src_alt)

  call get_param(param_file, mdl, "MEKE_VISC_DRAG", cs%visc_drag, &

                 "If true, use the vertvisc_type to calculate the bottom "//&

                 "drag acting on MEKE.", default=.true.)

  call get_param(param_file, mdl, "MEKE_KHTH_FAC", meke%KhTh_fac, &

                 "A factor that maps MEKE%Kh to KhTh.", units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_KHTR_FAC", meke%KhTr_fac, &

                 "A factor that maps MEKE%Kh to KhTr.", units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", cs%KhMEKE_Fac, &

                 "A factor that maps MEKE%Kh to Kh for MEKE itself.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_OLD_LSCALE", cs%use_old_lscale, &

                 "If true, use the old formula for length scale which is "//&

                 "a function of grid spacing and deformation radius.",  &

                 default=.false.)

  call get_param(param_file, mdl, "MEKE_MIN_LSCALE", cs%use_min_lscale, &

                 "If true, use a strict minimum of provided length scales "//&

                 "rather than harmonic mean.",  &

                 default=.false.)

  call get_param(param_file, mdl, "MEKE_LSCALE_MAX_VAL", cs%lscale_maxval, &

                 "The ceiling on the value of the MEKE length scale when MEKE_MIN_LSCALE=True.  "//&

                 "The default is the distance from the equator to the pole on Earth, as "//&

                 "estimated by enlightenment era scientists, but should probably scale with RAD_EARTH.", &

                 units="m", default=1.0e7, scale=us%m_to_L, do_not_log=.not.cs%use_min_lscale)

  call get_param(param_file, mdl, "MEKE_RD_MAX_SCALE", cs%Rd_as_max_scale, &

                 "If true, the length scale used by MEKE is the minimum of "//&

                 "the deformation radius or grid-spacing. Only used if "//&

                 "MEKE_OLD_LSCALE=True", default=.false.)

  call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", cs%viscosity_coeff_Ku, &

                 "If non-zero, is the scaling coefficient in the expression for "//&

                 "viscosity used to parameterize harmonic lateral momentum mixing by "//&

                 "unresolved eddies represented by MEKE. Can be negative to "//&

                 "represent backscatter from the unresolved eddies.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", cs%viscosity_coeff_Au, &

                 "If non-zero, is the scaling coefficient in the expression for "//&

                 "viscosity used to parameterize biharmonic lateral momentum mixing by "//&

                 "unresolved eddies represented by MEKE. Can be negative to "//&

                 "represent backscatter from the unresolved eddies.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", cs%Lfixed, &

                 "If positive, is a fixed length contribution to the expression "//&

                 "for mixing length used in MEKE-derived diffusivity.", &

                 units="m", default=0.0, scale=us%m_to_L)

  call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", cs%fixed_total_depth, &

                 "If true, use the nominal bathymetric depth as the estimate of the "//&

                 "time-varying ocean depth.  Otherwise base the depth on the total ocean mass "//&

                 "per unit area.", default=.true.)

  call get_param(param_file, mdl, "MEKE_TOTAL_DEPTH_RHO", cs%rho_fixed_total_depth, &

                 "A density used to translate the nominal bathymetric depth into an estimate "//&

                 "of the total ocean mass per unit area when MEKE_FIXED_TOTAL_DEPTH is true.", &

                 units="kg m-3", default=gv%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &

                 do_not_log=(gv%Boussinesq.or.(.not.cs%fixed_total_depth)))

  call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", cs%aDeform, &

                 "If positive, is a coefficient weighting the deformation scale "//&

                 "in the expression for mixing length used in MEKE-derived diffusivity.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ALPHA_RHINES", cs%aRhines, &

                 "If positive, is a coefficient weighting the Rhines scale "//&

                 "in the expression for mixing length used in MEKE-derived diffusivity.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ALPHA_EADY", cs%aEady, &

                 "If positive, is a coefficient weighting the Eady length scale "//&

                 "in the expression for mixing length used in MEKE-derived diffusivity.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ALPHA_FRICT", cs%aFrict, &

                 "If positive, is a coefficient weighting the frictional arrest scale "//&

                 "in the expression for mixing length used in MEKE-derived diffusivity.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ALPHA_GRID", cs%aGrid, &

                 "If positive, is a coefficient weighting the grid-spacing as a scale "//&

                 "in the expression for mixing length used in MEKE-derived diffusivity.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_COLD_START", coldstart, &

                 "If true, initialize EKE to zero. Otherwise a local equilibrium solution "//&

                 "is used as an initial condition for EKE.", default=.false.)

  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", meke%backscatter_Ro_c, &

                 "The coefficient in the Rossby number function for scaling the biharmonic "//&

                 "frictional energy source. Setting to non-zero enables the Rossby number function.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", meke%backscatter_Ro_pow, &

                 "The power in the Rossby number function for scaling the biharmonic "//&

                 "frictional energy source.", units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ADVECTION_FACTOR", cs%MEKE_advection_factor, &

                 "A scale factor in front of advection of eddy energy. Zero turns advection off. "//&

                 "Using unity would be normal but other values could accommodate a mismatch "//&

                 "between the advecting barotropic flow and the vertical structure of MEKE.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "MEKE_ADVECTION_BUG", cs%MEKE_advection_bug, &

                 "If true, recover a bug in the calculation of the barotropic transport for "//&

                 "the advection of MEKE.  With the bug, only the transports in the deepest "//&

                 "layer are used.", default=.false., do_not_log=(cs%MEKE_advection_factor<=0.))

  call get_param(param_file, mdl, "MEKE_TOPOGRAPHIC_BETA", cs%MEKE_topographic_beta, &

                 "A scale factor to determine how much topographic beta is weighed in " //&

                 "computing beta in the expression of Rhines scale. Use 1 if full "//&

                 "topographic beta effect is considered; use 0 if it's completely ignored.", &

                 units="nondim", default=0.0)

  call get_param(param_file, mdl, "SQG_USE_MEKE", cs%sqg_use_MEKE, &

                 "If true, the eddy scale of MEKE is used for the SQG vertical structure ",&

                 default=.false.)

  ! Nonlocal module parameters

  call get_param(param_file, mdl, "CDRAG", cdrag, &

                 "CDRAG is the drag coefficient relating the magnitude of the velocity "//&

                 "field to the bottom stress.", units="nondim", default=0.003)

  call get_param(param_file, mdl, "MEKE_CDRAG", cs%cdrag, &

                 "Drag coefficient relating the magnitude of the velocity "//&

                 "field to the bottom stress in MEKE.", units="nondim", default=cdrag, scale=us%L_to_m*gv%m_to_H)

  call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)

  call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)

  if (cs%viscosity_coeff_Ku/=0. .and. .not. laplacian) call mom_error(fatal, &

                 "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.")

  if (cs%viscosity_coeff_Au/=0. .and. .not. biharmonic) call mom_error(fatal, &

                 "BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.")

  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)

  ! Identify if any lateral diffusive processes are active

  cs%kh_flux_enabled = .false.

  if ((cs%MEKE_KH >= 0.0)  .or. (cs%KhMEKE_FAC > 0.0) .or. (cs%MEKE_advection_factor > 0.0)) &

    cs%kh_flux_enabled = .true.

! Register fields for output from this module.

  cs%diag => diag

  cs%id_MEKE = register_diag_field('ocean_model', 'MEKE', diag%axesT1, time, &

     'Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=us%L_T_to_m_s**2)

  if (.not. allocated(meke%MEKE)) cs%id_MEKE = -1

  cs%id_Kh = register_diag_field('ocean_model', 'MEKE_KH', diag%axesT1, time, &

     'MEKE derived diffusivity', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

  if (.not. allocated(meke%Kh)) cs%id_Kh = -1

  cs%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, time, &

     'MEKE derived lateral viscosity', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

  if (.not. allocated(meke%Ku)) cs%id_Ku = -1

  cs%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, time, &

     'MEKE derived lateral biharmonic viscosity', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)

  if (.not. allocated(meke%Au)) cs%id_Au = -1

  cs%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, time, &

     'MEKE derived eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)

  if (.not. allocated(meke%MEKE)) cs%id_Ue = -1

  cs%id_Ub = register_diag_field('ocean_model', 'MEKE_Ub', diag%axesT1, time, &

     'MEKE derived bottom eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)

  if (.not. allocated(meke%MEKE)) cs%id_Ub = -1

  cs%id_Ut = register_diag_field('ocean_model', 'MEKE_Ut', diag%axesT1, time, &

     'MEKE derived barotropic eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)

  if (.not. allocated(meke%MEKE)) cs%id_Ut = -1

  cs%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, time, &

     'MEKE energy source', 'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  cs%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, time, &

     'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  cs%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, time, &

     '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)

  if (cs%MEKE_K4 >= 0.) &

    cs%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', &

        diag%axesT1, time, 'MEKE energy source from the biharmonic of MEKE', &

        'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  if (cs%MEKE_GMcoeff >= 0.) &

    cs%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', &

        diag%axesT1, time, 'MEKE energy source from the thickness mixing (GM scheme)', &

        'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  if (cs%MEKE_FrCoeff >= 0.) &

    cs%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', &

        diag%axesT1, time, 'MEKE energy source from the Laplacian of resolved flows', &

        'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  if (cs%MEKE_bhFrCoeff >= 0.) &

    cs%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', &

        diag%axesT1, time, 'MEKE energy source from the biharmonic of resolved flows', &

        'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)

  cs%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, time, &

     'MEKE decay rate', 's-1', conversion=us%s_to_T)

  cs%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, time, &

     'MEKE energy available from thickness mixing', &

     'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

  if (.not. allocated(meke%GM_src)) cs%id_GM_src = -1

  cs%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, time, &

     'MEKE energy available from momentum', &

     'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

  if (.not. allocated(meke%mom_src)) cs%id_mom_src = -1

  cs%id_mom_src_bh = register_diag_field('ocean_model', 'MEKE_mom_src_bh',diag%axesT1, time, &

     'MEKE energy available from the biharmonic dissipation of momentum', &

     'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

  if (.not. allocated(meke%mom_src_bh)) cs%id_mom_src_bh = -1

  cs%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, time, &

     'MEKE energy lost to GME backscatter', &

     'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

  if (.not. allocated(meke%GME_snk)) cs%id_GME_snk = -1

  cs%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, time, &

     'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)

  cs%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, time, &

     'Rhines length scale used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)

  cs%id_Leady = register_diag_field('ocean_model', 'MEKE_Leady', diag%axesT1, time, &

     'Eady length scale used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)

  cs%id_gamma_b = register_diag_field('ocean_model', 'MEKE_gamma_b', diag%axesT1, time, &

     'Ratio of bottom-projected eddy velocity to column-mean eddy velocity', 'nondim')

  cs%id_gamma_t = register_diag_field('ocean_model', 'MEKE_gamma_t', diag%axesT1, time, &

     'Ratio of barotropic eddy velocity to column-mean eddy velocity', 'nondim')

  if (cs%kh_flux_enabled) then

    cs%id_KhMEKE_u = register_diag_field('ocean_model', 'KHMEKE_u', diag%axesCu1, time, &

     'Zonal diffusivity of MEKE', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

    cs%id_KhMEKE_v = register_diag_field('ocean_model', 'KHMEKE_v', diag%axesCv1, time, &

     'Meridional diffusivity of MEKE', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

  endif

  if (cs%MEKE_equilibrium_restoring) then

    cs%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, time, &

     'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=us%L_T_to_m_s**2)

  endif

  cs%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=clock_routine)

  ! Detect whether this instance of MEKE_init() is at the beginning of a run

  ! or after a restart. If at the beginning, we will initialize MEKE to a local

  ! equilibrium.

  cs%initialize = .not.query_initialized(meke%MEKE, "MEKE", restart_cs)

  if (coldstart) cs%initialize = .false.

  if (cs%initialize) call mom_error(warning, &

                       "MEKE_init: Initializing MEKE with a local equilibrium balance.")

  if (allocated(meke%Le)) then

    if (.not.query_initialized(meke%Le, "MEKE_Le", restart_cs)) then

      !$OMP parallel do default(shared)

      do j=js,je ; do i=is,ie

        meke%Le(i,j) = sqrt(g%areaT(i,j))

      enddo ; enddo

    endif

  endif

  ! Set up group passes.  In the case of a restart, these fields need a halo update now.

  if (allocated(meke%MEKE)) then

    call create_group_pass(cs%pass_MEKE, meke%MEKE, g%Domain)

    if (allocated(meke%Kh_diff)) call create_group_pass(cs%pass_MEKE, meke%Kh_diff, g%Domain)

    if (.not.cs%initialize) call do_group_pass(cs%pass_MEKE, g%Domain)

  endif

  if (allocated(meke%Kh)) call create_group_pass(cs%pass_Kh, meke%Kh, g%Domain)

  if (allocated(meke%Ku)) call create_group_pass(cs%pass_Kh, meke%Ku, g%Domain)

  if (allocated(meke%Au)) call create_group_pass(cs%pass_Kh, meke%Au, g%Domain)

  if (allocated(meke%Le)) call create_group_pass(cs%pass_Kh, meke%Le, g%Domain)

  if (allocated(meke%Kh) .or. allocated(meke%Ku) .or. allocated(meke%Au) &

      .or. allocated(meke%Le)) &

    call do_group_pass(cs%pass_Kh, g%Domain)

end function meke_init

!> Initializer for the variant of MEKE that uses ML to predict eddy kinetic energy

subroutine ml_meke_init(diag, G, US, Time, param_file, dbcomms_CS, CS)

  type(diag_ctrl), target, intent(inout) :: diag       !< Diagnostics structure.

  type(ocean_grid_type),         intent(inout) :: G           !< The ocean's grid structure.

  type(unit_scale_type),         intent(in)    :: US          !< A dimensional unit scaling type

  type(time_type),               intent(in)    :: Time        !< The current model time.

  type(param_file_type),         intent(in)    :: param_file  !< Parameter file parser structure.

  type(dbcomms_cs_type),         intent(in)    :: dbcomms_CS  !< Control structure for database communication

  type(meke_cs),                 intent(inout) :: CS          !< Control structure for this module

  character(len=200)  :: inputdir, backend, model_filename

  integer :: db_return_code, batch_size

  character(len=40) :: mdl = "MOM_ML_MEKE"

  ! Store pointers in control structure

  write(cs%key_suffix, '(A,I6.6)') '_', pe_here()

  ! Put some basic information into the database

  db_return_code = 0

  db_return_code = cs%client%put_tensor("meta"//cs%key_suffix, &

    REAL([G%isd_global, G%idg_offset, G%jsd_global, G%jdg_offset]),[4]) + db_return_code

  db_return_code = cs%client%put_tensor("geolat"//cs%key_suffix, g%geoLatT, shape(g%geoLatT)) + db_return_code

  db_return_code = cs%client%put_tensor("geolon"//cs%key_suffix, g%geoLonT, shape(g%geoLonT)) + db_return_code

  db_return_code = cs%client%put_tensor("EKE_shape"//cs%key_suffix, shape(g%geolonT), [2]) + db_return_code

  if (cs%client%SR_error_parser(db_return_code)) call mom_error(fatal, "Putting metadata into the database failed")

  call read_param(param_file, "INPUTDIR", inputdir)

  inputdir = slasher(inputdir)

  call get_param(param_file, mdl, "BATCH_SIZE", batch_size, "Batch size to use for inference", default=1)

  call get_param(param_file, mdl, "EKE_BACKEND", backend, &

                 "The computational backend to use for EKE inference (CPU or GPU)", default="GPU")

  call get_param(param_file, mdl, "EKE_MODEL", model_filename, &

                 "Filename of the a saved pyTorch model to use", fail_if_missing = .true.)

  call get_param(param_file, mdl, "EKE_MAX", cs%eke_max, &

                 "Maximum value of EKE allowed when inferring EKE", &

                 units="m2 s-2", default=2., scale=us%m_s_to_L_T**2)

  ! Set the machine learning model

  if (dbcomms_cs%colocated) then

    if (modulo(pe_here(),dbcomms_cs%colocated_stride) == 0) then

      db_return_code = cs%client%set_model_from_file(cs%model_key, trim(inputdir)//trim(model_filename), &

                                                  "TORCH", backend, batch_size=batch_size)

    endif

  else

    if (is_root_pe()) then

      db_return_code = cs%client%set_model_from_file(cs%model_key, trim(inputdir)//trim(model_filename), &

                                                  "TORCH", backend, batch_size=batch_size)

    endif

  endif

  if (cs%client%SR_error_parser(db_return_code)) then

    call mom_error(fatal, "MEKE: set_model failed")

  endif

  call get_param(param_file, mdl, "ONLINE_ANALYSIS", cs%online_analysis, &

               "If true, post EKE used in MOM6 to the database for analysis", default=.true.)

  ! Set various clock ids

  cs%id_client_init   = cpu_clock_id('(ML_MEKE client init)', grain=clock_routine)

  cs%id_put_tensor    = cpu_clock_id('(ML_MEKE put tensor)', grain=clock_routine)

  cs%id_run_model     = cpu_clock_id('(ML_MEKE run model)', grain=clock_routine)

  cs%id_unpack_tensor = cpu_clock_id('(ML_MEKE unpack tensor )', grain=clock_routine)

  ! Diagnostics for ML_MEKE

  cs%id_mke = register_diag_field('ocean_model', 'MEKE_MKE', diag%axesT1, time, &

     'Surface mean (resolved) kinetic energy used in MEKE', 'm2 s-2', conversion=us%L_T_to_m_s**2)

  cs%id_slope_z= register_diag_field('ocean_model', 'MEKE_slope_z', diag%axesT1, time, &

     'Vertically averaged isopyncal slope magnitude used in MEKE', 'nondim', conversion=us%Z_to_L)

  cs%id_slope_x= register_diag_field('ocean_model', 'MEKE_slope_x', diag%axesCui, time, &

     'Isopycnal slope in the x-direction used in MEKE', 'nondim', conversion=us%Z_to_L)

  cs%id_slope_y= register_diag_field('ocean_model', 'MEKE_slope_y', diag%axesCvi, time, &

     'Isopycnal slope in the y-direction used in MEKE', 'nondim', conversion=us%Z_to_L)

  cs%id_rv = register_diag_field('ocean_model', 'MEKE_RV', diag%axesT1, time, &

     'Surface relative vorticity used in MEKE', 's-1', conversion=us%s_to_T)

end subroutine ml_meke_init

!> Calculate the various features used for the machine learning prediction

subroutine ml_meke_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, features_array)

  type(ocean_grid_type),                     intent(inout) :: G  !< Ocean grid

  type(verticalgrid_type),                   intent(in)    :: GV !< Ocean vertical grid structure

  type(unit_scale_type),                     intent(in)    :: US !< A dimensional unit scaling type

  type(meke_cs),                             intent(in)    :: CS !< Control structure for MEKE

  real, dimension(SZI_(G),SZJ_(G)),          intent(in   ) :: Rd_dx_h !< Rossby radius of deformation over

                                                                 !! the grid length scale [nondim]

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u  !< Zonal velocity [L T-1 ~> m s-1]

  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v  !< Meridional velocity [L T-1 ~> m s-1]

  type(thermo_var_ptrs),                     intent(in)    :: tv !< Type containing thermodynamic variables

  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h  !< Layer thickness [H ~> m or kg m-2].

  real,                                      intent(in)    :: dt !< Model(baroclinic) time-step [T ~> s].

  real(kind=real32), dimension(SIZE(h),num_features), intent(  out) :: features_array

                                                                 !< The array of features needed for machine

                                                                 !! learning inference, with different units

                                                                 !! for the various subarrays [various]

  real, dimension(SZI_(G),SZJ_(G)) :: mke      ! Surface kinetic energy per unit mass [L2 T-2 ~> m2 s-2]

  real, dimension(SZI_(G),SZJ_(G)) :: slope_z  ! Vertically averaged isoneutral slopes [Z L-1 ~> nondim]

  real, dimension(SZIB_(G),SZJB_(G)) :: rv_z   ! Surface relative vorticity [T-1 ~> s-1]

  real, dimension(SZIB_(G),SZJB_(G)) :: rv_z_t ! Surface relative vorticity interpolated to tracer points [T-1 ~> s-1]

  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)) :: h_u ! Thickness at u point [H ~> m or kg m-2]

  real, dimension(SZI_(G),SZJB_(G), SZK_(G)) :: h_v ! Thickness at v point [H ~> m or kg m-2]

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: slope_x ! Isoneutral slope at U point [Z L-1 ~> nondim]

  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: slope_y ! Isoneutral slope at V point [Z L-1 ~> nondim]

  real, dimension(SZIB_(G),SZJ_(G)) :: slope_x_vert_avg ! Isoneutral slope at U point [Z L-1 ~> nondim]

  real, dimension(SZI_(G),SZJB_(G)) :: slope_y_vert_avg ! Isoneutral slope at V point [Z L-1 ~> nondim]

  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) ::  e ! The interface heights relative to mean sea level [Z ~> m].

  real :: slope_t  ! Slope interpolated to thickness points [Z L-1 ~> nondim]

  real :: u_t, v_t ! u and v interpolated to thickness points [L T-1 ~> m s-1]

  real :: dvdx, dudy ! Components of relative vorticity [T-1 ~> s-1]

  real :: a_e, a_w, a_n, a_s ! Fractional areas of neighboring cells for interpolating velocities [nondim]

  real :: Idenom    ! A normalizing factor in calculating weighted averages of areas [L-2 ~> m-2]

  real :: sum_area  ! A sum of adjacent cell areas [L2 ~> m2]

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  ! Calculate various features for used to infer eddy kinetic energy

  ! Linear interpolation to estimate thickness at a velocity points

  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1

    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

    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

  enddo ; enddo ; enddo

  call find_eta(h, tv, g, gv, us, e, halo_size=2)

  ! Note the hard-coded dimenisional constant in the following line.

  call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*1.e-7*gv%m2_s_to_HZ_T, .false., slope_x, slope_y)

  call pass_vector(slope_x, slope_y, g%Domain)

  do j=js-1,je+1 ; do i=is-1,ie+1

    slope_x_vert_avg(i,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), gv%H_subroundoff)

    slope_y_vert_avg(i,j) = vertical_average_interface(slope_y(i,j,:), h_v(i,j,:), gv%H_subroundoff)

  enddo ; enddo

  slope_z(:,:) = 0.

  call pass_vector(slope_x_vert_avg, slope_y_vert_avg, g%Domain)

  do j=js,je ; do i=is,ie

    ! Calculate weights for interpolation from velocity points to h points

    sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)

    if (sum_area>0.0) then

      idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)

      a_w = g%areaCu(i-1,j) * idenom

      a_e = g%areaCu(i,j) * idenom

    else

      a_w = 0.0 ; a_e = 0.0

    endif

    sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)

    if (sum_area>0.0) then

      idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)

      a_s = g%areaCv(i,j-1) * idenom

      a_n = g%areaCv(i,j) * idenom

    else

      a_s = 0.0 ; a_n = 0.0

    endif

    ! Calculate mean kinetic energy

    u_t = (a_e*u(i,j,1)) + (a_w*u(i-1,j,1))

    v_t = (a_n*v(i,j,1)) + (a_s*v(i,j-1,1))

    mke(i,j) = 0.5*( (u_t*u_t) + (v_t*v_t) )

    ! Calculate the magnitude of the slope

    slope_t = slope_x_vert_avg(i,j)*a_e+slope_x_vert_avg(i-1,j)*a_w

    slope_z(i,j) = sqrt(slope_t*slope_t)

    slope_t = slope_y_vert_avg(i,j)*a_n+slope_y_vert_avg(i,j-1)*a_s

    slope_z(i,j) = 0.5*(slope_z(i,j) + sqrt(slope_t*slope_t))*g%mask2dT(i,j)

  enddo ; enddo

  call pass_var(slope_z, g%Domain)

  ! Calculate relative vorticity

  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1

    dvdx = ((v(i+1,j,1)*g%dyCv(i+1,j)) - (v(i,j,1)*g%dyCv(i,j)))

    dudy = ((u(i,j+1,1)*g%dxCu(i,j+1)) - (u(i,j,1)*g%dxCu(i,j)))

    ! Assumed no slip

    rv_z(i,j) = (2.0-g%mask2dBu(i,j)) * (dvdx - dudy) * g%IareaBu(i,j)

  enddo ; enddo

  ! Interpolate RV to t-point, revisit this calculation to include metrics

  do j=js,je ; do i=is,ie

    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))

  enddo ; enddo

  ! Construct the feature array

  features_array(:,mke_idx) = pack(mke,.true.)

  features_array(:,slope_z_idx) = pack(slope_z,.true.)

  features_array(:,rd_dx_z_idx) = pack(rd_dx_h,.true.)

  features_array(:,rv_idx) = pack(rv_z_t,.true.)

  if (cs%id_rv>0) call post_data(cs%id_rv, rv_z, cs%diag)

  if (cs%id_mke>0) call post_data(cs%id_mke, mke, cs%diag)

  if (cs%id_slope_z>0) call post_data(cs%id_slope_z, slope_z, cs%diag)

  if (cs%id_slope_x>0) call post_data(cs%id_slope_x, slope_x, cs%diag)

  if (cs%id_slope_y>0) call post_data(cs%id_slope_y, slope_y, cs%diag)

end subroutine ml_meke_calculate_features

!> Use the machine learning interface to predict EKE

subroutine predict_meke(G, US, CS, npts, Time, features_array, MEKE)

  type(ocean_grid_type),                                 intent(inout) :: G  !< Ocean grid

  type(unit_scale_type),                                 intent(in)    :: US   !< A dimensional unit scaling type

  type(meke_cs),                                         intent(in   ) :: CS !< Control structure for MEKE

  integer,                                               intent(in   ) :: npts !< Number of T-grid cells on the local

                                                                               !! domain

  type(time_type),                                       intent(in   ) :: Time !< The current model time

  real(kind=real32), dimension(npts,num_features),       intent(in   ) :: features_array

                                                                          !< The array of features needed for machine

                                                                          !! learning inference, with different units

                                                                          !! for the various subarrays [various]

  real, dimension(SZI_(G),SZJ_(G)),                      intent(  out) :: MEKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2]

  ! Local variables

  integer :: db_return_code

  character(len=255), dimension(1) :: model_out, model_in

  character(len=255) :: time_suffix

  real(kind=real32), dimension(SIZE(MEKE)) :: meke_vec ! A one-dimensional array of the natural log of eddy kinetic

                                                       ! energy in mks units [m2 s-2]

  real, dimension(size(MEKE,1),size(MEKE,2)) :: ln_MEKE  ! the natural log of eddy kinetic energy

                                                       ! in mks units [m2 s-2]

  real, dimension(size(MEKE,1),size(MEKE,2)) :: MEKE_mks  ! The eddy kinetic energy in mks units [m2 s-2]

  integer :: i, j, is, ie, js, je

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

!> Use the database client to call a machine learning model to predict eddy kinetic energy

  call cpu_clock_begin(cs%id_put_tensor)

  db_return_code = cs%client%put_tensor("features"//cs%key_suffix, features_array, shape(features_array))

  call cpu_clock_end(cs%id_put_tensor)

  ! Run the ML model to predict EKE and return the result

  model_out(1) = "EKE"//cs%key_suffix

  model_in(1) = "features"//cs%key_suffix

  call cpu_clock_begin(cs%id_run_model)

  db_return_code = cs%client%run_model(cs%model_key, model_in, model_out)

  call cpu_clock_end(cs%id_run_model)

  if (cs%client%SR_error_parser(db_return_code)) then

    call mom_error(fatal, "MEKE: run_model failed")

  endif

  call cpu_clock_begin(cs%id_unpack_tensor)

  db_return_code = cs%client%unpack_tensor( model_out(1), meke_vec, shape(meke_vec) )

  call cpu_clock_end(cs%id_unpack_tensor)

  ln_meke = reshape(meke_vec, shape(meke))

  ! Zero out the halos.  These will usually be reset by the pass_var in a few lines.

  meke_mks(:,:) = 0.0

  do j=js,je ; do i=is,ie

    meke_mks(i,j) = min(exp(ln_meke(i,j)), us%L_T_to_m_s**2*cs%eke_max)

  enddo ; enddo

  call pass_var(meke_mks, g%Domain, halo=1)

  if (cs%online_analysis) then

    write(time_suffix,"(F16.0)") time_type_to_real(time)

    db_return_code = cs%client%put_tensor(trim("EKE_")//trim(adjustl(time_suffix))//cs%key_suffix, &

                                          meke_mks, shape(meke))

  endif

  ! Copy MEKE_mks into the argument in rescaled units.

  ! MEKE(:,:) = 0.0  ! This would fill in the wider halos of this intent(out) array.

  do j=js-1,je+1 ; do i=is-1,ie+1

    meke(i,j) = us%m_s_to_L_T**2 * meke_mks(i,j)

  enddo ; enddo

end subroutine predict_meke

!> Compute average of interface quantities weighted by the thickness of the surrounding

!! layers [arbitrary]

real function vertical_average_interface(h, w, h_min)

  real, dimension(:), intent(in) :: h  !< Layer Thicknesses [H ~> m or kg m-2]

  real, dimension(:), intent(in) :: w  !< Quantity to average [arbitrary]

  real, intent(in) :: h_min !< The vanishingly small layer thickness [H ~> m or kg m-2]

  real :: htot  ! Twice the sum of the layer thicknesses interpolated to interior interfaces [H ~> m or kg m-2]

  real :: inv_htot ! The inverse of htot  [H-1 ~> m-1 or m2 kg-1]

  integer :: k, nk

  nk = size(h)

  htot = h_min

  do k=2,nk

    htot = htot + (h(k-1)+h(k))

  enddo

  inv_htot = 1./htot

  vertical_average_interface = 0.

  do k=2,nk

    vertical_average_interface = vertical_average_interface + (w(k)*(h(k-1)+h(k)))*inv_htot

  enddo

end function vertical_average_interface

!> Allocates memory and register restart fields for the MOM_MEKE module.

subroutine meke_alloc_register_restart(HI, US, param_file, MEKE, restart_CS)

! Arguments

  type(hor_index_type),  intent(in)    :: hi         !< Horizontal index structure

  type(unit_scale_type), intent(in)    :: us         !< A dimensional unit scaling type

  type(param_file_type), intent(in)    :: param_file !< Parameter file parser structure.

  type(meke_type),       intent(inout) :: meke       !< MEKE fields

  type(mom_restart_cs),  intent(inout) :: restart_cs !< MOM restart control struct

  ! Local variables

  real :: meke_gmcoeff, meke_frcoeff, meke_bhfrcoeff, meke_gmecoeff  ! Coefficients for various terms [nondim]

  real :: meke_khcoeff, meke_visccoeff_ku, meke_visccoeff_au  ! Coefficients for various terms [nondim]

  logical :: use_kh_in_meke

  logical :: usemeke

  logical :: sqg_use_meke

  integer :: isd, ied, jsd, jed

! Determine whether this module will be used

  usemeke = .false. ; call read_param(param_file,"USE_MEKE",usemeke)

! Read these parameters to determine what should be in the restarts

  meke_gmcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",meke_gmcoeff)

  meke_frcoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",meke_frcoeff)

  meke_bhfrcoeff = -1. ; call read_param(param_file,"MEKE_bhFRCOEFF",meke_bhfrcoeff)

  meke_gmecoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",meke_gmecoeff)

  meke_khcoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",meke_khcoeff)

  meke_visccoeff_ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",meke_visccoeff_ku)

  meke_visccoeff_au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",meke_visccoeff_au)

  use_kh_in_meke = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", use_kh_in_meke)

  sqg_use_meke = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_meke)

  if (.not. usemeke) return

! Allocate memory

  call mom_mesg("MEKE_alloc_register_restart: allocating and registering", 5)

  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed

  allocate(meke%MEKE(isd:ied,jsd:jed), source=0.0)

  call register_restart_field(meke%MEKE, "MEKE", .false., restart_cs, &

           longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=us%L_T_to_m_s**2)

  if (meke_gmcoeff>=0.) allocate(meke%GM_src(isd:ied,jsd:jed), source=0.0)

  if (meke_frcoeff>=0. .or. meke_bhfrcoeff>=0. .or. meke_gmecoeff>=0.) &

    allocate(meke%mom_src(isd:ied,jsd:jed), source=0.0)

  if (meke_bhfrcoeff >= 0.) &

    allocate(meke%mom_src_bh(isd:ied,jsd:jed), source=0.0)

  if (meke_frcoeff<0.) meke_frcoeff = 0.

  if (meke_bhfrcoeff<0.) meke_bhfrcoeff = 0.

  if (meke_gmecoeff>=0.) allocate(meke%GME_snk(isd:ied,jsd:jed), source=0.0)

  if (meke_khcoeff>=0.) then

    allocate(meke%Kh(isd:ied,jsd:jed), source=0.0)

    call register_restart_field(meke%Kh, "MEKE_Kh", .false., restart_cs, &

             longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy", &

             units="m2 s-1", conversion=us%L_to_m**2*us%s_to_T)

  endif

  allocate(meke%Rd_dx_h(isd:ied,jsd:jed), source=0.0)

  if (meke_visccoeff_ku/=0.) then

    allocate(meke%Ku(isd:ied,jsd:jed), source=0.0)

    call register_restart_field(meke%Ku, "MEKE_Ku", .false., restart_cs, &

             longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", &

             units="m2 s-1", conversion=us%L_to_m**2*us%s_to_T)

  endif

  if (sqg_use_meke) then

    allocate(meke%Le(isd:ied,jsd:jed), source=0.0)

    call register_restart_field(meke%Le, "MEKE_Le", .false., restart_cs, &

             longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", &

             units="m", conversion=us%L_to_m)

  endif

  if (use_kh_in_meke) then

    allocate(meke%Kh_diff(isd:ied,jsd:jed), source=0.0)

    call register_restart_field(meke%Kh_diff, "MEKE_Kh_diff", .false., restart_cs, &

             longname="Copy of thickness diffusivity for diffusing MEKE", &

             units="m2 s-1", conversion=us%L_to_m**2*us%s_to_T)

  endif

  if (meke_visccoeff_au/=0.) then

    allocate(meke%Au(isd:ied,jsd:jed), source=0.0)

    call register_restart_field(meke%Au, "MEKE_Au", .false., restart_cs, &

             longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy", &

             units="m4 s-1", conversion=us%L_to_m**4*us%s_to_T)

  endif

end subroutine meke_alloc_register_restart

!> Deallocates any variables allocated in MEKE_alloc_register_restart.

subroutine meke_end(MEKE)

  type(meke_type), intent(inout) :: meke !< A structure with MEKE-related fields.

  ! NOTE: MEKE will always be allocated by MEKE_init, even if MEKE is disabled.

  !  So these must all be conditional, even though MEKE%MEKE and MEKE%Rd_dx_h

  !  are always allocated (when MEKE is enabled)

  if (allocated(meke%Au)) deallocate(meke%Au)

  if (allocated(meke%Kh_diff)) deallocate(meke%Kh_diff)

  if (allocated(meke%Ku)) deallocate(meke%Ku)

  if (allocated(meke%Rd_dx_h)) deallocate(meke%Rd_dx_h)

  if (allocated(meke%Kh)) deallocate(meke%Kh)

  if (allocated(meke%GME_snk)) deallocate(meke%GME_snk)

  if (allocated(meke%mom_src)) deallocate(meke%mom_src)

  if (allocated(meke%mom_src_bh)) deallocate(meke%mom_src_bh)

  if (allocated(meke%GM_src)) deallocate(meke%GM_src)

  if (allocated(meke%MEKE)) deallocate(meke%MEKE)

  if (allocated(meke%Le)) deallocate(meke%Le)

end subroutine meke_end

!> \namespace mom_meke

!!

!! \section section_MEKE The Mesoscale Eddy Kinetic Energy (MEKE) framework

!!

!! The MEKE framework accounts for the mean potential energy removed by

!! the first order closures used to parameterize mesoscale eddies.

!! It requires closure at the second order, namely dissipation and transport

!! of eddy energy.

!!

!! Monitoring the sub-grid scale eddy energy budget provides a means to predict

!! a sub-grid eddy-velocity scale which can be used in the lower order closures.

!!

!! \subsection section_MEKE_equations MEKE equations

!!

!! The eddy kinetic energy equation is:

!! \f[ \partial_{\tilde{t}} E =

!!   \overbrace{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v

!!             }^\text{sources}

!! - \overbrace{ ( \lambda + C_d | U_d | \gamma_b^2 ) E

!!             }^\text{local dissipation}

!! + \overbrace{ \nabla \cdot ( ( \kappa_E + \gamma_M \kappa_M ) \nabla E

!!                              - \kappa_4 \nabla^3 E )

!!             }^\text{smoothing}

!! \f]

!! where \f$ E \f$ is the eddy kinetic energy (variable <code>MEKE</code>) with units of

!! m<sup>2</sup>s<sup>-2</sup>,

!! and \f$\tilde{t} = a t\f$ is a scaled time. The non-dimensional factor

!! \f$ a\geq 1 \f$ is used to accelerate towards equilibrium.

!!

!! The MEKE equation is two-dimensional and obtained by depth averaging the

!! the three-dimensional eddy energy equation. In the following expressions

!! \f$ \left< \phi \right> = \frac{1}{H} \int^\eta_{-D} \phi \, dz \f$ maps

!! three dimensional terms into the two-dimensional quantities needed.

!!

!! \subsubsection section_MEKE_source_terms MEKE source terms

!!

!! The source term \f$ \dot{E}_b \f$ is a constant background source

!! of energy intended to avoid the limit \f$E\rightarrow 0\f$.

!!

!! The "GM" source term

!! \f[ \dot{E}_\eta = - \left< \overline{w^\prime b^\prime} \right>

!! = \left< \kappa_h N^2S^2 \right>

!! \approx \left< \kappa_h g\prime |\nabla_\sigma \eta|^2 \right>\f]

!! equals the mean potential energy removed by the Gent-McWilliams closure,

!! and is excluded/included in the MEKE budget by the efficiency parameter

!! \f$ \gamma_\eta \in [0,1] \f$.

!!

!! The "frictional" source term

!! \f[ \dot{E}_{v} = \left<  \partial_i u_j \tau_{ij} \right> \f]

!! equals the mean kinetic energy removed by lateral viscous fluxes, and

!! is excluded/included in the MEKE budget by the efficiency parameter

!! \f$ \gamma_v \in [0,1] \f$.

!!

!! \subsubsection section_MEKE_dissipation_terms MEKE dissipation terms

!!

!! The local dissipation of \f$ E \f$ is parameterized through a linear

!! damping, \f$\lambda\f$, and bottom drag, \f$ C_d | U_d | \gamma_b^2 \f$.

!! The \f$ \gamma_b \f$ accounts for the weak projection of the column-mean

!! eddy velocity to the bottom. In other words, the bottom velocity is

!! estimated as \f$ \gamma_b U_e \f$.

!! The bottom drag coefficient, \f$ C_d \f$ is the same as that used in the bottom

!! friction in the mean model equations.

!!

!! The bottom drag velocity scale, \f$ U_d \f$, has contributions from the

!! resolved state and \f$ E \f$:

!! \f[ U_d = \sqrt{ U_b^2 + |u|^2_{z=-D} + |\gamma_b U_e|^2 } .\f]

!! where the eddy velocity scale, \f$ U_e \f$, is given by:

!! \f[ U_e = \sqrt{ 2 E } .\f]

!! \f$ U_b \f$ is a constant background bottom velocity scale and is

!! typically not used (i.e. set to zero).

!!

!! Following \cite jansen2015, the projection of eddy energy on to the bottom

!! is given by the ratio of bottom energy to column mean energy:

!! \f[

!! \gamma_b^2  = \frac{E_b}{E} = \gamma_{d0}

!!    + \left( 1 + c_{b} \frac{L_d}{L_f} \right)^{-\frac{4}{5}}

!! ,

!! \f]

!! \f[

!! \gamma_b^2  \leftarrow  \max{\left( \gamma_b^2, \gamma_{min}^2 \right)}

!! .

!! \f]

!!

!! \subsection section_MEKE_smoothing MEKE smoothing terms

!!

!! \f$ E \f$ is laterally diffused by a diffusivity \f$ \kappa_E + \gamma_M

!! \kappa_M \f$ where \f$ \kappa_E \f$ is a constant diffusivity and the term

!! \f$ \gamma_M \kappa_M \f$ is a "self diffusion" using the diffusivity

!! calculated in the section \ref section_MEKE_diffusivity.

!! \f$ \kappa_4 \f$ is a constant bi-harmonic diffusivity.

!!

!! \subsection section_MEKE_diffusivity Diffusivity derived from MEKE

!!

!! The predicted eddy velocity scale, \f$ U_e \f$, can be combined with a

!! mixing length scale to form a diffusivity.

!! The primary use of a MEKE derived diffusivity is for use in thickness

!! diffusion (module mom_thickness_diffuse) and optionally in along

!! isopycnal mixing of tracers (module mom_tracer_hor_diff).

!! The original form used (enabled with MEKE_OLD_LSCALE=True):

!!

!! \f[  \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \f]

!!

!! where \f$ A_\Delta \f$ is the area of the grid cell.

!! Following \cite jansen2015, we now use

!!

!! \f[  \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \f]

!!

!! where \f$ \gamma_\kappa \in [0,1] \f$ is a non-dimensional factor and,

!! following \cite jansen2015, \f$\gamma_t^2\f$ is the ratio of barotropic

!! eddy energy to column mean eddy energy given by

!! \f[

!! \gamma_t^2  = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}}

!! ,

!! \f]

!! \f[

!! \gamma_t^2  \leftarrow  \max{\left( \gamma_t^2, \gamma_{min}^2 \right)}

!! .

!! \f]

!!

!! The length-scale is a configurable combination of multiple length scales:

!!

!! \f[

!! l_M = \left(

!!       \frac{\alpha_d}{L_d}

!!     + \frac{\alpha_f}{L_f}

!!     + \frac{\alpha_R}{L_R}

!!     + \frac{\alpha_e}{L_e}

!!     + \frac{\alpha_\Delta}{L_\Delta}

!!     + \frac{\delta[L_c]}{L_c}

!!       \right)^{-1}

!! \f]

!!

!! where

!!

!! \f{eqnarray*}{

!! L_d & = & \sqrt{\frac{c_g^2}{f^2+2\beta c_g}} \sim \frac{ c_g }{f} \\\\

!! L_R & = & \sqrt{\frac{U_e}{\beta^*}} \\\\

!! L_e & = & \frac{U_e}{|S| N} \\\\

!! L_f & = & \frac{H}{c_d} \\\\

!! L_\Delta & = & \sqrt{A_\Delta} .

!! \f}

!!

!! \f$L_c\f$ is a constant and \f$\delta[L_c]\f$ is the impulse function so that the term

!! \f$\frac{\delta[L_c]}{L_c}\f$ evaluates to \f$\frac{1}{L_c}\f$ when \f$L_c\f$ is non-zero

!! but is dropped if \f$L_c=0\f$.

!!

!! \f$\beta^*\f$ is the effective \f$\beta\f$ that combines both the planetary vorticity

!! gradient (i.e. \f$\beta=\nabla f\f$) and the topographic \f$\beta\f$ effect,

!! with the latter weighed by a weighting constant, \f$c_\beta\f$, that varies

!! from 0 to 1, so that \f$c_\beta=0\f$ means the topographic \f$\beta\f$ effect is ignored,

!! while \f$c_\beta=1\f$ means it is fully considered. The new \f$\beta^*\f$ therefore

!! takes the form of

!!

!! \f[

!! \beta^* = \sqrt{( \partial_xf - c_\beta\frac{f}{D}\partial_xD )^2 +

!!           ( \partial_yf - c_\beta\frac{f}{D}\partial_yD )^2}

!! \f]

!! where \f$D\f$ is water column depth at T points.

!!

!! \subsection section_MEKE_viscosity Viscosity derived from MEKE

!!

!! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be

!! used to form a harmonic eddy viscosity,

!!

!! \f[  \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta }  \f]

!!

!! as well as a biharmonic eddy viscosity,

!!

!! \f[  \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 }  \f]

!!

!! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance

!!

!! Note that in steady-state (or when \f$ a>>1 \f$) and there is no

!! diffusion of \f$ E \f$ then

!! \f[ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +

!!               \gamma_v \dot{E}_v }{ \lambda + C_d|U_d|\gamma_b^2 } . \f]

!!

!! In the linear drag limit, where

!! \f$ U_e << \min(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$, the equilibrium becomes

!! \f$ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +

!!               \gamma_v \dot{E}_v }{ \lambda + C_d \sqrt{ U_b^2 + |u|^2_{z=-D} } } \f$.

!!

!! In the nonlinear drag limit, where \f$ U_e >> \max(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$,

!! the equilibrium becomes

!! \f$ \overline{E} \approx \left( \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +

!!               \gamma_v \dot{E}_v }{ \sqrt{2} C_d \gamma_b^3 } \right)^\frac{2}{3} \f$.

!!

!! \subsubsection section_MEKE_module_parameters MEKE module parameters

!!

!! | Symbol                | Module parameter |

!! | ------                | --------------- |

!! | -                     | <code>USE_MEKE</code> |

!! | \f$ a \f$             | <code>MEKE_DTSCALE</code> |

!! | \f$ \dot{E}_b \f$     | <code>MEKE_BGSRC</code> |

!! | \f$ \gamma_\eta \f$   | <code>MEKE_GMCOEFF</code> |

!! | \f$ \gamma_v \f$      | <code>MEKE_FrCOEFF</code> |

!! | \f$ \lambda \f$       | <code>MEKE_DAMPING</code> |

!! | \f$ U_b \f$           | <code>MEKE_USCALE</code> |

!! | \f$ \gamma_{d0} \f$   | <code>MEKE_CD_SCALE</code> |

!! | \f$ c_{b} \f$         | <code>MEKE_CB</code> |

!! | \f$ c_{t} \f$         | <code>MEKE_CT</code> |

!! | \f$ \kappa_E \f$      | <code>MEKE_KH</code> |

!! | \f$ \kappa_4 \f$      | <code>MEKE_K4</code> |

!! | \f$ \gamma_\kappa \f$ | <code>MEKE_KHCOEFF</code> |

!! | \f$ \gamma_M \f$      | <code>MEKE_KHMEKE_FAC</code> |

!! | \f$ \gamma_u \f$      | <code>MEKE_VISCOSITY_COEFF_KU</code> |

!! | \f$ \gamma_4 \f$      | <code>MEKE_VISCOSITY_COEFF_AU</code> |

!! | \f$ \gamma_{min}^2 \f$| <code>MEKE_MIN_GAMMA2</code> |

!! | \f$ \alpha_d \f$      | <code>MEKE_ALPHA_DEFORM</code> |

!! | \f$ \alpha_f \f$      | <code>MEKE_ALPHA_FRICT</code> |

!! | \f$ \alpha_R \f$      | <code>MEKE_ALPHA_RHINES</code> |

!! | \f$ \alpha_e \f$      | <code>MEKE_ALPHA_EADY</code> |

!! | \f$ \alpha_\Delta \f$ | <code>MEKE_ALPHA_GRID</code> |

!! | \f$ L_c \f$           | <code>MEKE_FIXED_MIXING_LENGTH</code> |

!! | \f$ c_\beta \f$       | <code>MEKE_TOPOGRAPHIC_BETA</code> |

!! | -                     | <code>MEKE_KHTH_FAC</code> |

!! | -                     | <code>MEKE_KHTR_FAC</code> |

!!

!! | Symbol                | Model parameter |

!! | ------                | --------------- |

!! | \f$ C_d \f$           | <code>CDRAG</code> |

!!

!! \subsection section_MEKE_references References

!!

!! Jansen, M. F., A. J. Adcroft, R. Hallberg, and I. M. Held, 2015: Parameterization of eddy fluxes based on a

!! mesoscale energy budget. Ocean Modelling, 92, 28--41, http://doi.org/10.1016/j.ocemod.2015.05.007 .

!!

!! Marshall, D. P., and A. J. Adcroft, 2010: Parameterization of ocean eddies: Potential vorticity mixing, energetics

!! and Arnold first stability theorem. Ocean Modelling, 32, 188--204, http://doi.org/10.1016/j.ocemod.2010.02.001 .

end module mom_meke