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