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