MOM6
MOM_MEKE.F90
1 !> Implements the Mesoscale Eddy Kinetic Energy framework
2 !! with topographic beta effect included in computing beta in Rhines scale
3 
4 module mom_meke
5 
6 ! This file is part of MOM6. See LICENSE.md for the license.
7 
8 use mom_debugging, only : hchksum, uvchksum
9 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
10 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
11 use mom_diag_mediator, only : diag_ctrl, time_type
12 use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
13 use mom_error_handler, only : mom_error, fatal, warning, note, mom_mesg
14 use mom_file_parser, only : read_param, get_param, log_version, param_file_type
15 use mom_grid, only : ocean_grid_type
16 use mom_hor_index, only : hor_index_type
17 use mom_io, only : vardesc, var_desc
18 use mom_restart, only : mom_restart_cs, register_restart_field, query_initialized
19 use mom_unit_scaling, only : unit_scale_type
20 use mom_variables, only : vertvisc_type
21 use mom_verticalgrid, only : verticalgrid_type
22 use mom_meke_types, only : meke_type
23 
24 implicit none ; private
25 
26 #include <MOM_memory.h>
27 
28 public step_forward_meke, meke_init, meke_alloc_register_restart, meke_end
29 
30 !> Control structure that contains MEKE parameters and diagnostics handles
31 type, public :: meke_cs ; private
32  ! Parameters
33  real, dimension(:,:), pointer :: equilibrium_value => null() !< The equilbrium value
34  !! of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2]
35  real :: meke_frcoeff !< Efficiency of conversion of ME into MEKE [nondim]
36  real :: meke_gmcoeff !< Efficiency of conversion of PE into MEKE [nondim]
37  real :: meke_gmecoeff !< Efficiency of conversion of MEKE into ME by GME [nondim]
38  real :: meke_damping !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1].
39  real :: meke_cd_scale !< The ratio of the bottom eddy velocity to the column mean
40  !! eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1
41  !! to account for the surface intensification of MEKE.
42  real :: meke_cb !< Coefficient in the \f$\gamma_{bot}\f$ expression [nondim]
43  real :: meke_min_gamma!< Minimum value of gamma_b^2 allowed [nondim]
44  real :: meke_ct !< Coefficient in the \f$\gamma_{bt}\f$ expression [nondim]
45  logical :: visc_drag !< If true use the vertvisc_type to calculate bottom drag.
46  logical :: meke_geometric !< If true, uses the GM coefficient formulation from the GEOMETRIC
47  !! framework (Marshall et al., 2012)
48  real :: meke_geometric_alpha !< The nondimensional coefficient governing the efficiency of the
49  !! GEOMETRIC thickness diffusion.
50  logical :: meke_equilibrium_alt !< If true, use an alternative calculation for the
51  !! equilibrium value of MEKE.
52  logical :: meke_equilibrium_restoring !< If true, restore MEKE back to its equilibrium value,
53  !! which is calculated at each time step.
54  logical :: gm_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
55  !! than the streamfunction for the MEKE GM source term.
56  logical :: rd_as_max_scale !< If true the length scale can not exceed the
57  !! first baroclinic deformation radius.
58  logical :: use_old_lscale !< Use the old formula for mixing length scale.
59  logical :: use_min_lscale !< Use simple minimum for mixing length scale.
60  real :: cdrag !< The bottom drag coefficient for MEKE [nondim].
61  real :: meke_bgsrc !< Background energy source for MEKE [L2 T-3 ~> W kg-1] (= m2 s-3).
62  real :: meke_dtscale !< Scale factor to accelerate time-stepping [nondim]
63  real :: meke_khcoeff !< Scaling factor to convert MEKE into Kh [nondim]
64  real :: meke_uscale !< MEKE velocity scale for bottom drag [L T-1 ~> m s-1]
65  real :: meke_kh !< Background lateral diffusion of MEKE [L2 T-1 ~> m2 s-1]
66  real :: meke_k4 !< Background bi-harmonic diffusivity (of MEKE) [L4 T-1 ~> m4 s-1]
67  real :: khmeke_fac !< A factor relating MEKE%Kh to the diffusivity used for
68  !! MEKE itself [nondim].
69  real :: viscosity_coeff_ku !< The scaling coefficient in the expression for
70  !! viscosity used to parameterize lateral harmonic momentum mixing
71  !! by unresolved eddies represented by MEKE.
72  real :: viscosity_coeff_au !< The scaling coefficient in the expression for
73  !! viscosity used to parameterize lateral biharmonic momentum mixing
74  !! by unresolved eddies represented by MEKE.
75  real :: lfixed !< Fixed mixing length scale [L ~> m].
76  real :: adeform !< Weighting towards deformation scale of mixing length [nondim]
77  real :: arhines !< Weighting towards Rhines scale of mixing length [nondim]
78  real :: africt !< Weighting towards frictional arrest scale of mixing length [nondim]
79  real :: aeady !< Weighting towards Eady scale of mixing length [nondim]
80  real :: agrid !< Weighting towards grid scale of mixing length [nondim]
81  real :: meke_advection_factor !< A scaling in front of the advection of MEKE [nondim]
82  real :: meke_topographic_beta !< Weight for how much topographic beta is considered
83  !! when computing beta in Rhines scale [nondim]
84  real :: meke_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1].
85 
86  logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled.
87  logical :: initialize !< If True, invokes a steady state solver to calculate MEKE.
88  logical :: debug !< If true, write out checksums of data for debugging
89 
90  type(diag_ctrl), pointer :: diag => null() !< A type that regulates diagnostics output
91  !>@{ Diagnostic handles
92  integer :: id_meke = -1, id_ue = -1, id_kh = -1, id_src = -1
93  integer :: id_ub = -1, id_ut = -1
94  integer :: id_gm_src = -1, id_mom_src = -1, id_gme_snk = -1, id_decay = -1
95  integer :: id_khmeke_u = -1, id_khmeke_v = -1, id_ku = -1, id_au = -1
96  integer :: id_le = -1, id_gamma_b = -1, id_gamma_t = -1
97  integer :: id_lrhines = -1, id_leady = -1
98  integer :: id_meke_equilibrium = -1
99  !>@}
100 
101  ! Infrastructure
102  integer :: id_clock_pass !< Clock for group pass calls
103  type(group_pass_type) :: pass_meke !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff
104  type(group_pass_type) :: pass_kh !< Group halo pass handle for MEKE%Kh, MEKE%Ku, and/or MEKE%Au
105 end type meke_cs
106 
107 contains
108 
109 !> Integrates forward-in-time the MEKE eddy energy equation.
110 !! See \ref section_MEKE_equations.
111 subroutine step_forward_meke(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, hv)
112  type(meke_type), pointer :: meke !< MEKE data.
113  type(ocean_grid_type), intent(inout) :: g !< Ocean grid.
114  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure.
115  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
116  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
117  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: sn_u !< Eady growth rate at u-points [T-1 ~> s-1].
118  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: sn_v !< Eady growth rate at v-points [T-1 ~> s-1].
119  type(vertvisc_type), intent(in) :: visc !< The vertical viscosity type.
120  real, intent(in) :: dt !< Model(baroclinic) time-step [T ~> s].
121  type(meke_cs), pointer :: cs !< MEKE control structure.
122  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: hu !< Accumlated zonal mass flux [H L2 ~> m3 or kg].
123  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: hv !< Accumlated meridional mass flux [H L2 ~> m3 or kg]
124 
125  ! Local variables
126  real, dimension(SZI_(G),SZJ_(G)) :: &
127  mass, & ! The total mass of the water column [R Z ~> kg m-2].
128  i_mass, & ! The inverse of mass [R-1 Z-1 ~> m2 kg-1].
129  src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3).
130  meke_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1].
131  drag_rate_visc, & ! Near-bottom velocity contribution to bottom dratg [L T-1 ~> m s-1]
132  drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
133  drag_rate_j15, & ! The MEKE spindown timescale due to bottom drag with the Jansen 2015 scheme.
134  ! Unfortunately, as written the units seem inconsistent. [T-1 ~> s-1].
135  del2meke, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2].
136  del4meke, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
137  lmixscale, & ! Eddy mixing length [L ~> m].
138  barotrfac2, & ! Ratio of EKE_barotropic / EKE [nondim]
139  bottomfac2, & ! Ratio of EKE_bottom / EKE [nondim]
140  tmp ! Temporary variable for diagnostic computation
141 
142  real, dimension(SZIB_(G),SZJ_(G)) :: &
143  meke_uflux, & ! The zonal advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m-2 s-3].
144  ! In one place, MEKE_uflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
145  kh_u, & ! The zonal diffusivity that is actually used [L2 T-1 ~> m2 s-1].
146  barohu, & ! Depth integrated accumulated zonal mass flux [R Z L2 ~> kg].
147  drag_vel_u ! A (vertical) viscosity associated with bottom drag at
148  ! u-points [Z T-1 ~> m s-1].
149  real, dimension(SZI_(G),SZJB_(G)) :: &
150  meke_vflux, & ! The meridional advective and diffusive flux of MEKE with units of [R Z L4 T-3 ~> kg m-2 s-3].
151  ! In one place, MEKE_vflux is used as temporary work space with units of [L2 T-2 ~> m2 s-2].
152  kh_v, & ! The meridional diffusivity that is actually used [L2 T-1 ~> m2 s-1].
153  barohv, & ! Depth integrated accumulated meridional mass flux [R Z L2 ~> kg].
154  drag_vel_v ! A (vertical) viscosity associated with bottom drag at
155  ! v-points [Z T-1 ~> m s-1].
156  real :: kh_here ! The local horizontal viscosity [L2 T-1 ~> m2 s-1]
157  real :: inv_kh_max ! The inverse of the local horizontal viscosity [T L-2 ~> s m-2]
158  real :: k4_here ! The local horizontal biharmonic viscosity [L4 T-1 ~> m4 s-1]
159  real :: inv_k4_max ! The inverse of the local horizontal biharmonic viscosity [T L-4 ~> s m-4]
160  real :: cdrag2
161  real :: advfac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]
162  real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].
163  real :: ldamping ! The MEKE damping rate [T-1 ~> s-1].
164  real :: rho0 ! A density used to convert mass to distance [R ~> kg m-3].
165  real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate)
166  real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split).
167  logical :: use_drag_rate ! Flag to indicate drag_rate is finite
168  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
169 
170  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
171  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
172 
173  if (.not.associated(cs)) call mom_error(fatal, &
174  "MOM_MEKE: Module must be initialized before it is used.")
175  if (.not.associated(meke)) call mom_error(fatal, &
176  "MOM_MEKE: MEKE must be initialized before it is used.")
177 
178  if ((cs%MEKE_damping > 0.0) .or. (cs%MEKE_Cd_scale > 0.0) .or. (cs%MEKE_Cb>0.) &
179  .or. cs%visc_drag) then
180  use_drag_rate = .true.
181  else
182  use_drag_rate = .false.
183  endif
184 
185  ! Only integrate the MEKE equations if MEKE is required.
186  if (.not.associated(meke%MEKE)) then
187 ! call MOM_error(FATAL, "MOM_MEKE: MEKE%MEKE is not associated!")
188  return
189  endif
190 
191  if (cs%debug) then
192  if (associated(meke%mom_src)) &
193  call hchksum(meke%mom_src, 'MEKE mom_src', g%HI, scale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
194  if (associated(meke%GME_snk)) &
195  call hchksum(meke%GME_snk, 'MEKE GME_snk', g%HI, scale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
196  if (associated(meke%GM_src)) &
197  call hchksum(meke%GM_src, 'MEKE GM_src', g%HI, scale=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
198  if (associated(meke%MEKE)) call hchksum(meke%MEKE, 'MEKE MEKE', g%HI, scale=us%L_T_to_m_s**2)
199  call uvchksum("MEKE SN_[uv]", sn_u, sn_v, g%HI, scale=us%s_to_T, &
200  scalar_pair=.true.)
201  call uvchksum("MEKE h[uv]", hu, hv, g%HI, haloshift=1, &
202  scale=gv%H_to_m*(us%L_to_m**2))
203  endif
204 
205  sdt = dt*cs%MEKE_dtScale ! Scaled dt to use for time-stepping
206  rho0 = gv%Rho0
207  mass_neglect = gv%H_to_RZ * gv%H_subroundoff
208  cdrag2 = cs%cdrag**2
209 
210  ! With a depth-dependent (and possibly strong) damping, it seems
211  ! advisable to use Strang splitting between the damping and diffusion.
212  sdt_damp = sdt ; if (cs%MEKE_KH >= 0.0 .or. cs%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt
213 
214  ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]
215  if (cs%MEKE_advection_factor>0.) then
216  do j=js,je ; do i=is-1,ie
217  barohu(i,j) = 0.
218  enddo ; enddo
219  do k=1,nz
220  do j=js,je ; do i=is-1,ie
221  barohu(i,j) = hu(i,j,k) * gv%H_to_RZ
222  enddo ; enddo
223  enddo
224  do j=js-1,je ; do i=is,ie
225  barohv(i,j) = 0.
226  enddo ; enddo
227  do k=1,nz
228  do j=js-1,je ; do i=is,ie
229  barohv(i,j) = hv(i,j,k) * gv%H_to_RZ
230  enddo ; enddo
231  enddo
232  endif
233 
234  if (cs%MEKE_Cd_scale == 0.0 .and. .not. cs%visc_drag) then
235  !$OMP parallel do default(shared) private(ldamping)
236  do j=js,je ; do i=is,ie
237  drag_rate(i,j) = 0. ; drag_rate_j15(i,j) = 0.
238  enddo ; enddo
239  endif
240 
241  ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
242  if (cs%visc_drag) then
243  !$OMP parallel do default(shared)
244  do j=js,je ; do i=is-1,ie
245  drag_vel_u(i,j) = 0.0
246  if ((g%mask2dCu(i,j) > 0.0) .and. (visc%bbl_thick_u(i,j) > 0.0)) &
247  drag_vel_u(i,j) = visc%Kv_bbl_u(i,j) / visc%bbl_thick_u(i,j)
248  enddo ; enddo
249  !$OMP parallel do default(shared)
250  do j=js-1,je ; do i=is,ie
251  drag_vel_v(i,j) = 0.0
252  if ((g%mask2dCv(i,j) > 0.0) .and. (visc%bbl_thick_v(i,j) > 0.0)) &
253  drag_vel_v(i,j) = visc%Kv_bbl_v(i,j) / visc%bbl_thick_v(i,j)
254  enddo ; enddo
255 
256  !$OMP parallel do default(shared)
257  do j=js,je ; do i=is,ie
258  drag_rate_visc(i,j) = (0.25*g%IareaT(i,j) * us%Z_to_L * &
259  ((g%areaCu(i-1,j)*drag_vel_u(i-1,j) + &
260  g%areaCu(i,j)*drag_vel_u(i,j)) + &
261  (g%areaCv(i,j-1)*drag_vel_v(i,j-1) + &
262  g%areaCv(i,j)*drag_vel_v(i,j)) ) )
263  enddo ; enddo
264  else
265  !$OMP parallel do default(shared)
266  do j=js,je ; do i=is,ie
267  drag_rate_visc(i,j) = 0.
268  enddo ; enddo
269  endif
270 
271  !$OMP parallel do default(shared)
272  do j=js-1,je+1
273  do i=is-1,ie+1 ; mass(i,j) = 0.0 ; enddo
274  do k=1,nz ; do i=is-1,ie+1
275  mass(i,j) = mass(i,j) + g%mask2dT(i,j) * (gv%H_to_RZ * h(i,j,k)) ! [R Z ~> kg m-2]
276  enddo ; enddo
277  do i=is-1,ie+1
278  i_mass(i,j) = 0.0
279  if (mass(i,j) > 0.0) i_mass(i,j) = 1.0 / mass(i,j) ! [R-1 Z-1 ~> m2 kg-1]
280  enddo
281  enddo
282 
283  if (cs%initialize) then
284  call meke_equilibrium(cs, meke, g, gv, us, sn_u, sn_v, drag_rate_visc, i_mass)
285  cs%initialize = .false.
286  endif
287 
288  ! Calculates bottomFac2, barotrFac2 and LmixScale
289  call meke_lengthscales(cs, meke, g, gv, us, sn_u, sn_v, meke%MEKE, bottomfac2, barotrfac2, lmixscale)
290  if (cs%debug) then
291  if (cs%visc_drag) &
292  call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, g%HI, &
293  scale=us%Z_to_m*us%s_to_T, scalar_pair=.true.)
294  call hchksum(mass, 'MEKE mass',g%HI,haloshift=1, scale=us%RZ_to_kg_m2)
295  call hchksum(drag_rate_visc, 'MEKE drag_rate_visc', g%HI, scale=us%L_T_to_m_s)
296  call hchksum(bottomfac2, 'MEKE bottomFac2', g%HI)
297  call hchksum(barotrfac2, 'MEKE barotrFac2', g%HI)
298  call hchksum(lmixscale, 'MEKE LmixScale', g%HI,scale=us%L_to_m)
299  endif
300 
301  ! Aggregate sources of MEKE (background, frictional and GM)
302  !$OMP parallel do default(shared)
303  do j=js,je ; do i=is,ie
304  src(i,j) = cs%MEKE_BGsrc
305  enddo ; enddo
306 
307  if (associated(meke%mom_src)) then
308  !$OMP parallel do default(shared)
309  do j=js,je ; do i=is,ie
310  src(i,j) = src(i,j) - cs%MEKE_FrCoeff*i_mass(i,j)*meke%mom_src(i,j)
311  enddo ; enddo
312  endif
313 
314  if (associated(meke%GME_snk)) then
315  !$OMP parallel do default(shared)
316  do j=js,je ; do i=is,ie
317  src(i,j) = src(i,j) - cs%MEKE_GMECoeff*i_mass(i,j)*meke%GME_snk(i,j)
318  enddo ; enddo
319  endif
320 
321  if (associated(meke%GM_src)) then
322  if (cs%GM_src_alt) then
323  !$OMP parallel do default(shared)
324  do j=js,je ; do i=is,ie
325  src(i,j) = src(i,j) - cs%MEKE_GMcoeff*meke%GM_src(i,j) / &
326  (gv%Rho0 * max(1.0*us%m_to_Z, g%bathyT(i,j)))
327  enddo ; enddo
328  else
329  !$OMP parallel do default(shared)
330  do j=js,je ; do i=is,ie
331  src(i,j) = src(i,j) - cs%MEKE_GMcoeff*i_mass(i,j)*meke%GM_src(i,j)
332  enddo ; enddo
333  endif
334  endif
335 
336  if (cs%MEKE_equilibrium_restoring) then
337  call meke_equilibrium_restoring(cs, g, us, sn_u, sn_v)
338  do j=js,je ; do i=is,ie
339  src(i,j) = src(i,j) - cs%MEKE_restoring_rate*(meke%MEKE(i,j) - cs%equilibrium_value(i,j))
340  enddo ; enddo
341  endif
342 
343  ! Increase EKE by a full time-steps worth of source
344  !$OMP parallel do default(shared)
345  do j=js,je ; do i=is,ie
346  meke%MEKE(i,j) = (meke%MEKE(i,j) + sdt*src(i,j))*g%mask2dT(i,j)
347  enddo ; enddo
348 
349  if (use_drag_rate) then
350  ! Calculate a viscous drag rate (includes BBL contributions from mean flow and eddies)
351  !$OMP parallel do default(shared)
352  do j=js,je ; do i=is,ie
353  drag_rate(i,j) = (us%L_to_Z*rho0 * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
354  cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )
355  enddo ; enddo
356  endif
357 
358  ! First stage of Strang splitting
359  !$OMP parallel do default(shared)
360  do j=js,je ; do i=is,ie
361  ldamping = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)
362  if (meke%MEKE(i,j) < 0.) ldamping = 0.
363  ! notice that the above line ensures a damping only if MEKE is positive,
364  ! while leaving MEKE unchanged if it is negative
365  meke%MEKE(i,j) = meke%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
366  meke_decay(i,j) = ldamping*g%mask2dT(i,j)
367  enddo ; enddo
368 
369  if (cs%kh_flux_enabled .or. cs%MEKE_K4 >= 0.0) then
370  ! Update MEKE in the halos for lateral or bi-harmonic diffusion
371  call cpu_clock_begin(cs%id_clock_pass)
372  call do_group_pass(cs%pass_MEKE, g%Domain)
373  call cpu_clock_end(cs%id_clock_pass)
374  endif
375 
376  if (cs%MEKE_K4 >= 0.0) then
377  ! Calculate Laplacian of MEKE using MEKE_uflux and MEKE_vflux as temporary work space.
378  !$OMP parallel do default(shared)
379  do j=js-1,je+1 ; do i=is-2,ie+1
380  ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
381  meke_uflux(i,j) = ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * g%mask2dCu(i,j)) * &
382  (meke%MEKE(i+1,j) - meke%MEKE(i,j))
383  ! This would have units of [R Z L2 T-2 ~> kg s-2]
384  ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
385  ! ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
386  ! (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
387  enddo ; enddo
388  !$OMP parallel do default(shared)
389  do j=js-2,je+1 ; do i=is-1,ie+1
390  ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
391  meke_vflux(i,j) = ((g%dx_Cv(i,j)*g%IdyCv(i,j)) * g%mask2dCv(i,j)) * &
392  (meke%MEKE(i,j+1) - meke%MEKE(i,j))
393  ! This would have units of [R Z L2 T-2 ~> kg s-2]
394  ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &
395  ! ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
396  ! (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
397  enddo ; enddo
398 
399  !$OMP parallel do default(shared)
400  do j=js-1,je+1 ; do i=is-1,ie+1 ! del2MEKE has units [T-2 ~> s-2].
401  del2meke(i,j) = g%IareaT(i,j) * &
402  ((meke_uflux(i,j) - meke_uflux(i-1,j)) + (meke_vflux(i,j) - meke_vflux(i,j-1)))
403  enddo ; enddo
404 
405  ! Bi-harmonic diffusion of MEKE
406  !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)
407  do j=js,je ; do i=is-1,ie
408  k4_here = cs%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
409  ! Limit Kh to avoid CFL violations.
410  inv_k4_max = 64.0 * sdt * ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &
411  max(g%IareaT(i,j), g%IareaT(i+1,j)))**2
412  if (k4_here*inv_k4_max > 0.3) k4_here = 0.3 / inv_k4_max
413 
414  ! Here the units of MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
415  meke_uflux(i,j) = ((k4_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &
416  ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
417  (del2meke(i+1,j) - del2meke(i,j))
418  enddo ; enddo
419  !$OMP parallel do default(shared) private(K4_here,Inv_K4_max)
420  do j=js-1,je ; do i=is,ie
421  k4_here = cs%MEKE_K4 ! [L4 T-1 ~> m4 s-1]
422  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
423  if (k4_here*inv_k4_max > 0.3) k4_here = 0.3 / inv_k4_max
424 
425  ! Here the units of MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
426  meke_vflux(i,j) = ((k4_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &
427  ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
428  (del2meke(i,j+1) - del2meke(i,j))
429  enddo ; enddo
430  ! Store change in MEKE arising from the bi-harmonic in del4MEKE [L2 T-2 ~> m2 s-2].
431  !$OMP parallel do default(shared)
432  do j=js,je ; do i=is,ie
433  del4meke(i,j) = (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &
434  ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &
435  (meke_vflux(i,j-1) - meke_vflux(i,j)))
436  enddo ; enddo
437  endif !
438 
439  if (cs%kh_flux_enabled) then
440  ! Lateral diffusion of MEKE
441  kh_here = max(0., cs%MEKE_Kh)
442  !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)
443  do j=js,je ; do i=is-1,ie
444  ! Limit Kh to avoid CFL violations.
445  if (associated(meke%Kh)) &
446  kh_here = max(0., cs%MEKE_Kh) + &
447  cs%KhMEKE_Fac*0.5*(meke%Kh(i,j)+meke%Kh(i+1,j))
448  if (associated(meke%Kh_diff)) &
449  kh_here = max(0.,cs%MEKE_Kh) + &
450  cs%KhMEKE_Fac*0.5*(meke%Kh_diff(i,j)+meke%Kh_diff(i+1,j))
451  inv_kh_max = 2.0*sdt * ((g%dy_Cu(i,j)*g%IdxCu(i,j)) * &
452  max(g%IareaT(i,j),g%IareaT(i+1,j)))
453  if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max
454  kh_u(i,j) = kh_here
455 
456  ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
457  meke_uflux(i,j) = ((kh_here * (g%dy_Cu(i,j)*g%IdxCu(i,j))) * &
458  ((2.0*mass(i,j)*mass(i+1,j)) / ((mass(i,j)+mass(i+1,j)) + mass_neglect)) ) * &
459  (meke%MEKE(i,j) - meke%MEKE(i+1,j))
460  enddo ; enddo
461  !$OMP parallel do default(shared) firstprivate(Kh_here) private(Inv_Kh_max)
462  do j=js-1,je ; do i=is,ie
463  if (associated(meke%Kh)) &
464  kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac * 0.5*(meke%Kh(i,j)+meke%Kh(i,j+1))
465  if (associated(meke%Kh_diff)) &
466  kh_here = max(0.,cs%MEKE_Kh) + cs%KhMEKE_Fac * 0.5*(meke%Kh_diff(i,j)+meke%Kh_diff(i,j+1))
467  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)))
468  if (kh_here*inv_kh_max > 0.25) kh_here = 0.25 / inv_kh_max
469  kh_v(i,j) = kh_here
470 
471  ! Here the units of MEKE_uflux and MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
472  meke_vflux(i,j) = ((kh_here * (g%dx_Cv(i,j)*g%IdyCv(i,j))) * &
473  ((2.0*mass(i,j)*mass(i,j+1)) / ((mass(i,j)+mass(i,j+1)) + mass_neglect)) ) * &
474  (meke%MEKE(i,j) - meke%MEKE(i,j+1))
475  enddo ; enddo
476  if (cs%MEKE_advection_factor>0.) then
477  advfac = cs%MEKE_advection_factor / sdt ! [T-1 ~> s-1]
478  !$OMP parallel do default(shared)
479  do j=js,je ; do i=is-1,ie
480  ! Here the units of the quantities added to MEKE_uflux are [R Z L4 T-3 ~> kg m2 s-3].
481  if (barohu(i,j)>0.) then
482  meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i,j)*advfac
483  elseif (barohu(i,j)<0.) then
484  meke_uflux(i,j) = meke_uflux(i,j) + barohu(i,j)*meke%MEKE(i+1,j)*advfac
485  endif
486  enddo ; enddo
487  !$OMP parallel do default(shared)
488  do j=js-1,je ; do i=is,ie
489  ! Here the units of the quantities added to MEKE_vflux are [R Z L4 T-3 ~> kg m2 s-3].
490  if (barohv(i,j)>0.) then
491  meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j)*advfac
492  elseif (barohv(i,j)<0.) then
493  meke_vflux(i,j) = meke_vflux(i,j) + barohv(i,j)*meke%MEKE(i,j+1)*advfac
494  endif
495  enddo ; enddo
496  endif
497 
498  !$OMP parallel do default(shared)
499  do j=js,je ; do i=is,ie
500  meke%MEKE(i,j) = meke%MEKE(i,j) + (sdt*(g%IareaT(i,j)*i_mass(i,j))) * &
501  ((meke_uflux(i-1,j) - meke_uflux(i,j)) + &
502  (meke_vflux(i,j-1) - meke_vflux(i,j)))
503  enddo ; enddo
504  endif ! MEKE_KH>0
505 
506  ! Add on bi-harmonic tendency
507  if (cs%MEKE_K4 >= 0.0) then
508  !$OMP parallel do default(shared)
509  do j=js,je ; do i=is,ie
510  meke%MEKE(i,j) = meke%MEKE(i,j) + del4meke(i,j)
511  enddo ; enddo
512  endif
513 
514  ! Second stage of Strang splitting
515  if (cs%MEKE_KH >= 0.0 .or. cs%MEKE_K4 >= 0.0) then
516  if (sdt>sdt_damp) then
517  ! Recalculate the drag rate, since MEKE has changed.
518  if (use_drag_rate) then
519  !$OMP parallel do default(shared)
520  do j=js,je ; do i=is,ie
521  drag_rate(i,j) = (us%L_to_Z*rho0 * i_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
522  cdrag2 * ( max(0.0, 2.0*bottomfac2(i,j)*meke%MEKE(i,j)) + cs%MEKE_Uscale**2 ) )
523  enddo ; enddo
524  !$OMP parallel do default(shared)
525  do j=js,je ; do i=is,ie
526  ldamping = cs%MEKE_damping + drag_rate(i,j) * bottomfac2(i,j)
527  if (meke%MEKE(i,j) < 0.) ldamping = 0.
528  ! notice that the above line ensures a damping only if MEKE is positive,
529  ! while leaving MEKE unchanged if it is negative
530  meke%MEKE(i,j) = meke%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
531  meke_decay(i,j) = ldamping*g%mask2dT(i,j)
532  enddo ; enddo
533  endif
534  endif
535  endif ! MEKE_KH>=0
536 
537  ! do j=js,je ; do i=is,ie
538  ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0)
539  ! enddo ; enddo
540 
541  call cpu_clock_begin(cs%id_clock_pass)
542  call do_group_pass(cs%pass_MEKE, g%Domain)
543  call cpu_clock_end(cs%id_clock_pass)
544 
545  ! Calculate diffusivity for main model to use
546  if (cs%MEKE_KhCoeff>0.) then
547  if (.not.cs%MEKE_GEOMETRIC) then
548  if (cs%use_old_lscale) then
549  if (cs%Rd_as_max_scale) then
550  !$OMP parallel do default(shared)
551  do j=js,je ; do i=is,ie
552  meke%Kh(i,j) = (cs%MEKE_KhCoeff * &
553  sqrt(2.*max(0.,barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j)) ) * &
554  min(meke%Rd_dx_h(i,j), 1.0)
555  enddo ; enddo
556  else
557  !$OMP parallel do default(shared)
558  do j=js,je ; do i=is,ie
559  meke%Kh(i,j) = cs%MEKE_KhCoeff * &
560  sqrt(2.*max(0., barotrfac2(i,j)*meke%MEKE(i,j))*g%areaT(i,j))
561  enddo ; enddo
562  endif
563  else
564  !$OMP parallel do default(shared)
565  do j=js,je ; do i=is,ie
566  meke%Kh(i,j) = cs%MEKE_KhCoeff * &
567  sqrt(2.*max(0., barotrfac2(i,j)*meke%MEKE(i,j))) * lmixscale(i,j)
568  enddo ; enddo
569  endif
570  endif
571  endif
572 
573  ! Calculate viscosity for the main model to use
574  if (cs%viscosity_coeff_Ku /=0.) then
575  do j=js,je ; do i=is,ie
576  meke%Ku(i,j) = cs%viscosity_coeff_Ku * sqrt(2.*max(0.,meke%MEKE(i,j))) * lmixscale(i,j)
577  enddo ; enddo
578  endif
579 
580  if (cs%viscosity_coeff_Au /=0.) then
581  do j=js,je ; do i=is,ie
582  meke%Au(i,j) = cs%viscosity_coeff_Au * sqrt(2.*max(0.,meke%MEKE(i,j))) * lmixscale(i,j)**3
583  enddo ; enddo
584  endif
585 
586  if (associated(meke%Kh) .or. associated(meke%Ku) .or. associated(meke%Au)) then
587  call cpu_clock_begin(cs%id_clock_pass)
588  call do_group_pass(cs%pass_Kh, g%Domain)
589  call cpu_clock_end(cs%id_clock_pass)
590  endif
591 
592  ! Offer fields for averaging.
593  if (any([cs%id_Ue, cs%id_Ub, cs%id_Ut] > 0)) &
594  tmp(:,:) = 0.
595  if (cs%id_MEKE>0) call post_data(cs%id_MEKE, meke%MEKE, cs%diag)
596  if (cs%id_Ue>0) then
597  do j=js,je ; do i=is,ie
598  tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j)))
599  enddo ; enddo
600  call post_data(cs%id_Ue, tmp, cs%diag)
601  endif
602  if (cs%id_Ub>0) then
603  do j=js,je ; do i=is,ie
604  tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j) * bottomfac2(i,j)))
605  enddo ; enddo
606  call post_data(cs%id_Ub, tmp, cs%diag)
607  endif
608  if (cs%id_Ut>0) then
609  do j=js,je ; do i=is,ie
610  tmp(i,j) = sqrt(max(0., 2. * meke%MEKE(i,j) * barotrfac2(i,j)))
611  enddo ; enddo
612  call post_data(cs%id_Ut, tmp, cs%diag)
613  endif
614  if (cs%id_Kh>0) call post_data(cs%id_Kh, meke%Kh, cs%diag)
615  if (cs%id_Ku>0) call post_data(cs%id_Ku, meke%Ku, cs%diag)
616  if (cs%id_Au>0) call post_data(cs%id_Au, meke%Au, cs%diag)
617  if (cs%id_KhMEKE_u>0) call post_data(cs%id_KhMEKE_u, kh_u, cs%diag)
618  if (cs%id_KhMEKE_v>0) call post_data(cs%id_KhMEKE_v, kh_v, cs%diag)
619  if (cs%id_src>0) call post_data(cs%id_src, src, cs%diag)
620  if (cs%id_decay>0) call post_data(cs%id_decay, meke_decay, cs%diag)
621  if (cs%id_GM_src>0) call post_data(cs%id_GM_src, meke%GM_src, cs%diag)
622  if (cs%id_mom_src>0) call post_data(cs%id_mom_src, meke%mom_src, cs%diag)
623  if (cs%id_GME_snk>0) call post_data(cs%id_GME_snk, meke%GME_snk, cs%diag)
624  if (cs%id_Le>0) call post_data(cs%id_Le, lmixscale, cs%diag)
625  if (cs%id_gamma_b>0) then
626  do j=js,je ; do i=is,ie
627  bottomfac2(i,j) = sqrt(bottomfac2(i,j))
628  enddo ; enddo
629  call post_data(cs%id_gamma_b, bottomfac2, cs%diag)
630  endif
631  if (cs%id_gamma_t>0) then
632  do j=js,je ; do i=is,ie
633  barotrfac2(i,j) = sqrt(barotrfac2(i,j))
634  enddo ; enddo
635  call post_data(cs%id_gamma_t, barotrfac2, cs%diag)
636  endif
637 
638 end subroutine step_forward_meke
639 
640 !> Calculates the equilibrium solutino where the source depends only on MEKE diffusivity
641 !! and there is no lateral diffusion of MEKE.
642 !! Results is in MEKE%MEKE.
643 subroutine meke_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass)
644  type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
645  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure.
646  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
647  type(meke_cs), pointer :: CS !< MEKE control structure.
648  type(meke_type), pointer :: MEKE !< A structure with MEKE data.
649  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
650  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
651  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: drag_rate_visc !< Mean flow velocity contribution
652  !! to the MEKE drag rate [L T-1 ~> m s-1]
653  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_mass !< Inverse of column mass [R-1 Z-1 ~> m2 kg-1].
654  ! Local variables
655  real :: beta ! Combined topograpic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
656  real :: SN ! The local Eady growth rate [T-1 ~> s-1]
657  real :: bottomFac2, barotrFac2 ! Vertical structure factors [nondim]
658  real :: LmixScale, LRhines, LEady ! Various mixing length scales [L ~> m]
659  real :: I_H, KhCoeff
660  real :: Kh ! A lateral diffusivity [L2 T-1 ~> m2 s-1]
661  real :: Ubg2 ! Background (tidal?) velocity squared [L2 T-2 ~> m2 s-2]
662  real :: cd2
663  real :: drag_rate ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
664  real :: src ! The sum of MEKE sources [L2 T-3 ~> W kg-1]
665  real :: ldamping ! The MEKE damping rate [T-1 ~> s-1].
666  real :: EKE, EKEmin, EKEmax, EKEerr ! [L2 T-2 ~> m2 s-2]
667  real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1]
668  real :: FatH ! Coriolis parameter at h points; to compute topographic beta [T-1 ~> s-1]
669  real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
670  integer :: i, j, is, ie, js, je, n1, n2
671  real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2].
672  logical :: useSecant, debugIteration
673 
674  real :: Lgrid, Ldeform, Lfrict
675 
676  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
677 
678  debugiteration = .false.
679  khcoeff = cs%MEKE_KhCoeff
680  ubg2 = cs%MEKE_Uscale**2
681  cd2 = cs%cdrag**2
682  tolerance = 1.0e-12*us%m_s_to_L_T**2
683 
684 !$OMP do
685  do j=js,je ; do i=is,ie
686  ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
687  ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
688  sn = min(sn_u(i,j), sn_u(i-1,j), sn_v(i,j), sn_v(i,j-1))
689 
690  if (cs%MEKE_equilibrium_alt) then
691  meke%MEKE(i,j) = (cs%MEKE_GEOMETRIC_alpha * sn * us%Z_to_m*g%bathyT(i,j))**2 / cd2
692  else
693  fath = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
694  (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1))) ! Coriolis parameter at h points
695 
696  ! Since zero-bathymetry cells are masked, this avoids calculations on land
697  if (cs%MEKE_topographic_beta == 0. .or. g%bathyT(i,j) == 0.) then
698  beta_topo_x = 0. ; beta_topo_y = 0.
699  else
700  !### Consider different combinations of these estimates of topographic beta, and the use
701  ! of the water column thickness instead of the bathymetric depth.
702  beta_topo_x = -cs%MEKE_topographic_beta * fath * 0.5 * ( &
703  (g%bathyT(i+1,j)-g%bathyT(i,j)) * g%IdxCu(i,j) &
704  / max(g%bathyT(i+1,j),g%bathyT(i,j), gv%H_subroundoff) &
705  + (g%bathyT(i,j)-g%bathyT(i-1,j)) * g%IdxCu(i-1,j) &
706  / max(g%bathyT(i,j),g%bathyT(i-1,j), gv%H_subroundoff) )
707  beta_topo_y = -cs%MEKE_topographic_beta * fath * 0.5 * ( &
708  (g%bathyT(i,j+1)-g%bathyT(i,j)) * g%IdyCv(i,j) &
709  / max(g%bathyT(i,j+1),g%bathyT(i,j), gv%H_subroundoff) + &
710  (g%bathyT(i,j)-g%bathyT(i,j-1)) * g%IdyCv(i,j-1) &
711  / max(g%bathyT(i,j),g%bathyT(i,j-1), gv%H_subroundoff) )
712  endif
713  beta = sqrt((g%dF_dx(i,j) + beta_topo_x)**2 + &
714  (g%dF_dy(i,j) + beta_topo_y)**2 )
715 
716  i_h = us%L_to_Z*gv%Rho0 * i_mass(i,j)
717 
718  if (khcoeff*sn*i_h>0.) then
719  ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E
720  ekemin = 0. ! Use the trivial root as the left bracket
721  resmin = 0. ! Need to detect direction of left residual
722  ekemax = 0.01*us%m_s_to_L_T**2 ! First guess at right bracket
723  usesecant = .false. ! Start using a bisection method
724 
725  ! First find right bracket for which resid<0
726  resid = 1.0*us%m_to_L**2*us%T_to_s**3 ; n1 = 0
727  do while (resid>0.)
728  n1 = n1 + 1
729  eke = ekemax
730  call meke_lengthscales_0d(cs, us, g%areaT(i,j), beta, g%bathyT(i,j), &
731  meke%Rd_dx_h(i,j), sn, eke, &
732  bottomfac2, barotrfac2, lmixscale, lrhines, leady)
733  ! TODO: Should include resolution function in Kh
734  kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)
735  src = kh * (sn * sn)
736  drag_rate = i_h * sqrt(drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )
737  ldamping = cs%MEKE_damping + drag_rate * bottomfac2
738  resid = src - ldamping * eke
739  ! if (debugIteration) then
740  ! write(0,*) n1, 'EKE=',EKE,'resid=',resid
741  ! write(0,*) 'EKEmin=',EKEmin,'ResMin=',ResMin
742  ! write(0,*) 'src=',src,'ldamping=',ldamping
743  ! write(0,*) 'gamma-b=',bottomFac2,'gamma-t=',barotrFac2
744  ! write(0,*) 'drag_visc=',drag_rate_visc(i,j),'Ubg2=',Ubg2
745  ! endif
746  if (resid>0.) then ! EKE is to the left of the root
747  ekemin = eke ! so we move the left bracket here
748  ekemax = 10. * eke ! and guess again for the right bracket
749  if (resid<resmin) usesecant = .true.
750  resmin = resid
751  if (ekemax > 2.e17*us%m_s_to_L_T**2) then
752  if (debugiteration) stop 'Something has gone very wrong'
753  debugiteration = .true.
754  resid = 1. ; n1 = 0
755  ekemin = 0. ; resmin = 0.
756  ekemax = 0.01*us%m_s_to_L_T**2
757  usesecant = .false.
758  endif
759  endif
760  enddo ! while(resid>0.) searching for right bracket
761  resmax = resid
762 
763  ! Bisect the bracket
764  n2 = 0 ; ekeerr = ekemax - ekemin
765  do while (ekeerr > tolerance)
766  n2 = n2 + 1
767  if (usesecant) then
768  eke = ekemin + (ekemax - ekemin) * (resmin / (resmin - resmax))
769  else
770  eke = 0.5 * (ekemin + ekemax)
771  endif
772  ekeerr = min( eke-ekemin, ekemax-eke )
773  ! TODO: Should include resolution function in Kh
774  kh = (khcoeff * sqrt(2.*barotrfac2*eke) * lmixscale)
775  src = kh * (sn * sn)
776  drag_rate = i_h * sqrt( drag_rate_visc(i,j)**2 + cd2 * ( 2.0*bottomfac2*eke + ubg2 ) )
777  ldamping = cs%MEKE_damping + drag_rate * bottomfac2
778  resid = src - ldamping * eke
779  if (usesecant .and. resid>resmin) usesecant = .false.
780  if (resid>0.) then ! EKE is to the left of the root
781  ekemin = eke ! so we move the left bracket here
782  if (resid<resmin) usesecant = .true.
783  resmin = resid ! Save this for the secant method
784  elseif (resid<0.) then ! EKE is to the right of the root
785  ekemax = eke ! so we move the right bracket here
786  resmax = resid ! Save this for the secant method
787  else
788  exit ! resid=0 => EKE is exactly at the root
789  endif
790  if (n2>200) stop 'Failing to converge?'
791  enddo ! while(EKEmax-EKEmin>tolerance)
792 
793  else
794  eke = 0.
795  endif
796  meke%MEKE(i,j) = eke
797  endif
798  enddo ; enddo
799 
800 end subroutine meke_equilibrium
801 
802 
803 !< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into
804 !! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value
805 subroutine meke_equilibrium_restoring(CS, G, US, SN_u, SN_v)
806  type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
807  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type.
808  type(meke_cs), pointer :: CS !< MEKE control structure.
809  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
810  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
811  ! Local variables
812  real :: SN ! The local Eady growth rate [T-1 ~> s-1]
813  integer :: i, j, is, ie, js, je ! local indices
814  real :: cd2 ! bottom drag
815 
816  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
817  cd2 = cs%cdrag**2
818 
819  if (.not. associated(cs%equilibrium_value)) allocate(cs%equilibrium_value(szi_(g),szj_(g)))
820  cs%equilibrium_value(:,:) = 0.0
821 
822 !$OMP do
823  do j=js,je ; do i=is,ie
824  ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.)
825  ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v
826  sn = min(sn_u(i,j), sn_u(i-1,j), sn_v(i,j), sn_v(i,j-1))
827  cs%equilibrium_value(i,j) = (cs%MEKE_GEOMETRIC_alpha * sn * us%Z_to_m*g%bathyT(i,j))**2 / cd2
828  enddo ; enddo
829 
830  if (cs%id_MEKE_equilibrium>0) call post_data(cs%id_MEKE_equilibrium, cs%equilibrium_value, cs%diag)
831 
832 end subroutine meke_equilibrium_restoring
833 
834 !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
835 !! functions that are ratios of either bottom or barotropic eddy energy to the
836 !! column eddy energy, respectively. See \ref section_MEKE_equations.
837 subroutine meke_lengthscales(CS, MEKE, G, GV, US, SN_u, SN_v, &
838  EKE, bottomFac2, barotrFac2, LmixScale)
839  type(meke_cs), pointer :: CS !< MEKE control structure.
840  type(meke_type), pointer :: MEKE !< MEKE data.
841  type(ocean_grid_type), intent(inout) :: G !< Ocean grid.
842  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure.
843  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
844  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1].
845  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1].
846  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
847  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomFac2 !< gamma_b^2
848  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrFac2 !< gamma_t^2
849  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
850  ! Local variables
851  real, dimension(SZI_(G),SZJ_(G)) :: LRhines, LEady ! Possible mixing length scales [L ~> m]
852  real :: beta ! Combined topograpic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1]
853  real :: SN ! The local Eady growth rate [T-1 ~> s-1]
854  real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1]
855  real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1]
856  integer :: i, j, is, ie, js, je
857 
858  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
859 
860 !$OMP do
861  do j=js,je ; do i=is,ie
862  if (.not.cs%use_old_lscale) then
863  if (cs%aEady > 0.) then
864  sn = 0.25 * ( (sn_u(i,j) + sn_u(i-1,j)) + (sn_v(i,j) + sn_v(i,j-1)) )
865  else
866  sn = 0.
867  endif
868  fath = 0.25* ( ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1) ) + &
869  ( g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1) ) ) ! Coriolis parameter at h points
870 
871  ! If bathyT is zero, then a division by zero FPE will be raised. In this
872  ! case, we apply Adcroft's rule of reciprocals and set the term to zero.
873  ! Since zero-bathymetry cells are masked, this should not affect values.
874  if (cs%MEKE_topographic_beta == 0. .or. g%bathyT(i,j) == 0.0) then
875  beta_topo_x = 0. ; beta_topo_y = 0.
876  else
877  !### Consider different combinations of these estimates of topographic beta, and the use
878  ! of the water column thickness instead of the bathymetric depth.
879  beta_topo_x = -cs%MEKE_topographic_beta * fath * 0.5 * ( &
880  (g%bathyT(i+1,j)-g%bathyT(i,j)) * g%IdxCu(i,j) &
881  / max(g%bathyT(i+1,j),g%bathyT(i,j), gv%H_subroundoff) &
882  + (g%bathyT(i,j)-g%bathyT(i-1,j)) * g%IdxCu(i-1,j) &
883  / max(g%bathyT(i,j),g%bathyT(i-1,j), gv%H_subroundoff) )
884  beta_topo_y = -cs%MEKE_topographic_beta * fath * 0.5 * ( &
885  (g%bathyT(i,j+1)-g%bathyT(i,j)) * g%IdyCv(i,j) &
886  / max(g%bathyT(i,j+1),g%bathyT(i,j), gv%H_subroundoff) + &
887  (g%bathyT(i,j)-g%bathyT(i,j-1)) * g%IdyCv(i,j-1) &
888  / max(g%bathyT(i,j),g%bathyT(i,j-1), gv%H_subroundoff) )
889  endif
890  beta = sqrt((g%dF_dx(i,j) + beta_topo_x)**2 + &
891  (g%dF_dy(i,j) + beta_topo_y)**2 )
892 
893  else
894  beta = 0.
895  endif
896  ! Returns bottomFac2, barotrFac2 and LmixScale
897  call meke_lengthscales_0d(cs, us, g%areaT(i,j), beta, g%bathyT(i,j), &
898  meke%Rd_dx_h(i,j), sn, meke%MEKE(i,j), &
899  bottomfac2(i,j), barotrfac2(i,j), lmixscale(i,j), &
900  lrhines(i,j), leady(i,j))
901  enddo ; enddo
902  if (cs%id_Lrhines>0) call post_data(cs%id_LRhines, lrhines, cs%diag)
903  if (cs%id_Leady>0) call post_data(cs%id_LEady, leady, cs%diag)
904 
905 end subroutine meke_lengthscales
906 
907 !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$
908 !! functions that are ratios of either bottom or barotropic eddy energy to the
909 !! column eddy energy, respectively. See \ref section_MEKE_equations.
910 subroutine meke_lengthscales_0d(CS, US, area, beta, depth, Rd_dx, SN, EKE, & ! Z_to_L, &
911  bottomFac2, barotrFac2, LmixScale, Lrhines, Leady)
912  type(meke_cs), pointer :: CS !< MEKE control structure.
913  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
914  real, intent(in) :: area !< Grid cell area [L2 ~> m2]
915  real, intent(in) :: beta !< Planetary beta = \f$ \nabla f\f$ [T-1 L-1 ~> s-1 m-1]
916  real, intent(in) :: depth !< Ocean depth [Z ~> m]
917  real, intent(in) :: Rd_dx !< Resolution Ld/dx [nondim].
918  real, intent(in) :: SN !< Eady growth rate [T-1 ~> s-1].
919  real, intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2].
920 ! real, intent(in) :: Z_to_L !< A conversion factor from depth units (Z) to
921 ! !! the units for lateral distances (L).
922  real, intent(out) :: bottomFac2 !< gamma_b^2
923  real, intent(out) :: barotrFac2 !< gamma_t^2
924  real, intent(out) :: LmixScale !< Eddy mixing length [L ~> m].
925  real, intent(out) :: Lrhines !< Rhines length scale [L ~> m].
926  real, intent(out) :: Leady !< Eady length scale [L ~> m].
927  ! Local variables
928  real :: Lgrid, Ldeform, Lfrict ! Length scales [L ~> m]
929  real :: Ue ! An eddy velocity [L T-1 ~> m s-1]
930 
931  ! Length scale for MEKE derived diffusivity
932  lgrid = sqrt(area) ! Grid scale
933  ldeform = lgrid * rd_dx ! Deformation scale
934  lfrict = (us%Z_to_L * depth) / cs%cdrag ! Frictional arrest scale
935  ! gamma_b^2 is the ratio of bottom eddy energy to mean column eddy energy
936  ! used in calculating bottom drag
937  bottomfac2 = cs%MEKE_CD_SCALE**2
938  if (lfrict*cs%MEKE_Cb>0.) bottomfac2 = bottomfac2 + 1./( 1. + cs%MEKE_Cb*(ldeform/lfrict) )**0.8
939  bottomfac2 = max(bottomfac2, cs%MEKE_min_gamma)
940  ! gamma_t^2 is the ratio of barotropic eddy energy to mean column eddy energy
941  ! used in the velocity scale for diffusivity
942  barotrfac2 = 1.
943  if (lfrict*cs%MEKE_Ct>0.) barotrfac2 = 1. / ( 1. + cs%MEKE_Ct*(ldeform/lfrict) )**0.25
944  barotrfac2 = max(barotrfac2, cs%MEKE_min_gamma)
945  if (cs%use_old_lscale) then
946  if (cs%Rd_as_max_scale) then
947  lmixscale = min(ldeform, lgrid) ! The smaller of Ld or dx
948  else
949  lmixscale = lgrid
950  endif
951  else
952  ue = sqrt( 2.0 * max( 0., barotrfac2*eke ) ) ! Barotropic eddy flow scale
953  lrhines = sqrt( ue / max( beta, 1.e-30*us%T_to_s*us%L_to_m ) ) ! Rhines scale
954  if (cs%aEady > 0.) then
955  leady = ue / max( sn, 1.e-15*us%T_to_s ) ! Bound Eady time-scale < 1e15 seconds
956  else
957  leady = 0.
958  endif
959  if (cs%use_min_lscale) then
960  lmixscale = 1.e7
961  if (cs%aDeform*ldeform > 0.) lmixscale = min(lmixscale,cs%aDeform*ldeform)
962  if (cs%aFrict *lfrict > 0.) lmixscale = min(lmixscale,cs%aFrict *lfrict)
963  if (cs%aRhines*lrhines > 0.) lmixscale = min(lmixscale,cs%aRhines*lrhines)
964  if (cs%aEady *leady > 0.) lmixscale = min(lmixscale,cs%aEady *leady)
965  if (cs%aGrid *lgrid > 0.) lmixscale = min(lmixscale,cs%aGrid *lgrid)
966  if (cs%Lfixed > 0.) lmixscale = min(lmixscale,cs%Lfixed)
967  else
968  lmixscale = 0.
969  if (cs%aDeform*ldeform > 0.) lmixscale = lmixscale + 1./(cs%aDeform*ldeform)
970  if (cs%aFrict *lfrict > 0.) lmixscale = lmixscale + 1./(cs%aFrict *lfrict)
971  if (cs%aRhines*lrhines > 0.) lmixscale = lmixscale + 1./(cs%aRhines*lrhines)
972  if (cs%aEady *leady > 0.) lmixscale = lmixscale + 1./(cs%aEady *leady)
973  if (cs%aGrid *lgrid > 0.) lmixscale = lmixscale + 1./(cs%aGrid *lgrid)
974  if (cs%Lfixed > 0.) lmixscale = lmixscale + 1./cs%Lfixed
975  if (lmixscale > 0.) lmixscale = 1. / lmixscale
976  endif
977  endif
978 
979 end subroutine meke_lengthscales_0d
980 
981 !> Initializes the MOM_MEKE module and reads parameters.
982 !! Returns True if module is to be used, otherwise returns False.
983 logical function meke_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS)
984  type(time_type), intent(in) :: time !< The current model time.
985  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
986  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
987  type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
988  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
989  type(meke_cs), pointer :: cs !< MEKE control structure.
990  type(meke_type), pointer :: meke !< MEKE-related fields.
991  type(mom_restart_cs), pointer :: restart_cs !< Restart control structure for MOM_MEKE.
992 
993  ! Local variables
994  real :: i_t_rescale ! A rescaling factor for time from the internal representation in this
995  ! run to the representation in a restart file.
996  real :: l_rescale ! A rescaling factor for length from the internal representation in this
997  ! run to the representation in a restart file.
998  real :: meke_restoring_timescale ! The timescale used to nudge MEKE toward its equilibrium value.
999  real :: cdrag ! The default bottom drag coefficient [nondim].
1000  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1001  logical :: laplacian, biharmonic, usevarmix, coldstart
1002  ! This include declares and sets the variable "version".
1003 # include "version_variable.h"
1004  character(len=40) :: mdl = "MOM_MEKE" ! This module's name.
1005 
1006  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1007  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1008 
1009  ! Determine whether this module will be used
1010  call get_param(param_file, mdl, "USE_MEKE", meke_init, default=.false., do_not_log=.true.)
1011  call log_version(param_file, mdl, version, "", all_default=.not.meke_init)
1012  call get_param(param_file, mdl, "USE_MEKE", meke_init, &
1013  "If true, turns on the MEKE scheme which calculates "// &
1014  "a sub-grid mesoscale eddy kinetic energy budget.", &
1015  default=.false.)
1016  if (.not. meke_init) return
1017 
1018  if (.not. associated(meke)) then
1019  ! The MEKE structure should have been allocated in MEKE_alloc_register_restart()
1020  call mom_error(warning, "MEKE_init called with NO associated "// &
1021  "MEKE-type structure.")
1022  return
1023  endif
1024  if (associated(cs)) then
1025  call mom_error(warning, &
1026  "MEKE_init called with an associated control structure.")
1027  return
1028  else ; allocate(cs) ; endif
1029 
1030  call mom_mesg("MEKE_init: reading parameters ", 5)
1031 
1032  ! Read all relevant parameters and write them to the model log.
1033  call get_param(param_file, mdl, "MEKE_DAMPING", cs%MEKE_damping, &
1034  "The local depth-independent MEKE dissipation rate.", &
1035  units="s-1", default=0.0, scale=us%T_to_s)
1036  call get_param(param_file, mdl, "MEKE_CD_SCALE", cs%MEKE_Cd_scale, &
1037  "The ratio of the bottom eddy velocity to the column mean "//&
1038  "eddy velocity, i.e. sqrt(2*MEKE). This should be less than 1 "//&
1039  "to account for the surface intensification of MEKE.", &
1040  units="nondim", default=0.)
1041  call get_param(param_file, mdl, "MEKE_CB", cs%MEKE_Cb, &
1042  "A coefficient in the expression for the ratio of bottom projected "//&
1043  "eddy energy and mean column energy (see Jansen et al. 2015).",&
1044  units="nondim", default=25.)
1045  call get_param(param_file, mdl, "MEKE_MIN_GAMMA2", cs%MEKE_min_gamma, &
1046  "The minimum allowed value of gamma_b^2.",&
1047  units="nondim", default=0.0001)
1048  call get_param(param_file, mdl, "MEKE_CT", cs%MEKE_Ct, &
1049  "A coefficient in the expression for the ratio of barotropic "//&
1050  "eddy energy and mean column energy (see Jansen et al. 2015).",&
1051  units="nondim", default=50.)
1052  call get_param(param_file, mdl, "MEKE_GMCOEFF", cs%MEKE_GMcoeff, &
1053  "The efficiency of the conversion of potential energy "//&
1054  "into MEKE by the thickness mixing parameterization. "//&
1055  "If MEKE_GMCOEFF is negative, this conversion is not "//&
1056  "used or calculated.", units="nondim", default=-1.0)
1057  call get_param(param_file, mdl, "MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &
1058  "If MEKE_GEOMETRIC is true, uses the GM coefficient formulation "//&
1059  "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
1060  call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &
1061  "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
1062  "thickness diffusion.", units="nondim", default=0.05)
1063  call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_ALT", cs%MEKE_equilibrium_alt, &
1064  "If true, use an alternative formula for computing the (equilibrium)"//&
1065  "initial value of MEKE.", default=.false.)
1066  call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", cs%MEKE_equilibrium_restoring, &
1067  "If true, restore MEKE back to its equilibrium value, which is calculated at "//&
1068  "each time step.", default=.false.)
1069  if (cs%MEKE_equilibrium_restoring) then
1070  call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", meke_restoring_timescale, &
1071  "The timescale used to nudge MEKE toward its equilibrium value.", units="s", &
1072  default=1e6, scale=us%T_to_s)
1073  cs%MEKE_restoring_rate = 1.0 / meke_restoring_timescale
1074  endif
1075 
1076  call get_param(param_file, mdl, "MEKE_FRCOEFF", cs%MEKE_FrCoeff, &
1077  "The efficiency of the conversion of mean energy into "//&
1078  "MEKE. If MEKE_FRCOEFF is negative, this conversion "//&
1079  "is not used or calculated.", units="nondim", default=-1.0)
1080  call get_param(param_file, mdl, "MEKE_GMECOEFF", cs%MEKE_GMECoeff, &
1081  "The efficiency of the conversion of MEKE into mean energy "//&
1082  "by GME. If MEKE_GMECOEFF is negative, this conversion "//&
1083  "is not used or calculated.", units="nondim", default=-1.0)
1084  call get_param(param_file, mdl, "MEKE_BGSRC", cs%MEKE_BGsrc, &
1085  "A background energy source for MEKE.", units="W kg-1", &
1086  default=0.0, scale=us%m_to_L**2*us%T_to_s**3)
1087  call get_param(param_file, mdl, "MEKE_KH", cs%MEKE_Kh, &
1088  "A background lateral diffusivity of MEKE. "//&
1089  "Use a negative value to not apply lateral diffusion to MEKE.", &
1090  units="m2 s-1", default=-1.0, scale=us%m_to_L**2*us%T_to_s)
1091  call get_param(param_file, mdl, "MEKE_K4", cs%MEKE_K4, &
1092  "A lateral bi-harmonic diffusivity of MEKE. "//&
1093  "Use a negative value to not apply bi-harmonic diffusion to MEKE.", &
1094  units="m4 s-1", default=-1.0, scale=us%m_to_L**4*us%T_to_s)
1095  call get_param(param_file, mdl, "MEKE_DTSCALE", cs%MEKE_dtScale, &
1096  "A scaling factor to accelerate the time evolution of MEKE.", &
1097  units="nondim", default=1.0)
1098  call get_param(param_file, mdl, "MEKE_KHCOEFF", cs%MEKE_KhCoeff, &
1099  "A scaling factor in the expression for eddy diffusivity "//&
1100  "which is otherwise proportional to the MEKE velocity- "//&
1101  "scale times an eddy mixing-length. This factor "//&
1102  "must be >0 for MEKE to contribute to the thickness/ "//&
1103  "and tracer diffusivity in the rest of the model.", &
1104  units="nondim", default=1.0)
1105  call get_param(param_file, mdl, "MEKE_USCALE", cs%MEKE_Uscale, &
1106  "The background velocity that is combined with MEKE to "//&
1107  "calculate the bottom drag.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
1108  call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", cs%GM_src_alt, &
1109  "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&
1110  "than the streamfunction for the MEKE GM source term.", default=.false.)
1111  call get_param(param_file, mdl, "MEKE_VISC_DRAG", cs%visc_drag, &
1112  "If true, use the vertvisc_type to calculate the bottom "//&
1113  "drag acting on MEKE.", default=.true.)
1114  call get_param(param_file, mdl, "MEKE_KHTH_FAC", meke%KhTh_fac, &
1115  "A factor that maps MEKE%Kh to KhTh.", units="nondim", &
1116  default=0.0)
1117  call get_param(param_file, mdl, "MEKE_KHTR_FAC", meke%KhTr_fac, &
1118  "A factor that maps MEKE%Kh to KhTr.", units="nondim", &
1119  default=0.0)
1120  call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", cs%KhMEKE_Fac, &
1121  "A factor that maps MEKE%Kh to Kh for MEKE itself.", &
1122  units="nondim", default=0.0)
1123  call get_param(param_file, mdl, "MEKE_OLD_LSCALE", cs%use_old_lscale, &
1124  "If true, use the old formula for length scale which is "//&
1125  "a function of grid spacing and deformation radius.", &
1126  default=.false.)
1127  call get_param(param_file, mdl, "MEKE_MIN_LSCALE", cs%use_min_lscale, &
1128  "If true, use a strict minimum of provided length scales "//&
1129  "rather than harmonic mean.", &
1130  default=.false.)
1131  call get_param(param_file, mdl, "MEKE_RD_MAX_SCALE", cs%Rd_as_max_scale, &
1132  "If true, the length scale used by MEKE is the minimum of "//&
1133  "the deformation radius or grid-spacing. Only used if "//&
1134  "MEKE_OLD_LSCALE=True", units="nondim", default=.false.)
1135  call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_KU", cs%viscosity_coeff_Ku, &
1136  "If non-zero, is the scaling coefficient in the expression for"//&
1137  "viscosity used to parameterize harmonic lateral momentum mixing by"//&
1138  "unresolved eddies represented by MEKE. Can be negative to"//&
1139  "represent backscatter from the unresolved eddies.", &
1140  units="nondim", default=0.0)
1141  call get_param(param_file, mdl, "MEKE_VISCOSITY_COEFF_AU", cs%viscosity_coeff_Au, &
1142  "If non-zero, is the scaling coefficient in the expression for"//&
1143  "viscosity used to parameterize biharmonic lateral momentum mixing by"//&
1144  "unresolved eddies represented by MEKE. Can be negative to"//&
1145  "represent backscatter from the unresolved eddies.", &
1146  units="nondim", default=0.0)
1147  call get_param(param_file, mdl, "MEKE_FIXED_MIXING_LENGTH", cs%Lfixed, &
1148  "If positive, is a fixed length contribution to the expression "//&
1149  "for mixing length used in MEKE-derived diffusivity.", &
1150  units="m", default=0.0, scale=us%m_to_L)
1151  call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", cs%aDeform, &
1152  "If positive, is a coefficient weighting the deformation scale "//&
1153  "in the expression for mixing length used in MEKE-derived diffusivity.", &
1154  units="nondim", default=0.0)
1155  call get_param(param_file, mdl, "MEKE_ALPHA_RHINES", cs%aRhines, &
1156  "If positive, is a coefficient weighting the Rhines scale "//&
1157  "in the expression for mixing length used in MEKE-derived diffusivity.", &
1158  units="nondim", default=0.0)
1159  call get_param(param_file, mdl, "MEKE_ALPHA_EADY", cs%aEady, &
1160  "If positive, is a coefficient weighting the Eady length scale "//&
1161  "in the expression for mixing length used in MEKE-derived diffusivity.", &
1162  units="nondim", default=0.0)
1163  call get_param(param_file, mdl, "MEKE_ALPHA_FRICT", cs%aFrict, &
1164  "If positive, is a coefficient weighting the frictional arrest scale "//&
1165  "in the expression for mixing length used in MEKE-derived diffusivity.", &
1166  units="nondim", default=0.0)
1167  call get_param(param_file, mdl, "MEKE_ALPHA_GRID", cs%aGrid, &
1168  "If positive, is a coefficient weighting the grid-spacing as a scale "//&
1169  "in the expression for mixing length used in MEKE-derived diffusivity.", &
1170  units="nondim", default=0.0)
1171  call get_param(param_file, mdl, "MEKE_COLD_START", coldstart, &
1172  "If true, initialize EKE to zero. Otherwise a local equilibrium solution "//&
1173  "is used as an initial condition for EKE.", default=.false.)
1174  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", meke%backscatter_Ro_c, &
1175  "The coefficient in the Rossby number function for scaling the biharmonic "//&
1176  "frictional energy source. Setting to non-zero enables the Rossby number function.", &
1177  units="nondim", default=0.0)
1178  call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_POW", meke%backscatter_Ro_pow, &
1179  "The power in the Rossby number function for scaling the biharmonic "//&
1180  "frictional energy source.", units="nondim", default=0.0)
1181  call get_param(param_file, mdl, "MEKE_ADVECTION_FACTOR", cs%MEKE_advection_factor, &
1182  "A scale factor in front of advection of eddy energy. Zero turns advection off. "//&
1183  "Using unity would be normal but other values could accommodate a mismatch "//&
1184  "between the advecting barotropic flow and the vertical structure of MEKE.", &
1185  units="nondim", default=0.0)
1186  call get_param(param_file, mdl, "MEKE_TOPOGRAPHIC_BETA", cs%MEKE_topographic_beta, &
1187  "A scale factor to determine how much topographic beta is weighed in " //&
1188  "computing beta in the expression of Rhines scale. Use 1 if full "//&
1189  "topographic beta effect is considered; use 0 if it's completely ignored.", &
1190  units="nondim", default=0.0)
1191 
1192  ! Nonlocal module parameters
1193  call get_param(param_file, mdl, "CDRAG", cdrag, &
1194  "CDRAG is the drag coefficient relating the magnitude of "//&
1195  "the velocity field to the bottom stress.", units="nondim", &
1196  default=0.003)
1197  call get_param(param_file, mdl, "MEKE_CDRAG", cs%cdrag, &
1198  "Drag coefficient relating the magnitude of the velocity "//&
1199  "field to the bottom stress in MEKE.", units="nondim", &
1200  default=cdrag)
1201  call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
1202  call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)
1203 
1204  if (cs%viscosity_coeff_Ku/=0. .and. .not. laplacian) call mom_error(fatal, &
1205  "LAPLACIAN must be true if MEKE_VISCOSITY_COEFF_KU is true.")
1206 
1207  if (cs%viscosity_coeff_Au/=0. .and. .not. biharmonic) call mom_error(fatal, &
1208  "BIHARMONIC must be true if MEKE_VISCOSITY_COEFF_AU is true.")
1209 
1210  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
1211 
1212  ! Identify if any lateral diffusive processes are active
1213  cs%kh_flux_enabled = .false.
1214  if ((cs%MEKE_KH >= 0.0) .or. (cs%KhMEKE_FAC > 0.0) .or. (cs%MEKE_advection_factor > 0.0)) &
1215  cs%kh_flux_enabled = .true.
1216 
1217 ! Register fields for output from this module.
1218  cs%diag => diag
1219  cs%id_MEKE = register_diag_field('ocean_model', 'MEKE', diag%axesT1, time, &
1220  'Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=us%L_T_to_m_s**2)
1221  if (.not. associated(meke%MEKE)) cs%id_MEKE = -1
1222  cs%id_Kh = register_diag_field('ocean_model', 'MEKE_KH', diag%axesT1, time, &
1223  'MEKE derived diffusivity', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1224  if (.not. associated(meke%Kh)) cs%id_Kh = -1
1225  cs%id_Ku = register_diag_field('ocean_model', 'MEKE_KU', diag%axesT1, time, &
1226  'MEKE derived lateral viscosity', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1227  if (.not. associated(meke%Ku)) cs%id_Ku = -1
1228  cs%id_Au = register_diag_field('ocean_model', 'MEKE_AU', diag%axesT1, time, &
1229  'MEKE derived lateral biharmonic viscosity', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)
1230  if (.not. associated(meke%Au)) cs%id_Au = -1
1231  cs%id_Ue = register_diag_field('ocean_model', 'MEKE_Ue', diag%axesT1, time, &
1232  'MEKE derived eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)
1233  if (.not. associated(meke%MEKE)) cs%id_Ue = -1
1234  cs%id_Ub = register_diag_field('ocean_model', 'MEKE_Ub', diag%axesT1, time, &
1235  'MEKE derived bottom eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)
1236  if (.not. associated(meke%MEKE)) cs%id_Ub = -1
1237  cs%id_Ut = register_diag_field('ocean_model', 'MEKE_Ut', diag%axesT1, time, &
1238  'MEKE derived barotropic eddy-velocity scale', 'm s-1', conversion=us%L_T_to_m_s)
1239  if (.not. associated(meke%MEKE)) cs%id_Ut = -1
1240  cs%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, time, &
1241  'MEKE energy source', 'm2 s-3', conversion=(us%L_T_to_m_s**2)*us%s_to_T)
1242  cs%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, time, &
1243  'MEKE decay rate', 's-1', conversion=us%s_to_T)
1244  cs%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, time, &
1245  'MEKE energy available from thickness mixing', &
1246  'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
1247  if (.not. associated(meke%GM_src)) cs%id_GM_src = -1
1248  cs%id_mom_src = register_diag_field('ocean_model', 'MEKE_mom_src',diag%axesT1, time, &
1249  'MEKE energy available from momentum', &
1250  'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
1251  if (.not. associated(meke%mom_src)) cs%id_mom_src = -1
1252  cs%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, time, &
1253  'MEKE energy lost to GME backscatter', &
1254  'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
1255  if (.not. associated(meke%GME_snk)) cs%id_GME_snk = -1
1256  cs%id_Le = register_diag_field('ocean_model', 'MEKE_Le', diag%axesT1, time, &
1257  'Eddy mixing length used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)
1258  cs%id_Lrhines = register_diag_field('ocean_model', 'MEKE_Lrhines', diag%axesT1, time, &
1259  'Rhines length scale used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)
1260  cs%id_Leady = register_diag_field('ocean_model', 'MEKE_Leady', diag%axesT1, time, &
1261  'Eady length scale used in the MEKE derived eddy diffusivity', 'm', conversion=us%L_to_m)
1262  cs%id_gamma_b = register_diag_field('ocean_model', 'MEKE_gamma_b', diag%axesT1, time, &
1263  'Ratio of bottom-projected eddy velocity to column-mean eddy velocity', 'nondim')
1264  cs%id_gamma_t = register_diag_field('ocean_model', 'MEKE_gamma_t', diag%axesT1, time, &
1265  'Ratio of barotropic eddy velocity to column-mean eddy velocity', 'nondim')
1266 
1267  if (cs%kh_flux_enabled) then
1268  cs%id_KhMEKE_u = register_diag_field('ocean_model', 'KHMEKE_u', diag%axesCu1, time, &
1269  'Zonal diffusivity of MEKE', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1270  cs%id_KhMEKE_v = register_diag_field('ocean_model', 'KHMEKE_v', diag%axesCv1, time, &
1271  'Meridional diffusivity of MEKE', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1272  endif
1273 
1274  if (cs%MEKE_equilibrium_restoring) then
1275  cs%id_MEKE_equilibrium = register_diag_field('ocean_model', 'MEKE_equilibrium', diag%axesT1, time, &
1276  'Equilibrated Mesoscale Eddy Kinetic Energy', 'm2 s-2', conversion=us%L_T_to_m_s**2)
1277  endif
1278 
1279  cs%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=clock_routine)
1280 
1281  ! Detect whether this instance of MEKE_init() is at the beginning of a run
1282  ! or after a restart. If at the beginning, we will initialize MEKE to a local
1283  ! equilibrium.
1284  cs%initialize = .not.query_initialized(meke%MEKE, "MEKE", restart_cs)
1285  if (coldstart) cs%initialize = .false.
1286  if (cs%initialize) call mom_error(warning, &
1287  "MEKE_init: Initializing MEKE with a local equilibrium balance.")
1288 
1289  ! Account for possible changes in dimensional scaling for variables that have been
1290  ! read from a restart file.
1291  i_t_rescale = 1.0
1292  if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
1293  i_t_rescale = us%s_to_T_restart / us%s_to_T
1294  l_rescale = 1.0
1295  if ((us%m_to_L_restart /= 0.0) .and. (us%m_to_L_restart /= us%m_to_L)) &
1296  l_rescale = us%m_to_L / us%m_to_L_restart
1297 
1298  if (l_rescale*i_t_rescale /= 1.0) then
1299  if (associated(meke%MEKE)) then ; if (query_initialized(meke%MEKE, "MEKE_MEKE", restart_cs)) then
1300  do j=js,je ; do i=is,ie
1301  meke%MEKE(i,j) = l_rescale*i_t_rescale * meke%MEKE(i,j)
1302  enddo ; enddo
1303  endif ; endif
1304  endif
1305  if (l_rescale**2*i_t_rescale /= 1.0) then
1306  if (associated(meke%Kh)) then ; if (query_initialized(meke%Kh, "MEKE_Kh", restart_cs)) then
1307  do j=js,je ; do i=is,ie
1308  meke%Kh(i,j) = l_rescale**2*i_t_rescale * meke%Kh(i,j)
1309  enddo ; enddo
1310  endif ; endif
1311  if (associated(meke%Ku)) then ; if (query_initialized(meke%Ku, "MEKE_Ku", restart_cs)) then
1312  do j=js,je ; do i=is,ie
1313  meke%Ku(i,j) = l_rescale**2*i_t_rescale * meke%Ku(i,j)
1314  enddo ; enddo
1315  endif ; endif
1316  if (associated(meke%Kh_diff)) then ; if (query_initialized(meke%Kh, "MEKE_Kh_diff", restart_cs)) then
1317  do j=js,je ; do i=is,ie
1318  meke%Kh_diff(i,j) = l_rescale**2*i_t_rescale * meke%Kh_diff(i,j)
1319  enddo ; enddo
1320  endif ; endif
1321  endif
1322  if (l_rescale**4*i_t_rescale /= 1.0) then
1323  if (associated(meke%Au)) then ; if (query_initialized(meke%Au, "MEKE_Au", restart_cs)) then
1324  do j=js,je ; do i=is,ie
1325  meke%Au(i,j) = l_rescale**4*i_t_rescale * meke%Au(i,j)
1326  enddo ; enddo
1327  endif ; endif
1328  endif
1329 
1330  ! Set up group passes. In the case of a restart, these fields need a halo update now.
1331  if (associated(meke%MEKE)) then
1332  call create_group_pass(cs%pass_MEKE, meke%MEKE, g%Domain)
1333  if (associated(meke%Kh_diff)) call create_group_pass(cs%pass_MEKE, meke%Kh_diff, g%Domain)
1334  if (.not.cs%initialize) call do_group_pass(cs%pass_MEKE, g%Domain)
1335  endif
1336  if (associated(meke%Kh)) call create_group_pass(cs%pass_Kh, meke%Kh, g%Domain)
1337  if (associated(meke%Ku)) call create_group_pass(cs%pass_Kh, meke%Ku, g%Domain)
1338  if (associated(meke%Au)) call create_group_pass(cs%pass_Kh, meke%Au, g%Domain)
1339 
1340  if (associated(meke%Kh) .or. associated(meke%Ku) .or. associated(meke%Au)) &
1341  call do_group_pass(cs%pass_Kh, g%Domain)
1342 
1343 end function meke_init
1344 
1345 !> Allocates memory and register restart fields for the MOM_MEKE module.
1346 subroutine meke_alloc_register_restart(HI, param_file, MEKE, restart_CS)
1347 ! Arguments
1348  type(hor_index_type), intent(in) :: hi !< Horizontal index structure
1349  type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
1350  type(meke_type), pointer :: meke !< A structure with MEKE-related fields.
1351  type(mom_restart_cs), pointer :: restart_cs !< Restart control structure for MOM_MEKE.
1352 ! Local variables
1353  type(vardesc) :: vd
1354  real :: meke_gmcoeff, meke_frcoeff, meke_gmecoeff, meke_khcoeff, meke_visccoeff_ku, meke_visccoeff_au
1355  logical :: use_kh_in_meke
1356  logical :: usemeke
1357  integer :: isd, ied, jsd, jed
1358 
1359 ! Determine whether this module will be used
1360  usemeke = .false.; call read_param(param_file,"USE_MEKE",usemeke)
1361 
1362 ! Read these parameters to determine what should be in the restarts
1363  meke_gmcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",meke_gmcoeff)
1364  meke_frcoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",meke_frcoeff)
1365  meke_gmecoeff =-1.; call read_param(param_file,"MEKE_GMECOEFF",meke_gmecoeff)
1366  meke_khcoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",meke_khcoeff)
1367  meke_visccoeff_ku =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",meke_visccoeff_ku)
1368  meke_visccoeff_au =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",meke_visccoeff_au)
1369  use_kh_in_meke = .false.; call read_param(param_file,"USE_KH_IN_MEKE", use_kh_in_meke)
1370 ! Allocate control structure
1371  if (associated(meke)) then
1372  call mom_error(warning, "MEKE_alloc_register_restart called with an associated "// &
1373  "MEKE type.")
1374  return
1375  else; allocate(meke); endif
1376 
1377  if (.not. usemeke) return
1378 
1379 ! Allocate memory
1380  call mom_mesg("MEKE_alloc_register_restart: allocating and registering", 5)
1381  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
1382  allocate(meke%MEKE(isd:ied,jsd:jed)) ; meke%MEKE(:,:) = 0.0
1383  vd = var_desc("MEKE", "m2 s-2", hor_grid='h', z_grid='1', &
1384  longname="Mesoscale Eddy Kinetic Energy")
1385  call register_restart_field(meke%MEKE, vd, .false., restart_cs)
1386  if (meke_gmcoeff>=0.) then
1387  allocate(meke%GM_src(isd:ied,jsd:jed)) ; meke%GM_src(:,:) = 0.0
1388  endif
1389  if (meke_frcoeff>=0. .or. meke_gmecoeff>=0.) then
1390  allocate(meke%mom_src(isd:ied,jsd:jed)) ; meke%mom_src(:,:) = 0.0
1391  endif
1392  if (meke_gmecoeff>=0.) then
1393  allocate(meke%GME_snk(isd:ied,jsd:jed)) ; meke%GME_snk(:,:) = 0.0
1394  endif
1395  if (meke_khcoeff>=0.) then
1396  allocate(meke%Kh(isd:ied,jsd:jed)) ; meke%Kh(:,:) = 0.0
1397  vd = var_desc("MEKE_Kh", "m2 s-1", hor_grid='h', z_grid='1', &
1398  longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy")
1399  call register_restart_field(meke%Kh, vd, .false., restart_cs)
1400  endif
1401  allocate(meke%Rd_dx_h(isd:ied,jsd:jed)) ; meke%Rd_dx_h(:,:) = 0.0
1402  if (meke_visccoeff_ku/=0.) then
1403  allocate(meke%Ku(isd:ied,jsd:jed)) ; meke%Ku(:,:) = 0.0
1404  vd = var_desc("MEKE_Ku", "m2 s-1", hor_grid='h', z_grid='1', &
1405  longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy")
1406  call register_restart_field(meke%Ku, vd, .false., restart_cs)
1407  endif
1408  if (use_kh_in_meke) then
1409  allocate(meke%Kh_diff(isd:ied,jsd:jed)) ; meke%Kh_diff(:,:) = 0.0
1410  vd = var_desc("MEKE_Kh_diff", "m2 s-1",hor_grid='h',z_grid='1', &
1411  longname="Copy of thickness diffusivity for diffusing MEKE")
1412  call register_restart_field(meke%Kh_diff, vd, .false., restart_cs)
1413  endif
1414 
1415  if (meke_visccoeff_au/=0.) then
1416  allocate(meke%Au(isd:ied,jsd:jed)) ; meke%Au(:,:) = 0.0
1417  vd = var_desc("MEKE_Au", "m4 s-1", hor_grid='h', z_grid='1', &
1418  longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy")
1419  call register_restart_field(meke%Au, vd, .false., restart_cs)
1420  endif
1421 
1422 end subroutine meke_alloc_register_restart
1423 
1424 !> Deallocates any variables allocated in MEKE_init or
1425 !! MEKE_alloc_register_restart.
1426 subroutine meke_end(MEKE, CS)
1427  type(meke_type), pointer :: meke !< A structure with MEKE-related fields.
1428  type(meke_cs), pointer :: cs !< The control structure for MOM_MEKE.
1429 
1430  if (associated(cs)) deallocate(cs)
1431 
1432  if (.not.associated(meke)) return
1433 
1434  if (associated(meke%MEKE)) deallocate(meke%MEKE)
1435  if (associated(meke%GM_src)) deallocate(meke%GM_src)
1436  if (associated(meke%mom_src)) deallocate(meke%mom_src)
1437  if (associated(meke%GME_snk)) deallocate(meke%GME_snk)
1438  if (associated(meke%Kh)) deallocate(meke%Kh)
1439  if (associated(meke%Kh_diff)) deallocate(meke%Kh_diff)
1440  if (associated(meke%Ku)) deallocate(meke%Ku)
1441  if (associated(meke%Au)) deallocate(meke%Au)
1442  deallocate(meke)
1443 
1444 end subroutine meke_end
1445 
1446 !> \namespace mom_meke
1447 !!
1448 !! \section section_MEKE The Mesoscale Eddy Kinetic Energy (MEKE) framework
1449 !!
1450 !! The MEKE framework accounts for the mean potential energy removed by
1451 !! the first order closures used to parameterize mesoscale eddies.
1452 !! It requires closure at the second order, namely dissipation and transport
1453 !! of eddy energy.
1454 !!
1455 !! Monitoring the sub-grid scale eddy energy budget provides a means to predict
1456 !! a sub-grid eddy-velocity scale which can be used in the lower order closures.
1457 !!
1458 !! \subsection section_MEKE_equations MEKE equations
1459 !!
1460 !! The eddy kinetic energy equation is:
1461 !! \f[ \partial_{\tilde{t}} E =
1462 !! \overbrace{ \dot{E}_b + \gamma_\eta \dot{E}_\eta + \gamma_v \dot{E}_v
1463 !! }^\text{sources}
1464 !! - \overbrace{ ( \lambda + C_d | U_d | \gamma_b^2 ) E
1465 !! }^\text{local dissipation}
1466 !! + \overbrace{ \nabla \cdot ( ( \kappa_E + \gamma_M \kappa_M ) \nabla E
1467 !! - \kappa_4 \nabla^3 E )
1468 !! }^\text{smoothing}
1469 !! \f]
1470 !! where \f$ E \f$ is the eddy kinetic energy (variable <code>MEKE</code>) with units of
1471 !! m<sup>2</sup>s<sup>-2</sup>,
1472 !! and \f$\tilde{t} = a t\f$ is a scaled time. The non-dimensional factor
1473 !! \f$ a\geq 1 \f$ is used to accelerate towards equilibrium.
1474 !!
1475 !! The MEKE equation is two-dimensional and obtained by depth averaging the
1476 !! the three-dimensional eddy energy equation. In the following expressions
1477 !! \f$ \left< \phi \right> = \frac{1}{H} \int^\eta_{-D} \phi \, dz \f$ maps
1478 !! three dimensional terms into the two-dimensional quantities needed.
1479 !!
1480 !! \subsubsection section_MEKE_source_terms MEKE source terms
1481 !!
1482 !! The source term \f$ \dot{E}_b \f$ is a constant background source
1483 !! of energy intended to avoid the limit \f$E\rightarrow 0\f$.
1484 !!
1485 !! The "GM" source term
1486 !! \f[ \dot{E}_\eta = - \left< \overline{w^\prime b^\prime} \right>
1487 !! = \left< \kappa_h N^2S^2 \right>
1488 !! \approx \left< \kappa_h g\prime |\nabla_\sigma \eta|^2 \right>\f]
1489 !! equals the mean potential energy removed by the Gent-McWilliams closure,
1490 !! and is excluded/included in the MEKE budget by the efficiency parameter
1491 !! \f$ \gamma_\eta \in [0,1] \f$.
1492 !!
1493 !! The "frictional" source term
1494 !! \f[ \dot{E}_{v} = \left< \partial_i u_j \tau_{ij} \right> \f]
1495 !! equals the mean kinetic energy removed by lateral viscous fluxes, and
1496 !! is excluded/included in the MEKE budget by the efficiency parameter
1497 !! \f$ \gamma_v \in [0,1] \f$.
1498 !!
1499 !! \subsubsection section_MEKE_dissipation_terms MEKE dissipation terms
1500 !!
1501 !! The local dissipation of \f$ E \f$ is parameterized through a linear
1502 !! damping, \f$\lambda\f$, and bottom drag, \f$ C_d | U_d | \gamma_b^2 \f$.
1503 !! The \f$ \gamma_b \f$ accounts for the weak projection of the column-mean
1504 !! eddy velocty to the bottom. In other words, the bottom velocity is
1505 !! estimated as \f$ \gamma_b U_e \f$.
1506 !! The bottom drag coefficient, \f$ C_d \f$ is the same as that used in the bottom
1507 !! friction in the mean model equations.
1508 !!
1509 !! The bottom drag velocity scale, \f$ U_d \f$, has contributions from the
1510 !! resolved state and \f$ E \f$:
1511 !! \f[ U_d = \sqrt{ U_b^2 + |u|^2_{z=-D} + |\gamma_b U_e|^2 } .\f]
1512 !! where the eddy velocity scale, \f$ U_e \f$, is given by:
1513 !! \f[ U_e = \sqrt{ 2 E } .\f]
1514 !! \f$ U_b \f$ is a constant background bottom velocity scale and is
1515 !! typically not used (i.e. set to zero).
1516 !!
1517 !! Following Jansen et al., 2015, the projection of eddy energy on to the bottom
1518 !! is given by the ratio of bottom energy to column mean energy:
1519 !! \f[
1520 !! \gamma_b^2 = \frac{E_b}{E} = \gamma_{d0}
1521 !! + \left( 1 + c_{b} \frac{L_d}{L_f} \right)^{-\frac{4}{5}}
1522 !! ,
1523 !! \f]
1524 !! \f[
1525 !! \gamma_b^2 \leftarrow \max{\left( \gamma_b^2, \gamma_{min}^2 \right)}
1526 !! .
1527 !! \f]
1528 !!
1529 !! \subsection section_MEKE_smoothing MEKE smoothing terms
1530 !!
1531 !! \f$ E \f$ is laterally diffused by a diffusivity \f$ \kappa_E + \gamma_M
1532 !! \kappa_M \f$ where \f$ \kappa_E \f$ is a constant diffusivity and the term
1533 !! \f$ \gamma_M \kappa_M \f$ is a "self diffusion" using the diffusivity
1534 !! calculated in the section \ref section_MEKE_diffusivity.
1535 !! \f$ \kappa_4 \f$ is a constant bi-harmonic diffusivity.
1536 !!
1537 !! \subsection section_MEKE_diffusivity Diffusivity derived from MEKE
1538 !!
1539 !! The predicted eddy velocity scale, \f$ U_e \f$, can be combined with a
1540 !! mixing length scale to form a diffusivity.
1541 !! The primary use of a MEKE derived diffusivity is for use in thickness
1542 !! diffusion (module mom_thickness_diffuse) and optionally in along
1543 !! isopycnal mixing of tracers (module mom_tracer_hor_diff).
1544 !! The original form used (enabled with MEKE_OLD_LSCALE=True):
1545 !!
1546 !! \f[ \kappa_M = \gamma_\kappa \sqrt{ \gamma_t^2 U_e^2 A_\Delta } \f]
1547 !!
1548 !! where \f$ A_\Delta \f$ is the area of the grid cell.
1549 !! Following Jansen et al., 2015, we now use
1550 !!
1551 !! \f[ \kappa_M = \gamma_\kappa l_M \sqrt{ \gamma_t^2 U_e^2 } \f]
1552 !!
1553 !! where \f$ \gamma_\kappa \in [0,1] \f$ is a non-dimensional factor and,
1554 !! following Jansen et al., 2015, \f$\gamma_t^2\f$ is the ratio of barotropic
1555 !! eddy energy to column mean eddy energy given by
1556 !! \f[
1557 !! \gamma_t^2 = \frac{E_t}{E} = \left( 1 + c_{t} \frac{L_d}{L_f} \right)^{-\frac{1}{4}}
1558 !! ,
1559 !! \f]
1560 !! \f[
1561 !! \gamma_t^2 \leftarrow \max{\left( \gamma_t^2, \gamma_{min}^2 \right)}
1562 !! .
1563 !! \f]
1564 !!
1565 !! The length-scale is a configurable combination of multiple length scales:
1566 !!
1567 !! \f[
1568 !! l_M = \left(
1569 !! \frac{\alpha_d}{L_d}
1570 !! + \frac{\alpha_f}{L_f}
1571 !! + \frac{\alpha_R}{L_R}
1572 !! + \frac{\alpha_e}{L_e}
1573 !! + \frac{\alpha_\Delta}{L_\Delta}
1574 !! + \frac{\delta[L_c]}{L_c}
1575 !! \right)^{-1}
1576 !! \f]
1577 !!
1578 !! where
1579 !!
1580 !! \f{eqnarray*}{
1581 !! L_d & = & \sqrt{\frac{c_g^2}{f^2+2\beta c_g}} \sim \frac{ c_g }{f} \\\\
1582 !! L_R & = & \sqrt{\frac{U_e}{\beta^*}} \\\\
1583 !! L_e & = & \frac{U_e}{|S| N} \\\\
1584 !! L_f & = & \frac{H}{c_d} \\\\
1585 !! L_\Delta & = & \sqrt{A_\Delta} .
1586 !! \f}
1587 !!
1588 !! \f$L_c\f$ is a constant and \f$\delta[L_c]\f$ is the impulse function so that the term
1589 !! \f$\frac{\delta[L_c]}{L_c}\f$ evaluates to \f$\frac{1}{L_c}\f$ when \f$L_c\f$ is non-zero
1590 !! but is dropped if \f$L_c=0\f$.
1591 !!
1592 !! \f$\beta^*\f$ is the effective \f$\beta\f$ that combines both the planetary vorticity
1593 !! gradient (i.e. \f$\beta=\nabla f\f$) and the topographic \f$\beta\f$ effect,
1594 !! with the latter weighed by a weighting constant, \f$c_\beta\f$, that varies
1595 !! from 0 to 1, so that \f$c_\beta=0\f$ means the topographic \f$\beta\f$ effect is ignored,
1596 !! while \f$c_\beta=1\f$ means it is fully considered. The new \f$\beta^*\f$ therefore
1597 !! takes the form of
1598 !!
1599 !! \f[
1600 !! \beta^* = \sqrt{( \partial_xf - c_\beta\frac{f}{D}\partial_xD )^2 +
1601 !! ( \partial_yf - c_\beta\frac{f}{D}\partial_yD )^2}
1602 !! \f]
1603 !! where \f$D\f$ is water column depth at T points.
1604 !!
1605 !! \subsection section_MEKE_viscosity Viscosity derived from MEKE
1606 !!
1607 !! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be
1608 !! used to form a harmonic eddy viscosity,
1609 !!
1610 !! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } \f]
1611 !!
1612 !! as well as a biharmonic eddy viscosity,
1613 !!
1614 !! \f[ \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 } \f]
1615 !!
1616 !! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance
1617 !!
1618 !! Note that in steady-state (or when \f$ a>>1 \f$) and there is no
1619 !! diffusion of \f$ E \f$ then
1620 !! \f[ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
1621 !! \gamma_v \dot{E}_v }{ \lambda + C_d|U_d|\gamma_b^2 } . \f]
1622 !!
1623 !! In the linear drag limit, where
1624 !! \f$ U_e << \min(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$, the equilibrium becomes
1625 !! \f$ \overline{E} \approx \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
1626 !! \gamma_v \dot{E}_v }{ \lambda + C_d \sqrt{ U_b^2 + |u|^2_{z=-D} } } \f$.
1627 !!
1628 !! In the nonlinear drag limit, where \f$ U_e >> \max(U_b, |u|_{z=-D}, C_d^{-1}\lambda) \f$,
1629 !! the equilibrium becomes
1630 !! \f$ \overline{E} \approx \left( \frac{ \dot{E}_b + \gamma_\eta \dot{E}_\eta +
1631 !! \gamma_v \dot{E}_v }{ \sqrt{2} C_d \gamma_b^3 } \right)^\frac{2}{3} \f$.
1632 !!
1633 !! \subsubsection section_MEKE_module_parameters MEKE module parameters
1634 !!
1635 !! | Symbol | Module parameter |
1636 !! | ------ | --------------- |
1637 !! | - | <code>USE_MEKE</code> |
1638 !! | \f$ a \f$ | <code>MEKE_DTSCALE</code> |
1639 !! | \f$ \dot{E}_b \f$ | <code>MEKE_BGSRC</code> |
1640 !! | \f$ \gamma_\eta \f$ | <code>MEKE_GMCOEFF</code> |
1641 !! | \f$ \gamma_v \f$ | <code>MEKE_FrCOEFF</code> |
1642 !! | \f$ \lambda \f$ | <code>MEKE_DAMPING</code> |
1643 !! | \f$ U_b \f$ | <code>MEKE_USCALE</code> |
1644 !! | \f$ \gamma_{d0} \f$ | <code>MEKE_CD_SCALE</code> |
1645 !! | \f$ c_{b} \f$ | <code>MEKE_CB</code> |
1646 !! | \f$ c_{t} \f$ | <code>MEKE_CT</code> |
1647 !! | \f$ \kappa_E \f$ | <code>MEKE_KH</code> |
1648 !! | \f$ \kappa_4 \f$ | <code>MEKE_K4</code> |
1649 !! | \f$ \gamma_\kappa \f$ | <code>MEKE_KHCOEFF</code> |
1650 !! | \f$ \gamma_M \f$ | <code>MEKE_KHMEKE_FAC</code> |
1651 !! | \f$ \gamma_u \f$ | <code>MEKE_VISCOSITY_COEFF_KU</code> |
1652 !! | \f$ \gamma_4 \f$ | <code>MEKE_VISCOSITY_COEFF_AU</code> |
1653 !! | \f$ \gamma_{min}^2 \f$| <code>MEKE_MIN_GAMMA2</code> |
1654 !! | \f$ \alpha_d \f$ | <code>MEKE_ALPHA_DEFORM</code> |
1655 !! | \f$ \alpha_f \f$ | <code>MEKE_ALPHA_FRICT</code> |
1656 !! | \f$ \alpha_R \f$ | <code>MEKE_ALPHA_RHINES</code> |
1657 !! | \f$ \alpha_e \f$ | <code>MEKE_ALPHA_EADY</code> |
1658 !! | \f$ \alpha_\Delta \f$ | <code>MEKE_ALPHA_GRID</code> |
1659 !! | \f$ L_c \f$ | <code>MEKE_FIXED_MIXING_LENGTH</code> |
1660 !! | \f$ c_\beta \f$ | <code>MEKE_TOPOGRAPHIC_BETA</code> |
1661 !! | - | <code>MEKE_KHTH_FAC</code> |
1662 !! | - | <code>MEKE_KHTR_FAC</code> |
1663 !!
1664 !! | Symbol | Model parameter |
1665 !! | ------ | --------------- |
1666 !! | \f$ C_d \f$ | <code>CDRAG</code> |
1667 !!
1668 !! \subsection section_MEKE_references References
1669 !!
1670 !! Jansen, M. F., A. J. Adcroft, R. Hallberg, and I. M. Held, 2015: Parameterization of eddy fluxes based on a
1671 !! mesoscale energy budget. Ocean Modelling, 92, 28--41, http://doi.org/10.1016/j.ocemod.2015.05.007 .
1672 !!
1673 !! Marshall, D. P., and A. J. Adcroft, 2010: Parameterization of ocean eddies: Potential vorticity mixing, energetics
1674 !! and Arnold first stability theorem. Ocean Modelling, 32, 188--204, http://doi.org/10.1016/j.ocemod.2010.02.001 .
1675 
1676 end module mom_meke
1677 
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:75
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:208
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_meke_types::meke_type
This type is used to exchange information related to the MEKE calculations.
Definition: MOM_MEKE_types.F90:8
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_meke::meke_cs
Control structure that contains MEKE parameters and diagnostics handles.
Definition: MOM_MEKE.F90:31
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_meke
Implements the Mesoscale Eddy Kinetic Energy framework with topographic beta effect included in compu...
Definition: MOM_MEKE.F90:4
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:16
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:111
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_io::vardesc
Type for describing a variable, typically a tracer.
Definition: MOM_io.F90:53
mom_restart::query_initialized
Indicate whether a field has been read from a restart file.
Definition: MOM_restart.F90:127
mom_domains::create_group_pass
Set up a group of halo updates.
Definition: MOM_domains.F90:84
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240
mom_file_parser::read_param
An overloaded interface to read various types of parameters.
Definition: MOM_file_parser.F90:90