MOM6
MOM_lateral_mixing_coeffs.F90
1 !> Variable mixing coefficients
2 module mom_lateral_mixing_coeffs
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_debugging, only : hchksum, uvchksum
7 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
8 use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, post_data
9 use mom_diag_mediator, only : diag_ctrl, time_type, query_averaging_enabled
10 use mom_domains, only : create_group_pass, do_group_pass
11 use mom_domains, only : group_pass_type, pass_var, pass_vector
12 use mom_file_parser, only : get_param, log_version, param_file_type
13 use mom_interface_heights, only : find_eta
14 use mom_isopycnal_slopes, only : calc_isoneutral_slopes
15 use mom_grid, only : ocean_grid_type
16 use mom_unit_scaling, only : unit_scale_type
17 use mom_variables, only : thermo_var_ptrs
18 use mom_verticalgrid, only : verticalgrid_type
19 use mom_wave_speed, only : wave_speed, wave_speed_cs, wave_speed_init
20 
21 implicit none ; private
22 
23 #include <MOM_memory.h>
24 
25 !> Variable mixing coefficients
26 type, public :: varmix_cs
27  logical :: use_variable_mixing !< If true, use the variable mixing.
28  logical :: resoln_scaled_kh !< If true, scale away the Laplacian viscosity
29  !! when the deformation radius is well resolved.
30  logical :: resoln_scaled_khth !< If true, scale away the thickness diffusivity
31  !! when the deformation radius is well resolved.
32  logical :: resoln_scaled_khtr !< If true, scale away the tracer diffusivity
33  !! when the deformation radius is well resolved.
34  logical :: interpolate_res_fn !< If true, interpolate the resolution function
35  !! to the velocity points from the thickness
36  !! points; otherwise interpolate the wave
37  !! speed and calculate the resolution function
38  !! independently at each point.
39  logical :: use_stored_slopes !< If true, stores isopycnal slopes in this structure.
40  logical :: resoln_use_ebt !< If true, uses the equivalent barotropic wave speed instead
41  !! of first baroclinic wave for calculating the resolution fn.
42  logical :: khth_use_ebt_struct !< If true, uses the equivalent barotropic structure
43  !! as the vertical structure of thickness diffusivity.
44  logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first
45  !! baroclinic wave speed and populate CS%cg1.
46  !! This parameter is set depending on other parameters.
47  logical :: calculate_rd_dx !< If true, calculates Rd/dx and populate CS%Rd_dx_h.
48  !! This parameter is set depending on other parameters.
49  logical :: calculate_res_fns !< If true, calculate all the resolution factors.
50  !! This parameter is set depending on other parameters.
51  logical :: calculate_eady_growth_rate !< If true, calculate all the Eady growth rate.
52  !! This parameter is set depending on other parameters.
53  real, dimension(:,:), pointer :: &
54  sn_u => null(), & !< S*N at u-points [s-1]
55  sn_v => null(), & !< S*N at v-points [s-1]
56  l2u => null(), & !< Length scale^2 at u-points [m2]
57  l2v => null(), & !< Length scale^2 at v-points [m2]
58  cg1 => null(), & !< The first baroclinic gravity wave speed [m s-1].
59  res_fn_h => null(), & !< Non-dimensional function of the ratio the first baroclinic
60  !! deformation radius to the grid spacing at h points [nondim].
61  res_fn_q => null(), & !< Non-dimensional function of the ratio the first baroclinic
62  !! deformation radius to the grid spacing at q points [nondim].
63  res_fn_u => null(), & !< Non-dimensional function of the ratio the first baroclinic
64  !! deformation radius to the grid spacing at u points [nondim].
65  res_fn_v => null(), & !< Non-dimensional function of the ratio the first baroclinic
66  !! deformation radius to the grid spacing at v points [nondim].
67  beta_dx2_h => null(), & !< The magnitude of the gradient of the Coriolis parameter
68  !! times the grid spacing squared at h points [m T-1 ~> m s-1].
69  beta_dx2_q => null(), & !< The magnitude of the gradient of the Coriolis parameter
70  !! times the grid spacing squared at q points [m T-1 ~> m s-1].
71  beta_dx2_u => null(), & !< The magnitude of the gradient of the Coriolis parameter
72  !! times the grid spacing squared at u points [m T-1 ~> m s-1].
73  beta_dx2_v => null(), & !< The magnitude of the gradient of the Coriolis parameter
74  !! times the grid spacing squared at v points [m T-1 ~> m s-1].
75  f2_dx2_h => null(), & !< The Coriolis parameter squared times the grid
76  !! spacing squared at h [m2 T-2 ~> m2 s-2].
77  f2_dx2_q => null(), & !< The Coriolis parameter squared times the grid
78  !! spacing squared at q [m2 T-2 ~> m2 s-2].
79  f2_dx2_u => null(), & !< The Coriolis parameter squared times the grid
80  !! spacing squared at u [m2 T-2 ~> m2 s-2].
81  f2_dx2_v => null(), & !< The Coriolis parameter squared times the grid
82  !! spacing squared at v [m2 T-2 ~> m2 s-2].
83  rd_dx_h => null() !< Deformation radius over grid spacing [nondim]
84 
85  real, dimension(:,:,:), pointer :: &
86  slope_x => null(), & !< Zonal isopycnal slope [nondim]
87  slope_y => null(), & !< Meridional isopycnal slope [nondim]
88  n2_u => null(), & !< Brunt-Vaisala frequency at u-points [s-2]
89  n2_v => null(), & !< Brunt-Vaisala frequency at v-points [s-2]
90  ebt_struct => null() !< Vertical structure function to scale diffusivities with [nondim]
91  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
92  laplac3_const_u !< Laplacian metric-dependent constants [nondim]
93 
94  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
95  laplac3_const_v !< Laplacian metric-dependent constants [nondim]
96 
97  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
98  kh_u_qg !< QG Leith GM coefficient at u-points [m2 s-1]
99 
100  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
101  kh_v_qg !< QG Leith GM coefficient at v-points [m2 s-1]
102 
103  ! Parameters
104  logical :: use_visbeck !< Use Visbeck formulation for thickness diffusivity
105  integer :: varmix_ktop !< Top layer to start downward integrals
106  real :: visbeck_l_scale !< Fixed length scale in Visbeck formula
107  real :: res_coef_khth !< A non-dimensional number that determines the function
108  !! of resolution, used for thickness and tracer mixing, as:
109  !! F = 1 / (1 + (Res_coef_khth*Ld/dx)^Res_fn_power)
110  real :: res_coef_visc !< A non-dimensional number that determines the function
111  !! of resolution, used for lateral viscosity, as:
112  !! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power)
113  real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [m2 s-1]
114  integer :: res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any
115  !! positive integer power may be used, but even powers
116  !! and especially 2 are coded to be more efficient.
117  integer :: res_fn_power_visc !< The power of dx/Ld in the Kh resolution function. Any
118  !! positive integer power may be used, but even powers
119  !! and especially 2 are coded to be more efficient.
120  real :: visbeck_s_max !< Upper bound on slope used in Eady growth rate [nondim].
121 
122  ! Leith parameters
123  logical :: use_qg_leith_gm !< If true, uses the QG Leith viscosity as the GM coefficient
124  logical :: use_beta_in_qg_leith !< If true, includes the beta term in the QG Leith GM coefficient
125 
126  ! Diagnostics
127  !>@{
128  !! Diagnostic identifier
129  integer :: id_sn_u=-1, id_sn_v=-1, id_l2u=-1, id_l2v=-1, id_res_fn = -1
130  integer :: id_n2_u=-1, id_n2_v=-1, id_s2_u=-1, id_s2_v=-1
131  integer :: id_rd_dx=-1, id_kh_u_qg = -1, id_kh_v_qg = -1
132  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
133  !! timing of diagnostic output.
134  !>@}
135 
136  type(wave_speed_cs), pointer :: wave_speed_csp => null() !< Wave speed control structure
137  type(group_pass_type) :: pass_cg1 !< For group halo pass
138  logical :: debug !< If true, write out checksums of data for debugging
139 end type varmix_cs
140 
141 public varmix_init, calc_slope_functions, calc_resoln_function
142 public calc_qg_leith_viscosity
143 
144 contains
145 
146 !> Calculates and stores the non-dimensional resolution functions
147 subroutine calc_resoln_function(h, tv, G, GV, US, CS)
148  type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
149  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
150  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
151  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
152  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
153  type(varmix_cs), pointer :: cs !< Variable mixing coefficients
154 
155  ! Local variables
156  ! Depending on the power-function being used, dimensional rescaling may be limited, so some
157  ! of the following variables have units that depend on that power.
158  real :: cg1_q ! The gravity wave speed interpolated to q points [m T-1 ~> m s-1] or [m s-1].
159  real :: cg1_u ! The gravity wave speed interpolated to u points [m T-1 ~> m s-1] or [m s-1].
160  real :: cg1_v ! The gravity wave speed interpolated to v points [m T-1 ~> m s-1] or [m s-1].
161  real :: dx_term ! A term in the denominator [m2 T-2 ~> m2 s-2] or [m2 s-2]
162  integer :: power_2
163  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
164  integer :: i, j, k
165  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
166  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
167 
168  if (.not. associated(cs)) call mom_error(fatal, "calc_resoln_function:"// &
169  "Module must be initialized before it is used.")
170  if (cs%calculate_cg1) then
171  if (.not. associated(cs%cg1)) call mom_error(fatal, &
172  "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
173  if (cs%khth_use_ebt_struct) then
174  if (.not. associated(cs%ebt_struct)) call mom_error(fatal, &
175  "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
176  if (cs%Resoln_use_ebt) then
177  ! Both resolution fn and vertical structure are using EBT
178  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct)
179  else
180  ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode
181  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct, &
182  use_ebt_mode=.true.)
183  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
184  endif
185  call pass_var(cs%ebt_struct, g%Domain)
186  else
187  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
188  endif
189 
190  call create_group_pass(cs%pass_cg1, cs%cg1, g%Domain)
191  call do_group_pass(cs%pass_cg1, g%Domain)
192  endif
193 
194  ! Calculate and store the ratio between deformation radius and grid-spacing
195  ! at h-points [nondim].
196  if (cs%calculate_rd_dx) then
197  if (.not. associated(cs%Rd_dx_h)) call mom_error(fatal, &
198  "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
199 !$OMP parallel default(none) shared(is,ie,js,je,CS)
200 !$OMP do
201  do j=js-1,je+1 ; do i=is-1,ie+1
202  cs%Rd_dx_h(i,j) = us%T_to_s*cs%cg1(i,j) / &
203  (sqrt(cs%f2_dx2_h(i,j) + us%T_to_s*cs%cg1(i,j)*cs%beta_dx2_h(i,j)))
204  enddo ; enddo
205 !$OMP end parallel
206  if (query_averaging_enabled(cs%diag)) then
207  if (cs%id_Rd_dx > 0) call post_data(cs%id_Rd_dx, cs%Rd_dx_h, cs%diag)
208  endif
209  endif
210 
211  if (.not. cs%calculate_res_fns) return
212 
213  if (.not. associated(cs%Res_fn_h)) call mom_error(fatal, &
214  "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
215  if (.not. associated(cs%Res_fn_q)) call mom_error(fatal, &
216  "calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
217  if (.not. associated(cs%Res_fn_u)) call mom_error(fatal, &
218  "calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
219  if (.not. associated(cs%Res_fn_v)) call mom_error(fatal, &
220  "calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
221  if (.not. associated(cs%f2_dx2_h)) call mom_error(fatal, &
222  "calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
223  if (.not. associated(cs%f2_dx2_q)) call mom_error(fatal, &
224  "calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
225  if (.not. associated(cs%f2_dx2_u)) call mom_error(fatal, &
226  "calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
227  if (.not. associated(cs%f2_dx2_v)) call mom_error(fatal, &
228  "calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
229  if (.not. associated(cs%beta_dx2_h)) call mom_error(fatal, &
230  "calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
231  if (.not. associated(cs%beta_dx2_q)) call mom_error(fatal, &
232  "calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
233  if (.not. associated(cs%beta_dx2_u)) call mom_error(fatal, &
234  "calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
235  if (.not. associated(cs%beta_dx2_v)) call mom_error(fatal, &
236  "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")
237 
238  ! Do this calculation on the extent used in MOM_hor_visc.F90, and
239  ! MOM_tracer.F90 so that no halo update is needed.
240 
241 !$OMP parallel default(none) shared(is,ie,js,je,Ieq,Jeq,CS) &
242 !$OMP private(dx_term,cg1_q,power_2,cg1_u,cg1_v)
243  if (cs%Res_fn_power_visc >= 100) then
244 !$OMP do
245  do j=js-1,je+1 ; do i=is-1,ie+1
246  dx_term = cs%f2_dx2_h(i,j) + us%T_to_s*cs%cg1(i,j)*cs%beta_dx2_h(i,j)
247  if ((cs%Res_coef_visc * us%T_to_s*cs%cg1(i,j))**2 > dx_term) then
248  cs%Res_fn_h(i,j) = 0.0
249  else
250  cs%Res_fn_h(i,j) = 1.0
251  endif
252  enddo ; enddo
253 !$OMP do
254  do j=js-1,jeq ; do i=is-1,ieq
255  cg1_q = us%T_to_s * 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + &
256  (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
257  dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
258  if ((cs%Res_coef_visc * cg1_q)**2 > dx_term) then
259  cs%Res_fn_q(i,j) = 0.0
260  else
261  cs%Res_fn_q(i,j) = 1.0
262  endif
263  enddo ; enddo
264  elseif (cs%Res_fn_power_visc == 2) then
265 !$OMP do
266  do j=js-1,je+1 ; do i=is-1,ie+1
267  dx_term = cs%f2_dx2_h(i,j) + us%T_to_s*cs%cg1(i,j)*cs%beta_dx2_h(i,j)
268  cs%Res_fn_h(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * us%T_to_s*cs%cg1(i,j))**2)
269  enddo ; enddo
270 !$OMP do
271  do j=js-1,jeq ; do i=is-1,ieq
272  cg1_q = us%T_to_s * 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + &
273  (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
274  dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
275  cs%Res_fn_q(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cg1_q)**2)
276  enddo ; enddo
277  elseif (mod(cs%Res_fn_power_visc, 2) == 0) then
278  power_2 = cs%Res_fn_power_visc / 2
279 !$OMP do
280  do j=js-1,je+1 ; do i=is-1,ie+1
281  dx_term = (us%s_to_T**2*cs%f2_dx2_h(i,j) + cs%cg1(i,j)*us%s_to_T*cs%beta_dx2_h(i,j))**power_2
282  cs%Res_fn_h(i,j) = dx_term / &
283  (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**cs%Res_fn_power_visc)
284  enddo ; enddo
285 !$OMP do
286  do j=js-1,jeq ; do i=is-1,ieq
287  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + &
288  (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
289  dx_term = (us%s_to_T**2*cs%f2_dx2_q(i,j) + cg1_q * us%s_to_T*cs%beta_dx2_q(i,j))**power_2
290  cs%Res_fn_q(i,j) = dx_term / &
291  (dx_term + (cs%Res_coef_visc * cg1_q)**cs%Res_fn_power_visc)
292  enddo ; enddo
293  else
294 !$OMP do
295  do j=js-1,je+1 ; do i=is-1,ie+1
296  dx_term = (us%s_to_T*sqrt(cs%f2_dx2_h(i,j) + &
297  us%T_to_s*cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**cs%Res_fn_power_visc
298  cs%Res_fn_h(i,j) = dx_term / &
299  (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**cs%Res_fn_power_visc)
300  enddo ; enddo
301 !$OMP do
302  do j=js-1,jeq ; do i=is-1,ieq
303  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + &
304  (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
305  dx_term = (us%s_to_T*sqrt(cs%f2_dx2_q(i,j) + &
306  us%T_to_s*cg1_q * cs%beta_dx2_q(i,j)))**cs%Res_fn_power_visc
307  cs%Res_fn_q(i,j) = dx_term / &
308  (dx_term + (cs%Res_coef_visc * cg1_q)**cs%Res_fn_power_visc)
309  enddo ; enddo
310  endif
311 
312  if (cs%interpolate_Res_fn) then
313  do j=js,je ; do i=is-1,ieq
314  cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
315  enddo ; enddo
316  do j=js-1,jeq ; do i=is,ie
317  cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
318  enddo ; enddo
319  else ! .not.CS%interpolate_Res_fn
320  if (cs%Res_fn_power_khth >= 100) then
321 !$OMP do
322  do j=js,je ; do i=is-1,ieq
323  cg1_u = 0.5 * us%T_to_s * (cs%cg1(i,j) + cs%cg1(i+1,j))
324  dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
325  if ((cs%Res_coef_khth * cg1_u)**2 > dx_term) then
326  cs%Res_fn_u(i,j) = 0.0
327  else
328  cs%Res_fn_u(i,j) = 1.0
329  endif
330  enddo ; enddo
331 !$OMP do
332  do j=js-1,jeq ; do i=is,ie
333  cg1_v = 0.5 * us%T_to_s * (cs%cg1(i,j) + cs%cg1(i,j+1))
334  dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
335  if ((cs%Res_coef_khth * cg1_v)**2 > dx_term) then
336  cs%Res_fn_v(i,j) = 0.0
337  else
338  cs%Res_fn_v(i,j) = 1.0
339  endif
340  enddo ; enddo
341  elseif (cs%Res_fn_power_khth == 2) then
342 !$OMP do
343  do j=js,je ; do i=is-1,ieq
344  cg1_u = 0.5 * us%T_to_s * (cs%cg1(i,j) + cs%cg1(i+1,j))
345  dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
346  cs%Res_fn_u(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_u)**2)
347  enddo ; enddo
348 !$OMP do
349  do j=js-1,jeq ; do i=is,ie
350  cg1_v = 0.5 * us%T_to_s * (cs%cg1(i,j) + cs%cg1(i,j+1))
351  dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
352  cs%Res_fn_v(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_v)**2)
353  enddo ; enddo
354  elseif (mod(cs%Res_fn_power_khth, 2) == 0) then
355  power_2 = cs%Res_fn_power_khth / 2
356 !$OMP do
357  do j=js,je ; do i=is-1,ieq
358  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
359  dx_term = (us%s_to_T**2*cs%f2_dx2_u(i,j) + cg1_u * us%s_to_T*cs%beta_dx2_u(i,j))**power_2
360  cs%Res_fn_u(i,j) = dx_term / &
361  (dx_term + (cs%Res_coef_khth * cg1_u)**cs%Res_fn_power_khth)
362  enddo ; enddo
363 !$OMP do
364  do j=js-1,jeq ; do i=is,ie
365  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
366  dx_term = (us%s_to_T**2*cs%f2_dx2_v(i,j) + cg1_v * us%s_to_T*cs%beta_dx2_v(i,j))**power_2
367  cs%Res_fn_v(i,j) = dx_term / &
368  (dx_term + (cs%Res_coef_khth * cg1_v)**cs%Res_fn_power_khth)
369  enddo ; enddo
370  else
371 !$OMP do
372  do j=js,je ; do i=is-1,ieq
373  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
374  dx_term = (us%s_to_T*sqrt(cs%f2_dx2_u(i,j) + &
375  us%T_to_s*cg1_u * cs%beta_dx2_u(i,j)))**cs%Res_fn_power_khth
376  cs%Res_fn_u(i,j) = dx_term / &
377  (dx_term + (cs%Res_coef_khth * cg1_u)**cs%Res_fn_power_khth)
378  enddo ; enddo
379 !$OMP do
380  do j=js-1,jeq ; do i=is,ie
381  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
382  dx_term = (us%s_to_T*sqrt(cs%f2_dx2_v(i,j) + &
383  us%T_to_s*cg1_v * cs%beta_dx2_v(i,j)))**cs%Res_fn_power_khth
384  cs%Res_fn_v(i,j) = dx_term / &
385  (dx_term + (cs%Res_coef_khth * cg1_v)**cs%Res_fn_power_khth)
386  enddo ; enddo
387  endif
388  endif
389 !$OMP end parallel
390 
391  if (query_averaging_enabled(cs%diag)) then
392  if (cs%id_Res_fn > 0) call post_data(cs%id_Res_fn, cs%Res_fn_h, cs%diag)
393  endif
394 
395 end subroutine calc_resoln_function
396 
397 !> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
398 !! style scaling of diffusivity
399 subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS)
400  type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
401  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
402  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
403  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
404  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
405  real, intent(in) :: dt !< Time increment [s]
406  type(varmix_cs), pointer :: cs !< Variable mixing coefficients
407  ! Local variables
408  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
409  e ! The interface heights relative to mean sea level [Z ~> m].
410  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: n2_u ! Square of Brunt-Vaisala freq at u-points [s-2]
411  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: n2_v ! Square of Brunt-Vaisala freq at v-points [s-2]
412 
413  if (.not. associated(cs)) call mom_error(fatal, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//&
414  "Module must be initialized before it is used.")
415 
416  if (cs%calculate_Eady_growth_rate) then
417  call find_eta(h, tv, g, gv, us, e, halo_size=2)
418  if (cs%use_stored_slopes) then
419  call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, &
420  cs%slope_x, cs%slope_y, n2_u, n2_v, 1)
421  call calc_visbeck_coeffs(h, cs%slope_x, cs%slope_y, n2_u, n2_v, g, gv, cs)
422 ! call calc_slope_functions_using_just_e(h, G, CS, e, .false.)
423  else
424  !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y)
425  call calc_slope_functions_using_just_e(h, g, gv, us, cs, e, .true.)
426  endif
427  endif
428 
429  if (query_averaging_enabled(cs%diag)) then
430  if (cs%id_SN_u > 0) call post_data(cs%id_SN_u, cs%SN_u, cs%diag)
431  if (cs%id_SN_v > 0) call post_data(cs%id_SN_v, cs%SN_v, cs%diag)
432  if (cs%id_L2u > 0) call post_data(cs%id_L2u, cs%L2u, cs%diag)
433  if (cs%id_L2v > 0) call post_data(cs%id_L2v, cs%L2v, cs%diag)
434  if (cs%id_N2_u > 0) call post_data(cs%id_N2_u, cs%N2_u, cs%diag)
435  if (cs%id_N2_v > 0) call post_data(cs%id_N2_v, cs%N2_v, cs%diag)
436  endif
437 
438 end subroutine calc_slope_functions
439 
440 !> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al.
441 subroutine calc_visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, CS)
442  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
443  type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
444  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
445  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: slope_x !< Zonal isoneutral slope
446  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: N2_u !< Brunt-Vaisala frequency at u-points [s-2]
447  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: slope_y !< Meridional isoneutral slope
448  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: N2_v !< Brunt-Vaisala frequency at v-points [s-2]
449  type(varmix_cs), pointer :: CS !< Variable mixing coefficients
450 
451  ! Local variables
452  real :: S2 ! Interface slope squared [nondim]
453  real :: N2 ! Brunt-Vaisala frequency [s-1]
454  real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
455  real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
456  integer :: is, ie, js, je, nz
457  integer :: i, j, k, kb_max
458  real :: S2max, wNE, wSE, wSW, wNW
459  real :: H_u(SZIB_(G)), H_v(SZI_(G))
460  real :: S2_u(SZIB_(G), SZJ_(G))
461  real :: S2_v(SZI_(G), SZJB_(G))
462 
463  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
464  "Module must be initialized before it is used.")
465  if (.not. cs%calculate_Eady_growth_rate) return
466  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
467  "%SN_u is not associated with use_variable_mixing.")
468  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
469  "%SN_v is not associated with use_variable_mixing.")
470 
471  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
472 
473  s2max = cs%Visbeck_S_max**2
474 
475  !$OMP parallel do default(shared)
476  do j=js-1,je+1 ; do i=is-1,ie+1
477  cs%SN_u(i,j) = 0.0
478  cs%SN_v(i,j) = 0.0
479  enddo ; enddo
480 
481  ! To set the length scale based on the deformation radius, use wave_speed to
482  ! calculate the first-mode gravity wave speed and then blend the equatorial
483  ! and midlatitude deformation radii, using calc_resoln_function as a template.
484 
485  !$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
486  do j = js,je
487  do i=is-1,ie
488  cs%SN_u(i,j) = 0. ; h_u(i) = 0. ; s2_u(i,j) = 0.
489  enddo
490  do k=2,nz ; do i=is-1,ie
491  hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
492  hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
493  h_geom = sqrt( hdn * hup )
494  !H_geom = H_geom * sqrt(N2) ! WKB-ish
495  !H_geom = H_geom * N2 ! WKB-ish
496  wse = g%mask2dCv(i+1,j-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) )
497  wnw = g%mask2dCv(i ,j ) * ( (h(i ,j,k)*h(i ,j+1,k)) * (h(i ,j,k-1)*h(i ,j+1,k-1)) )
498  wne = g%mask2dCv(i+1,j ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) )
499  wsw = g%mask2dCv(i ,j-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) )
500  s2 = slope_x(i,j,k)**2 + &
501  ((wnw*slope_y(i,j,k)**2 + wse*slope_y(i+1,j-1,k)**2) + &
502  (wne*slope_y(i+1,j,k)**2 + wsw*slope_y(i,j-1,k)**2) ) / &
503  ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
504  if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
505 
506  n2 = max(0., n2_u(i,j,k))
507  cs%SN_u(i,j) = cs%SN_u(i,j) + sqrt( s2*n2 )*h_geom
508  s2_u(i,j) = s2_u(i,j) + s2*h_geom
509  h_u(i) = h_u(i) + h_geom
510  enddo ; enddo
511  do i=is-1,ie
512  if (h_u(i)>0.) then
513  cs%SN_u(i,j) = g%mask2dCu(i,j) * cs%SN_u(i,j) / h_u(i)
514  s2_u(i,j) = g%mask2dCu(i,j) * s2_u(i,j) / h_u(i)
515  else
516  cs%SN_u(i,j) = 0.
517  endif
518  enddo
519  enddo
520 
521  !$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
522  do j = js-1,je
523  do i=is,ie
524  cs%SN_v(i,j) = 0. ; h_v(i) = 0. ; s2_v(i,j) = 0.
525  enddo
526  do k=2,nz ; do i=is,ie
527  hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
528  hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
529  h_geom = sqrt( hdn * hup )
530  !H_geom = H_geom * sqrt(N2) ! WKB-ish
531  !H_geom = H_geom * N2 ! WKB-ish
532  wse = g%mask2dCu(i,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) )
533  wnw = g%mask2dCu(i-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) )
534  wne = g%mask2dCu(i,j+1) * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) )
535  wsw = g%mask2dCu(i-1,j) * ( (h(i,j ,k)*h(i-1,j ,k)) * (h(i,j ,k-1)*h(i-1,j ,k-1)) )
536  s2 = slope_y(i,j,k)**2 + &
537  ((wse*slope_x(i,j,k)**2 + wnw*slope_x(i-1,j+1,k)**2) + &
538  (wne*slope_x(i,j+1,k)**2 + wsw*slope_x(i-1,j,k)**2) ) / &
539  ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
540  if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
541 
542  n2 = max(0., n2_v(i,j,k))
543  cs%SN_v(i,j) = cs%SN_v(i,j) + sqrt( s2*n2 )*h_geom
544  s2_v(i,j) = s2_v(i,j) + s2*h_geom
545  h_v(i) = h_v(i) + h_geom
546  enddo ; enddo
547  do i=is,ie
548  if (h_v(i)>0.) then
549  cs%SN_v(i,j) = g%mask2dCv(i,j) * cs%SN_v(i,j) / h_v(i)
550  s2_v(i,j) = g%mask2dCv(i,j) * s2_v(i,j) / h_v(i)
551  else
552  cs%SN_v(i,j) = 0.
553  endif
554  enddo
555  enddo
556 
557 ! Offer diagnostic fields for averaging.
558  if (query_averaging_enabled(cs%diag)) then
559  if (cs%id_S2_u > 0) call post_data(cs%id_S2_u, s2_u, cs%diag)
560  if (cs%id_S2_v > 0) call post_data(cs%id_S2_v, s2_v, cs%diag)
561  endif
562 
563  if (cs%debug) then
564  call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, g%HI, haloshift=1)
565  call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", n2_u, n2_v, g%HI)
566  call uvchksum("calc_Visbeck_coeffs SN_[uv]", cs%SN_u, cs%SN_v, g%HI)
567  endif
568 
569 end subroutine calc_visbeck_coeffs
570 
571 !> The original calc_slope_function() that calculated slopes using
572 !! interface positions only, not accounting for density variations.
573 subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slopes)
574  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
575  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
576  type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
577  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
578  type(varmix_cs), pointer :: CS !< Variable mixing coefficients
579  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface position [Z ~> m]
580  logical, intent(in) :: calculate_slopes !< If true, calculate slopes internally
581  !! otherwise use slopes stored in CS
582  ! Local variables
583  real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics)
584  real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics)
585  real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
586  real :: h_neglect ! A thickness that is so small it is usually lost
587  ! in roundoff and can be neglected [H ~> m or kg m-2].
588  real :: S2 ! Interface slope squared [nondim]
589  real :: N2 ! Brunt-Vaisala frequency squared [T-2 ~> s-2]
590  real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
591  real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
592  real :: Z_to_L ! A conversion factor between from units for e to the
593  ! units for lateral distances.
594  real :: one_meter ! One meter in thickness units [H ~> m or kg m-2].
595  integer :: is, ie, js, je, nz
596  integer :: i, j, k, kb_max
597  real :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(G))
598  real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(G))
599 
600  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
601  "Module must be initialized before it is used.")
602  if (.not. cs%calculate_Eady_growth_rate) return
603  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
604  "%SN_u is not associated with use_variable_mixing.")
605  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
606  "%SN_v is not associated with use_variable_mixing.")
607 
608  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
609 
610  one_meter = 1.0 * gv%m_to_H
611  h_neglect = gv%H_subroundoff
612  h_cutoff = real(2*nz) * (gv%Angstrom_H + h_neglect)
613  z_to_l = us%Z_to_m
614 
615  ! To set the length scale based on the deformation radius, use wave_speed to
616  ! calculate the first-mode gravity wave speed and then blend the equatorial
617  ! and midlatitude deformation radii, using calc_resoln_function as a template.
618 
619  !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
620  do k=nz,cs%VarMix_Ktop,-1
621 
622  if (calculate_slopes) then
623  ! Calculate the interface slopes E_x and E_y and u- and v- points respectively
624  do j=js-1,je+1 ; do i=is-1,ie
625  e_x(i,j) = z_to_l*(e(i+1,j,k)-e(i,j,k))*g%IdxCu(i,j)
626  ! Mask slopes where interface intersects topography
627  if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
628  enddo ; enddo
629  do j=js-1,je ; do i=is-1,ie+1
630  e_y(i,j) = z_to_l*(e(i,j+1,k)-e(i,j,k))*g%IdyCv(i,j)
631  ! Mask slopes where interface intersects topography
632  if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
633  enddo ; enddo
634  else
635  do j=js-1,je+1 ; do i=is-1,ie
636  e_x(i,j) = cs%slope_x(i,j,k)
637  if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
638  enddo ; enddo
639  do j=js-1,je ; do i=is-1,ie+1
640  e_y(i,j) = cs%slope_y(i,j,k)
641  if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
642  enddo ; enddo
643  endif
644 
645  ! Calculate N*S*h from this layer and add to the sum
646  do j=js,je ; do i=is-1,ie
647  s2 = ( e_x(i,j)**2 + 0.25*( &
648  (e_y(i,j)**2+e_y(i+1,j-1)**2)+(e_y(i+1,j)**2+e_y(i,j-1)**2) ) )
649  hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
650  hup = 2.*h(i+1,j,k)*h(i+1,j,k-1) / (h(i+1,j,k) + h(i+1,j,k-1) + h_neglect)
651  h_geom = sqrt(hdn*hup)
652  n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
653  if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < h_cutoff) &
654  s2 = 0.0
655  s2n2_u_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
656  enddo ; enddo
657  do j=js-1,je ; do i=is,ie
658  s2 = ( e_y(i,j)**2 + 0.25*( &
659  (e_x(i,j)**2+e_x(i-1,j+1)**2)+(e_x(i,j+1)**2+e_x(i-1,j)**2) ) )
660  hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
661  hup = 2.*h(i,j+1,k)*h(i,j+1,k-1) / (h(i,j+1,k) + h(i,j+1,k-1) + h_neglect)
662  h_geom = sqrt(hdn*hup)
663  n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
664  if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < h_cutoff) &
665  s2 = 0.0
666  s2n2_v_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
667  enddo ; enddo
668 
669  enddo ! k
670  !$OMP parallel do default(shared)
671  do j=js,je
672  do i=is-1,ie ; cs%SN_u(i,j) = 0.0 ; enddo
673  do k=nz,cs%VarMix_Ktop,-1 ; do i=is-1,ie
674  cs%SN_u(i,j) = cs%SN_u(i,j) + s2n2_u_local(i,j,k)
675  enddo ; enddo
676  ! SN above contains S^2*N^2*H, convert to vertical average of S*N
677  do i=is-1,ie
678  !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom_Z ) ))
679  !The code below behaves better than the line above. Not sure why? AJA
680  if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
681  cs%SN_u(i,j) = g%mask2dCu(i,j) * us%s_to_T * sqrt( cs%SN_u(i,j) / &
682  (max(g%bathyT(i,j), g%bathyT(i+1,j))) )
683  else
684  cs%SN_u(i,j) = 0.0
685  endif
686  enddo
687  enddo
688  !$OMP parallel do default(shared)
689  do j=js-1,je
690  do i=is,ie ; cs%SN_v(i,j) = 0.0 ; enddo
691  do k=nz,cs%VarMix_Ktop,-1 ; do i=is,ie
692  cs%SN_v(i,j) = cs%SN_v(i,j) + s2n2_v_local(i,j,k)
693  enddo ; enddo
694  do i=is,ie
695  !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom_Z ) ))
696  !The code below behaves better than the line above. Not sure why? AJA
697  if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
698  cs%SN_v(i,j) = g%mask2dCv(i,j) * us%s_to_T * sqrt( cs%SN_v(i,j) / &
699  (max(g%bathyT(i,j), g%bathyT(i,j+1))) )
700  else
701  cs%SN_v(i,j) = 0.0
702  endif
703  enddo
704  enddo
705 
706 end subroutine calc_slope_functions_using_just_e
707 
708 !> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
709 subroutine calc_qg_leith_viscosity(CS, G, GV, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy)
710  type(varmix_cs), pointer :: cs !< Variable mixing coefficients
711  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
712  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
713 ! real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow [m s-1]
714 ! real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow [m s-1]
715  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
716  integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
717  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence
718  !! (d/dx(du/dx + dv/dy)) [m-1 s-1]
719  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence
720  !! (d/dy(du/dx + dv/dy)) [m-1 s-1]
721  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
722  !! (d/dx(dv/dx - du/dy)) [m-1 s-1]
723  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
724  !! (d/dy(dv/dx - du/dy)) [m-1 s-1]
725 ! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity
726  !! at h-points [m2 s-1]
727 ! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity
728  !! at q-points [m2 s-1]
729 ! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity
730  !! at h-points [m4 s-1]
731 ! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity
732  !! at q-points [m4 s-1]
733 
734  ! Local variables
735 ! real, dimension(SZIB_(G),SZJB_(G)) :: vort_xy, & ! Vertical vorticity (dv/dx - du/dy) [s-1]
736 ! dudy, & ! Meridional shear of zonal velocity [s-1]
737 ! dvdx ! Zonal shear of meridional velocity [s-1]
738  real, dimension(SZI_(G),SZJB_(G)) :: &
739 ! vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [m-1 s-1]
740 ! div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [m-1 s-1]
741  dslopey_dz, & ! z-derivative of y-slope at v-points [m-1]
742  h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
743  beta_v, & ! Beta at v-points [m-1 s-1]
744  grad_vort_mag_v, & ! mag. of vort. grad. at v-points [s-1]
745  grad_div_mag_v ! mag. of div. grad. at v-points [s-1]
746 
747  real, dimension(SZIB_(G),SZJ_(G)) :: &
748 ! vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [m-1 s-1]
749 ! div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [m-1 s-1]
750  dslopex_dz, & ! z-derivative of x-slope at u-points (m-1)
751  h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
752  beta_u, & ! Beta at u-points [m-1 s-1]
753  grad_vort_mag_u, & ! mag. of vort. grad. at u-points [s-1]
754  grad_div_mag_u ! mag. of div. grad. at u-points [s-1]
755 ! real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points [s-1]
756 ! real :: mod_Leith, DY_dxBu, DX_dyBu, vert_vort_mag
757  real :: h_at_slope_above, h_at_slope_below, ih, f
758  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq,nz
759  real :: inv_pi3
760 
761  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
762  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
763  nz = g%ke
764 
765  inv_pi3 = 1.0/((4.0*atan(1.0))**3)
766 
767  !### I believe this halo update to be unnecessary. -RWH
768  call pass_var(h, g%Domain)
769 
770  if ((k > 1) .and. (k < nz)) then
771 
772  ! Add in stretching term for the QG Leith vsicosity
773 ! if (CS%use_QG_Leith) then
774 
775  !### do j=js-1,je+1 ; do I=is-2,Ieq+1
776  do j=js-2,jeq+2 ; do i=is-2,ieq+1
777  h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / &
778  ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
779  + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + gv%H_subroundoff )
780  h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / &
781  ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
782  + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + gv%H_subroundoff )
783  ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_m )
784  dslopex_dz(i,j) = 2. * ( cs%slope_x(i,j,k) - cs%slope_x(i,j,k+1) ) * ih
785  h_at_u(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
786  enddo ; enddo
787 
788  !### do J=js-2,Jeq+1 ; do i=is-1,ie+1
789  do j=js-2,jeq+1 ; do i=is-2,ieq+2
790  h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / &
791  ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
792  + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + gv%H_subroundoff )
793  h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / &
794  ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
795  + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + gv%H_subroundoff )
796  ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_m )
797  dslopey_dz(i,j) = 2. * ( cs%slope_y(i,j,k) - cs%slope_y(i,j,k+1) ) * ih
798  h_at_v(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
799  enddo ; enddo
800 
801  !### do J=js-1,je ; do i=is-1,Ieq+1
802  do j=js-2,jeq+1 ; do i=is-1,ieq+1
803  f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
804  vort_xy_dx(i,j) = vort_xy_dx(i,j) - f * &
805  ( ( h_at_u(i,j) * dslopex_dz(i,j) + h_at_u(i-1,j+1) * dslopex_dz(i-1,j+1) ) &
806  + ( h_at_u(i-1,j) * dslopex_dz(i-1,j) + h_at_u(i,j+1) * dslopex_dz(i,j+1) ) ) / &
807  ( ( h_at_u(i,j) + h_at_u(i-1,j+1) ) + ( h_at_u(i-1,j) + h_at_u(i,j+1) ) + gv%H_subroundoff)
808  enddo ; enddo
809 
810  !### do j=js-1,Jeq+1 ; do I=is-1,ie
811  do j=js-1,jeq+1 ; do i=is-2,ieq+1
812  f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
813  !### I think that this should be vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * &
814  vort_xy_dy(i,j) = vort_xy_dx(i,j) - f * &
815  ( ( h_at_v(i,j) * dslopey_dz(i,j) + h_at_v(i+1,j-1) * dslopey_dz(i+1,j-1) ) &
816  + ( h_at_v(i,j-1) * dslopey_dz(i,j-1) + h_at_v(i+1,j) * dslopey_dz(i+1,j) ) ) / &
817  ( ( h_at_v(i,j) + h_at_v(i+1,j-1) ) + ( h_at_v(i,j-1) + h_at_v(i+1,j) ) + gv%H_subroundoff)
818  enddo ; enddo
819  endif ! k > 1
820 
821  !### I believe this halo update to be unnecessary. -RWH
822  call pass_vector(vort_xy_dy,vort_xy_dx,g%Domain)
823 
824  if (cs%use_QG_Leith_GM) then
825 
826  do j=js,je ; do i=is-1,ieq
827  grad_vort_mag_u(i,j) = sqrt(vort_xy_dy(i,j)**2 + (0.25*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j) &
828  + vort_xy_dx(i,j-1) + vort_xy_dx(i+1,j-1)))**2)
829  grad_div_mag_u(i,j) = sqrt(div_xx_dx(i,j)**2 + (0.25*(div_xx_dy(i,j) + div_xx_dy(i+1,j) &
830  + div_xx_dy(i,j-1) + div_xx_dy(i+1,j-1)))**2)
831  if (cs%use_beta_in_QG_Leith) then
832  beta_u(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i+1,j))**2) + &
833  (0.5*(g%dF_dy(i,j)+g%dF_dy(i+1,j))**2) )
834  cs%KH_u_QG(i,j,k) = min(grad_vort_mag_u(i,j) + grad_div_mag_u(i,j), beta_u(i,j)*3) &
835  * cs%Laplac3_const_u(i,j) * inv_pi3
836  else
837  cs%KH_u_QG(i,j,k) = (grad_vort_mag_u(i,j) + grad_div_mag_u(i,j)) &
838  * cs%Laplac3_const_u(i,j) * inv_pi3
839  endif
840  enddo ; enddo
841 
842  do j=js-1,jeq ; do i=is,ie
843  grad_vort_mag_v(i,j) = sqrt(vort_xy_dx(i,j)**2 + (0.25*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j) &
844  + vort_xy_dy(i,j+1) + vort_xy_dy(i-1,j+1)))**2)
845  grad_div_mag_v(i,j) = sqrt(div_xx_dy(i,j)**2 + (0.25*(div_xx_dx(i,j) + div_xx_dx(i-1,j) &
846  + div_xx_dx(i,j+1) + div_xx_dx(i-1,j+1)))**2)
847  if (cs%use_beta_in_QG_Leith) then
848  beta_v(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i,j+1))**2) + &
849  (0.5*(g%dF_dy(i,j)+g%dF_dy(i,j+1))**2) )
850  cs%KH_v_QG(i,j,k) = min(grad_vort_mag_v(i,j) + grad_div_mag_v(i,j), beta_v(i,j)*3) &
851  * cs%Laplac3_const_v(i,j) * inv_pi3
852  else
853  cs%KH_v_QG(i,j,k) = (grad_vort_mag_v(i,j) + grad_div_mag_v(i,j)) &
854  * cs%Laplac3_const_v(i,j) * inv_pi3
855  endif
856  enddo ; enddo
857  ! post diagnostics
858 
859  if (k==nz) then
860  if (cs%id_KH_v_QG > 0) call post_data(cs%id_KH_v_QG, cs%KH_v_QG, cs%diag)
861  if (cs%id_KH_u_QG > 0) call post_data(cs%id_KH_u_QG, cs%KH_u_QG, cs%diag)
862  endif
863  endif
864 
865 end subroutine calc_qg_leith_viscosity
866 
867 !> Initializes the variables mixing coefficients container
868 subroutine varmix_init(Time, G, GV, US, param_file, diag, CS)
869  type(time_type), intent(in) :: time !< Current model time
870  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
871  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
872  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
873  type(param_file_type), intent(in) :: param_file !< Parameter file handles
874  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
875  type(varmix_cs), pointer :: cs !< Variable mixing coefficients
876  ! Local variables
877  real :: khtr_slope_cff, khth_slope_cff, oneortwo, n2_filter_depth
878  real :: khtr_passivity_coeff
879  real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
880  ! default value is roughly (pi / (the age of the universe)).
881  logical :: gill_equatorial_ld, use_fgnv_streamfn, use_meke, in_use
882  real :: mle_front_length
883  real :: leith_lap_const ! The non-dimensional coefficient in the Leith viscosity
884  real :: grid_sp_u2, grid_sp_u3
885  real :: grid_sp_v2, grid_sp_v3 ! Intermediate quantities for Leith metrics
886 ! This include declares and sets the variable "version".
887 #include "version_variable.h"
888  character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
889  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
890  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
891  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
892  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
893  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
894  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
895 
896  if (associated(cs)) then
897  call mom_error(warning, "VarMix_init called with an associated "// &
898  "control structure.")
899  return
900  endif
901 
902  allocate(cs)
903  in_use = .false. ! Set to true to avoid deallocating
904  cs%diag => diag ! Diagnostics pointer
905  cs%calculate_cg1 = .false.
906  cs%calculate_Rd_dx = .false.
907  cs%calculate_res_fns = .false.
908  cs%calculate_Eady_growth_rate = .false.
909 
910  ! Read all relevant parameters and write them to the model log.
911  call log_version(param_file, mdl, version, "")
912  call get_param(param_file, mdl, "USE_VARIABLE_MIXING", cs%use_variable_mixing,&
913  "If true, the variable mixing code will be called. This "//&
914  "allows diagnostics to be created even if the scheme is "//&
915  "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
916  "this is set to true regardless of what is in the "//&
917  "parameter file.", default=.false.)
918  call get_param(param_file, mdl, "USE_VISBECK", cs%use_Visbeck,&
919  "If true, use the Visbeck et al. (1997) formulation for \n"//&
920  "thickness diffusivity.", default=.false.)
921  call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
922  "If true, the Laplacian lateral viscosity is scaled away "//&
923  "when the first baroclinic deformation radius is well "//&
924  "resolved.", default=.false.)
925  call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
926  "If true, the interface depth diffusivity is scaled away "//&
927  "when the first baroclinic deformation radius is well "//&
928  "resolved.", default=.false.)
929  call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
930  "If true, the epipycnal tracer diffusivity is scaled "//&
931  "away when the first baroclinic deformation radius is "//&
932  "well resolved.", default=.false.)
933  call get_param(param_file, mdl, "RESOLN_USE_EBT", cs%Resoln_use_ebt, &
934  "If true, uses the equivalent barotropic wave speed instead "//&
935  "of first baroclinic wave for calculating the resolution fn.",&
936  default=.false.)
937  call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
938  "If true, uses the equivalent barotropic structure "//&
939  "as the vertical structure of thickness diffusivity.",&
940  default=.false.)
941  call get_param(param_file, mdl, "KHTH_SLOPE_CFF", khth_slope_cff, &
942  "The nondimensional coefficient in the Visbeck formula "//&
943  "for the interface depth diffusivity", units="nondim", &
944  default=0.0)
945  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", khtr_slope_cff, &
946  "The nondimensional coefficient in the Visbeck formula "//&
947  "for the epipycnal tracer diffusivity", units="nondim", &
948  default=0.0)
949  call get_param(param_file, mdl, "USE_STORED_SLOPES", cs%use_stored_slopes,&
950  "If true, the isopycnal slopes are calculated once and "//&
951  "stored for re-use. This uses more memory but avoids calling "//&
952  "the equation of state more times than should be necessary.", &
953  default=.false.)
954  call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
955  "A miniscule frequency that is used to avoid division by 0. The default "//&
956  "value is roughly (pi / (the age of the universe)).", &
957  default=1.0e-17, units="s-1", scale=us%T_to_s)
958  call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
959  default=.false., do_not_log=.true.)
960  cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn
961  call get_param(param_file, mdl, "USE_MEKE", use_meke, &
962  default=.false., do_not_log=.true.)
963  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
964  cs%calculate_Eady_growth_rate = cs%calculate_Eady_growth_rate .or. use_meke
965  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
966  default=0., do_not_log=.true.)
967  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
968  call get_param(param_file, mdl, "MLE_FRONT_LENGTH", mle_front_length, &
969  default=0., do_not_log=.true.)
970  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
971 
972  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
973 
974  if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct) then
975  in_use = .true.
976  call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
977  "The depth below which N2 is monotonized to avoid stratification "//&
978  "artifacts from altering the equivalent barotropic mode structure.",&
979  units="m", default=2000.)
980  allocate(cs%ebt_struct(isd:ied,jsd:jed,g%ke)) ; cs%ebt_struct(:,:,:) = 0.0
981  endif
982 
983  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
984  cs%calculate_Eady_growth_rate = .true.
985  call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
986  "If non-zero, is an upper bound on slopes used in the "//&
987  "Visbeck formula for diffusivity. This does not affect the "//&
988  "isopycnal slope calculation used within thickness diffusion.", &
989  units="nondim", default=0.0)
990  endif
991 
992  if (cs%use_stored_slopes) then
993  in_use = .true.
994  allocate(cs%slope_x(isdb:iedb,jsd:jed,g%ke+1)) ; cs%slope_x(:,:,:) = 0.0
995  allocate(cs%slope_y(isd:ied,jsdb:jedb,g%ke+1)) ; cs%slope_y(:,:,:) = 0.0
996  allocate(cs%N2_u(isdb:iedb,jsd:jed,g%ke+1)) ; cs%N2_u(:,:,:) = 0.0
997  allocate(cs%N2_v(isd:ied,jsdb:jedb,g%ke+1)) ; cs%N2_v(:,:,:) = 0.0
998  call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
999  "A diapycnal diffusivity that is used to interpolate "//&
1000  "more sensible values of T & S into thin layers.", &
1001  units="m2 s-1", default=1.0e-6, scale=us%m_to_Z**2)
1002  endif
1003 
1004  if (cs%calculate_Eady_growth_rate) then
1005  in_use = .true.
1006  allocate(cs%SN_u(isdb:iedb,jsd:jed)) ; cs%SN_u(:,:) = 0.0
1007  allocate(cs%SN_v(isd:ied,jsdb:jedb)) ; cs%SN_v(:,:) = 0.0
1008  cs%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, time, &
1009  'Inverse eddy time-scale, S*N, at u-points', 's-1')
1010  cs%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, time, &
1011  'Inverse eddy time-scale, S*N, at v-points', 's-1')
1012  call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
1013  "The layer number at which to start vertical integration "//&
1014  "of S*N for purposes of finding the Eady growth rate.", &
1015  units="nondim", default=2)
1016  endif
1017 
1018  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1019  in_use = .true.
1020  call get_param(param_file, mdl, "VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1021  "The fixed length scale in the Visbeck formula.", units="m", &
1022  default=0.0)
1023  allocate(cs%L2u(isdb:iedb,jsd:jed)) ; cs%L2u(:,:) = 0.0
1024  allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = 0.0
1025  if (cs%Visbeck_L_scale<0) then
1026  do j=js,je ; do i=is-1,ieq
1027  cs%L2u(i,j) = cs%Visbeck_L_scale**2*g%areaCu(i,j)
1028  enddo; enddo
1029  do j=js-1,jeq ; do i=is,ie
1030  cs%L2v(i,j) = cs%Visbeck_L_scale**2*g%areaCv(i,j)
1031  enddo; enddo
1032  else
1033  cs%L2u(:,:) = cs%Visbeck_L_scale**2
1034  cs%L2v(:,:) = cs%Visbeck_L_scale**2
1035  endif
1036 
1037  cs%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, time, &
1038  'Length scale squared for mixing coefficient, at u-points', 'm2')
1039  cs%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, time, &
1040  'Length scale squared for mixing coefficient, at v-points', 'm2')
1041  endif
1042 
1043  if (cs%use_stored_slopes) then
1044  cs%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, time, &
1045  'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', 's-2')
1046  cs%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, time, &
1047  'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', 's-2')
1048  cs%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, time, &
1049  'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 's-2')
1050  cs%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, time, &
1051  'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 's-2')
1052  endif
1053 
1054  oneortwo = 1.0
1055  if (cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. cs%Resoln_scaled_KhTr) then
1056  cs%calculate_Rd_dx = .true.
1057  cs%calculate_res_fns = .true.
1058  allocate(cs%Res_fn_h(isd:ied,jsd:jed)) ; cs%Res_fn_h(:,:) = 0.0
1059  allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb)) ; cs%Res_fn_q(:,:) = 0.0
1060  allocate(cs%Res_fn_u(isdb:iedb,jsd:jed)) ; cs%Res_fn_u(:,:) = 0.0
1061  allocate(cs%Res_fn_v(isd:ied,jsdb:jedb)) ; cs%Res_fn_v(:,:) = 0.0
1062  allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%beta_dx2_q(:,:) = 0.0
1063  allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed)) ; cs%beta_dx2_u(:,:) = 0.0
1064  allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb)) ; cs%beta_dx2_v(:,:) = 0.0
1065  allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%f2_dx2_q(:,:) = 0.0
1066  allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed)) ; cs%f2_dx2_u(:,:) = 0.0
1067  allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb)) ; cs%f2_dx2_v(:,:) = 0.0
1068 
1069  cs%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, time, &
1070  'Resolution function for scaling diffusivities', 'nondim')
1071 
1072  call get_param(param_file, mdl, "KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1073  "A coefficient that determines how KhTh is scaled away if "//&
1074  "RESOLN_SCALED_... is true, as "//&
1075  "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1076  units="nondim", default=1.0)
1077  call get_param(param_file, mdl, "KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1078  "The power of dx/Ld in the Kh resolution function. Any "//&
1079  "positive integer may be used, although even integers "//&
1080  "are more efficient to calculate. Setting this greater "//&
1081  "than 100 results in a step-function being used.", &
1082  units="nondim", default=2)
1083  call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1084  "A coefficient that determines how Kh is scaled away if "//&
1085  "RESOLN_SCALED_... is true, as "//&
1086  "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1087  "This function affects lateral viscosity, Kh, and not KhTh.", &
1088  units="nondim", default=cs%Res_coef_khth)
1089  call get_param(param_file, mdl, "VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1090  "The power of dx/Ld in the Kh resolution function. Any "//&
1091  "positive integer may be used, although even integers "//&
1092  "are more efficient to calculate. Setting this greater "//&
1093  "than 100 results in a step-function being used. "//&
1094  "This function affects lateral viscosity, Kh, and not KhTh.", &
1095  units="nondim", default=cs%Res_fn_power_khth)
1096  call get_param(param_file, mdl, "INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1097  "If true, interpolate the resolution function to the "//&
1098  "velocity points from the thickness points; otherwise "//&
1099  "interpolate the wave speed and calculate the resolution "//&
1100  "function independently at each point.", default=.true.)
1101  if (cs%interpolate_Res_fn) then
1102  if (cs%Res_coef_visc /= cs%Res_coef_khth) call mom_error(fatal, &
1103  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1104  "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1105  if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth) call mom_error(fatal, &
1106  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1107  "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1108  endif
1109  !### Change the default of GILL_EQUATORIAL_LD to True.
1110  call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1111  "If true, uses Gill's definition of the baroclinic "//&
1112  "equatorial deformation radius, otherwise, if false, use "//&
1113  "Pedlosky's definition. These definitions differ by a factor "//&
1114  "of 2 in front of the beta term in the denominator. Gill's "//&
1115  "is the more appropriate definition.", default=.false.)
1116  if (gill_equatorial_ld) then
1117  oneortwo = 2.0
1118  endif
1119 
1120  do j=js-1,jeq ; do i=is-1,ieq
1121  cs%f2_dx2_q(i,j) = (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * &
1122  max(g%CoriolisBu(i,j)**2, absurdly_small_freq**2)
1123  cs%beta_dx2_q(i,j) = oneortwo * (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * (sqrt(0.5 * &
1124  ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1125  ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1126  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1127  ((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) ) ))
1128  enddo ; enddo
1129 
1130  do j=js,je ; do i=is-1,ieq
1131  cs%f2_dx2_u(i,j) = (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * &
1132  max(0.5* (g%CoriolisBu(i,j)**2+g%CoriolisBu(i,j-1)**2), absurdly_small_freq**2)
1133  cs%beta_dx2_u(i,j) = oneortwo * (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * (sqrt( &
1134  0.25*( (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2 + &
1135  ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1136  (((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2 + &
1137  ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) ) + &
1138  ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 ))
1139  enddo ; enddo
1140 
1141  do j=js-1,jeq ; do i=is,ie
1142  cs%f2_dx2_v(i,j) = (g%dxCv(i,j)**2 + g%dyCv(i,j)**2) * &
1143  max(0.5*(g%CoriolisBu(i,j)**2+g%CoriolisBu(i-1,j)**2), absurdly_small_freq**2)
1144  cs%beta_dx2_v(i,j) = oneortwo * (g%dxCv(i,j)**2 + g%dyCv(i,j)**2) * (sqrt( &
1145  ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1146  0.25*( (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1147  ((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2) + &
1148  (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2 + &
1149  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1150  enddo ; enddo
1151 
1152  endif
1153 
1154  ! Resolution %Rd_dx_h
1155  cs%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, time, &
1156  'Ratio between deformation radius and grid spacing', 'm m-1')
1157  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
1158 
1159  if (cs%calculate_Rd_dx) then
1160  cs%calculate_cg1 = .true. ! We will need %cg1
1161  allocate(cs%Rd_dx_h(isd:ied,jsd:jed)) ; cs%Rd_dx_h(:,:) = 0.0
1162  allocate(cs%beta_dx2_h(isd:ied,jsd:jed)); cs%beta_dx2_h(:,:) = 0.0
1163  allocate(cs%f2_dx2_h(isd:ied,jsd:jed)) ; cs%f2_dx2_h(:,:) = 0.0
1164  do j=js-1,je+1 ; do i=is-1,ie+1
1165  cs%f2_dx2_h(i,j) = (g%dxT(i,j)**2 + g%dyT(i,j)**2) * &
1166  max(0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1167  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)), &
1168  absurdly_small_freq**2)
1169  cs%beta_dx2_h(i,j) = oneortwo * (g%dxT(i,j)**2 + g%dyT(i,j)**2) * (sqrt(0.5 * &
1170  ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1171  ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1172  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1173  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1174  enddo ; enddo
1175  endif
1176 
1177  if (cs%calculate_cg1) then
1178  in_use = .true.
1179  allocate(cs%cg1(isd:ied,jsd:jed)); cs%cg1(:,:) = 0.0
1180  call wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, mono_n2_depth=n2_filter_depth)
1181  endif
1182 
1183  ! Leith parameters
1184  call get_param(param_file, mdl, "USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
1185  "If true, use the QG Leith viscosity as the GM coefficient.", &
1186  default=.false.)
1187 
1188  if (cs%Use_QG_Leith_GM) then
1189  call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
1190  "The nondimensional Laplacian Leith constant, \n"//&
1191  "often set to 1.0", units="nondim", default=0.0)
1192 
1193  call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
1194  "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1195  default=.true.)
1196 
1197  alloc_(cs%Laplac3_const_u(isdb:iedb,jsd:jed)) ; cs%Laplac3_const_u(:,:) = 0.0
1198  alloc_(cs%Laplac3_const_v(isd:ied,jsdb:jedb)) ; cs%Laplac3_const_v(:,:) = 0.0
1199  alloc_(cs%KH_u_QG(isdb:iedb,jsd:jed,g%ke)) ; cs%KH_u_QG(:,:,:) = 0.0
1200  alloc_(cs%KH_v_QG(isd:ied,jsdb:jedb,g%ke)) ; cs%KH_v_QG(:,:,:) = 0.0
1201  ! register diagnostics
1202 
1203  cs%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, time, &
1204  'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1')
1205  cs%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, time, &
1206  'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1')
1207 
1208  do j=jsq,jeq+1 ; do i=is-1,ieq
1209  ! Static factors in the Leith schemes
1210  grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
1211  grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
1212  cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
1213  enddo ; enddo
1214  do j=js-1,jeq ; do i=isq,ieq+1
1215  ! Static factors in the Leith schemes
1216  grid_sp_v2 = g%dyCv(i,j)*g%dxCu(i,j)
1217  grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
1218  cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
1219  enddo ; enddo
1220 
1221  if (.not. cs%use_stored_slopes) call mom_error(fatal, &
1222  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1223  "USE_STORED_SLOPES must be True when using QG Leith.")
1224  endif
1225 
1226  ! If nothing is being stored in this class then deallocate
1227  if (in_use) then
1228  cs%use_variable_mixing = .true.
1229  else
1230  deallocate(cs)
1231  return
1232  endif
1233 
1234 end subroutine varmix_init
1235 
1236 !> \namespace mom_lateral_mixing_coeffs
1237 !!
1238 !! This module provides a container for various factors used in prescribing diffusivities, that are
1239 !! a function of the state (in particular the stratification and isoneutral slopes).
1240 !!
1241 !! \section section_Resolution_Function The resolution function
1242 !!
1243 !! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
1244 !! The square of the resolution parameter is
1245 !!
1246 !! \f[
1247 !! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
1248 !! \f]
1249 !!
1250 !! where the grid spacing is calculated as
1251 !!
1252 !! \f[
1253 !! \Delta^2 = \Delta x^2 + \Delta y^2 .
1254 !! \f]
1255 !!
1256 !! \todo Check this reference to Bob on/off paper.
1257 !! The resolution function used in scaling diffusivities (Hallberg, 2010) is
1258 !!
1259 !! \f[
1260 !! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
1261 !! \f]
1262 !!
1263 !! The resolution function can be applied independently to thickness diffusion (module mom_thickness_diffuse),
1264 !! tracer diffusion (mom_tracer_hordiff) lateral viscosity (mom_hor_visc).
1265 !!
1266 !! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
1267 !! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007
1268 !!
1269 !! | Symbol | Module parameter |
1270 !! | ------ | --------------- |
1271 !! | - | <code>USE_VARIABLE_MIXING</code> |
1272 !! | - | <code>RESOLN_SCALED_KH</code> |
1273 !! | - | <code>RESOLN_SCALED_KHTH</code> |
1274 !! | - | <code>RESOLN_SCALED_KHTR</code> |
1275 !! | \f$ \alpha \f$ | <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) |
1276 !! | \f$ p \f$ | <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) |
1277 !! | \f$ \alpha \f$ | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) |
1278 !! | \f$ p \f$ | <code>VISC_RES_FN_POWER</code> (for lateral viscosity) |
1279 !! | - | <code>GILL_EQUATORIAL_LD</code> |
1280 !!
1281 !!
1282 !!
1283 !! \section section_Vicbeck Visbeck diffusivity
1284 !!
1285 !! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
1286 !! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse() but calculated in this module.
1287 !!
1288 !! \f[
1289 !! \kappa_h = \alpha_s L_s^2 S N
1290 !! \f]
1291 !!
1292 !! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
1293 !!
1294 !! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
1295 !! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
1296 !!
1297 !! | Symbol | Module parameter |
1298 !! | ------ | --------------- |
1299 !! | - | <code>USE_VARIABLE_MIXING</code> |
1300 !! | \f$ \alpha_s \f$ | <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)|
1301 !! | \f$ \alpha_s \f$ | <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)|
1302 !! | \f$ L_{s} \f$ | <code>VISBECK_L_SCALE</code> |
1303 !! | \f$ S_{max} \f$ | <code>VISBECK_MAX_SLOPE</code> |
1304 !!
1305 !!
1306 !! \section section_vertical_structure_khth Vertical structure function for KhTh
1307 !!
1308 !! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
1309 !! velocity mode. The structure function is stored in the control structure for thie module (varmix_cs) but is
1310 !! calculated using subroutines in mom_wave_speed.
1311 !!
1312 !! | Symbol | Module parameter |
1313 !! | ------ | --------------- |
1314 !! | - | <code>KHTH_USE_EBT_STRUCT</code> |
1315 
1316 end module mom_lateral_mixing_coeffs
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_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:82
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
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_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_isopycnal_slopes
Calculations of isoneutral slopes and stratification.
Definition: MOM_isopycnal_slopes.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:54
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_lateral_mixing_coeffs::varmix_cs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:26
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_wave_speed::wave_speed_cs
Control structure for MOM_wave_speed.
Definition: MOM_wave_speed.F90:28
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_interface_heights::find_eta
Calculates the heights of sruface or all interfaces from layer thicknesses.
Definition: MOM_interface_heights.F90:21
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_wave_speed
Routines for calculating baroclinic wave speeds.
Definition: MOM_wave_speed.F90:2
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_lateral_mixing_coeffs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:2
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_interface_heights
Functions for calculating interface heights, including free surface height.
Definition: MOM_interface_heights.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