MOM_lateral_mixing_coeffs.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Variable mixing coefficients
7
8use mom_debugging, only : hchksum, uvchksum
9use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
10use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, post_data
12use mom_domains, only : create_group_pass, do_group_pass
13use mom_domains, only : group_pass_type, pass_var, pass_vector
14use mom_eos, only : calculate_density_derivs, eos_domain
15use mom_file_parser, only : get_param, log_version, param_file_type
16use mom_interface_heights, only : find_eta, thickness_to_dz
18use mom_grid, only : ocean_grid_type
23use mom_open_boundary, only : ocean_obc_type, obc_none
24use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
25use mom_meke_types, only : meke_type
26
27implicit none ; private
28
29#include <MOM_memory.h>
30
31!> Variable mixing coefficients
32type, public :: varmix_cs
33 logical :: initialized = .false. !< True if this control structure has been initialized.
34 logical :: use_variable_mixing !< If true, use the variable mixing.
35 logical :: resoln_scaling_used !< If true, a resolution function is used somewhere to scale
36 !! away one of the viscosities or diffusivities when the
37 !! deformation radius is well resolved.
38 logical :: resoln_scaled_kh !< If true, scale away the Laplacian viscosity
39 !! when the deformation radius is well resolved.
40 logical :: resoln_scaled_khth !< If true, scale away the thickness diffusivity
41 !! when the deformation radius is well resolved.
42 logical :: depth_scaled_khth !< If true, KHTH is scaled away when the depth is
43 !! shallower than a reference depth.
44 logical :: resoln_scaled_khtr !< If true, scale away the tracer diffusivity
45 !! when the deformation radius is well resolved.
46 logical :: interpolate_res_fn !< If true, interpolate the resolution function
47 !! to the velocity points from the thickness
48 !! points; otherwise interpolate the wave
49 !! speed and calculate the resolution function
50 !! independently at each point.
51 logical :: use_stored_slopes !< If true, stores isopycnal slopes in this structure.
52 logical :: resoln_use_ebt !< If true, use the equivalent barotropic wave speed instead of the
53 !! first baroclinic wave speed for calculating the resolution function.
54 logical :: khth_use_ebt_struct !< If true, uses the equivalent barotropic structure
55 !! as the vertical structure of thickness diffusivity.
56 logical :: kdgl90_use_ebt_struct !< If true, uses the equivalent barotropic structure
57 !! as the vertical structure of diffusivity in the GL90 scheme.
58 logical :: kdgl90_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
59 !! as the vertical structure of diffusivity in the GL90 scheme.
60 logical :: khth_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
61 !! as the vertical structure of thickness diffusivity.
62 logical :: khtr_use_ebt_struct !< If true, uses the equivalent barotropic structure
63 !! as the vertical structure of tracer diffusivity.
64 logical :: khtr_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
65 !! as the vertical structure of tracer diffusivity.
66 logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first
67 !! baroclinic wave speed and populate CS%cg1.
68 !! This parameter is set depending on other parameters.
69 logical :: calculate_rd_dx !< If true, calculates Rd/dx and populate CS%Rd_dx_h.
70 !! This parameter is set depending on other parameters.
71 logical :: calculate_res_fns !< If true, calculate all the resolution factors.
72 !! This parameter is set depending on other parameters.
73 logical :: calculate_depth_fns !< If true, calculate all the depth factors.
74 !! This parameter is set depending on other parameters.
75 logical :: calculate_eady_growth_rate !< If true, calculate all the Eady growth rates.
76 !! This parameter is set depending on other parameters.
77 logical :: use_stanley_iso !< If true, use Stanley parameterization in MOM_isopycnal_slopes
78 logical :: use_simpler_eady_growth_rate !< If true, use a simpler method to calculate the
79 !! Eady growth rate that avoids division by layer thickness.
80 !! This parameter is set depending on other parameters.
81 logical :: full_depth_eady_growth_rate !< If true, calculate the Eady growth rate based on an
82 !! average that includes contributions from sea-level changes
83 !! in its denominator, rather than just the nominal depth of
84 !! the bathymetry. This only applies when using the model
85 !! interface heights as a proxy for isopycnal slopes.
86 logical :: obc_friendly !< If true, use only interior data for thickness weighting and
87 !! to calculate stratification and other fields at open boundary
88 !! condition faces.
89 logical :: res_fn_obc_bug !< If false, use only interior data for calculating the resolution
90 !! functions at open boundary condition faces and vertices.
91 real :: cropping_distance !< Distance from surface or bottom to filter out outcropped or
92 !! incropped interfaces for the Eady growth rate calc [Z ~> m]
93 real :: h_min_n2 !< The minimum vertical distance to use in the denominator of the
94 !! buoyancy frequency used in the slope calculation [H ~> m or kg m-2]
95
96 real, allocatable :: sn_u(:,:) !< S*N at u-points [T-1 ~> s-1]
97 real, allocatable :: sn_v(:,:) !< S*N at v-points [T-1 ~> s-1]
98 real, allocatable :: l2u(:,:) !< Length scale^2 at u-points [L2 ~> m2]
99 real, allocatable :: l2v(:,:) !< Length scale^2 at v-points [L2 ~> m2]
100 real, allocatable :: cg1(:,:) !< The first baroclinic gravity wave speed [L T-1 ~> m s-1].
101 real, allocatable :: res_fn_h(:,:) !< Non-dimensional function of the ratio the first baroclinic
102 !! deformation radius to the grid spacing at h points [nondim].
103 real, allocatable :: res_fn_q(:,:) !< Non-dimensional function of the ratio the first baroclinic
104 !! deformation radius to the grid spacing at q points [nondim].
105 real, allocatable :: res_fn_u(:,:) !< Non-dimensional function of the ratio the first baroclinic
106 !! deformation radius to the grid spacing at u points [nondim].
107 real, allocatable :: res_fn_v(:,:) !< Non-dimensional function of the ratio the first baroclinic
108 !! deformation radius to the grid spacing at v points [nondim].
109 real, allocatable :: depth_fn_u(:,:) !< Non-dimensional function of the ratio of the depth to
110 !! a reference depth (maximum 1) at u points [nondim]
111 real, allocatable :: depth_fn_v(:,:) !< Non-dimensional function of the ratio of the depth to
112 !! a reference depth (maximum 1) at v points [nondim]
113 real, allocatable :: beta_dx2_h(:,:) !< The magnitude of the gradient of the Coriolis parameter
114 !! times the grid spacing squared at h points [L T-1 ~> m s-1].
115 real, allocatable :: beta_dx2_q(:,:) !< The magnitude of the gradient of the Coriolis parameter
116 !! times the grid spacing squared at q points [L T-1 ~> m s-1].
117 real, allocatable :: beta_dx2_u(:,:) !< The magnitude of the gradient of the Coriolis parameter
118 !! times the grid spacing squared at u points [L T-1 ~> m s-1].
119 real, allocatable :: beta_dx2_v(:,:) !< The magnitude of the gradient of the Coriolis parameter
120 !! times the grid spacing squared at v points [L T-1 ~> m s-1].
121 real, allocatable :: f2_dx2_h(:,:) !< The Coriolis parameter squared times the grid
122 !! spacing squared at h [L2 T-2 ~> m2 s-2].
123 real, allocatable :: f2_dx2_q(:,:) !< The Coriolis parameter squared times the grid
124 !! spacing squared at q [L2 T-2 ~> m2 s-2].
125 real, allocatable :: f2_dx2_u(:,:) !< The Coriolis parameter squared times the grid
126 !! spacing squared at u [L2 T-2 ~> m2 s-2].
127 real, allocatable :: f2_dx2_v(:,:) !< The Coriolis parameter squared times the grid
128 !! spacing squared at v [L2 T-2 ~> m2 s-2].
129 real, allocatable :: rd_dx_h(:,:) !< Deformation radius over grid spacing [nondim]
130
131 real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim]
132 real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim]
133 real, allocatable :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim]
134 real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure to scale diffusivities with [nondim]
135 real, allocatable :: bs_struct(:,:,:) !< Vertical structure function used in backscatter [nondim]
136 real, allocatable :: khth_struct(:,:,:) !< Vertical structure function used in thickness diffusivity [nondim]
137 real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim]
138 real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [nondim]
139 real :: bs_ebt_power !< Power to raise EBT vertical structure to. Default 0.0.
140 real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 1.0
141 logical :: interpolated_sqg_struct !< If true, interpolate properties to velocity points and then
142 !! interpolate the buoyancy frequencies and layer thicknesses
143 !! back to tracer points when calculating the SQG vertical
144 !! structure.
145 logical :: bs_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure.
146
147 real, allocatable :: laplac3_const_u(:,:) !< Laplacian metric-dependent constants at u-points [L3 ~> m3]
148 real, allocatable :: laplac3_const_v(:,:) !< Laplacian metric-dependent constants at u-points [L3 ~> m3]
149 real, allocatable :: kh_u_qg(:,:,:) !< QG Leith GM coefficient at u-points [L2 T-1 ~> m2 s-1]
150 real, allocatable :: kh_v_qg(:,:,:) !< QG Leith GM coefficient at v-points [L2 T-1 ~> m2 s-1]
151
152 ! Parameters
153 logical :: use_visbeck !< Use Visbeck formulation for thickness diffusivity
154 integer :: varmix_ktop !< Top layer to start downward integrals
155 real :: visbeck_l_scale !< Fixed length scale in Visbeck formula [L ~> m], or if negative a scaling
156 !! factor [nondim] relating this length scale squared to the cell area
157 real :: eady_gr_d_scale !< Depth over which to average SN [Z ~> m]
158 real :: res_coef_khth !< A coefficient [nondim] that determines the function
159 !! of resolution, used for thickness and tracer mixing, as:
160 !! F = 1 / (1 + (Res_coef_khth*Ld/dx)^Res_fn_power)
161 real :: res_coef_visc !< A coefficient [nondim] that determines the function
162 !! of resolution, used for lateral viscosity, as:
163 !! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power)
164 real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m]
165 real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim]
166 real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
167 integer :: res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any
168 !! positive integer power may be used, but even powers
169 !! and especially 2 are coded to be more efficient.
170 integer :: res_fn_power_visc !< The power of dx/Ld in the Kh resolution function. Any
171 !! positive integer power may be used, but even powers
172 !! and especially 2 are coded to be more efficient.
173 real :: visbeck_s_max !< Upper bound on slope used in Eady growth rate [Z L-1 ~> nondim].
174
175 ! Leith parameters
176 logical :: use_qg_leith_gm !< If true, uses the QG Leith viscosity as the GM coefficient
177 logical :: use_beta_in_qg_leith !< If true, includes the beta term in the QG Leith GM coefficient
178
179 ! Diagnostics
180 !>@{
181 !! Diagnostic identifier
182 integer :: id_sn_u=-1, id_sn_v=-1, id_l2u=-1, id_l2v=-1, id_res_fn = -1
183 integer :: id_n2_u=-1, id_n2_v=-1, id_s2_u=-1, id_s2_v=-1
184 integer :: id_dzu=-1, id_dzv=-1, id_dzsxn=-1, id_dzsyn=-1
185 integer :: id_rd_dx=-1, id_kh_u_qg = -1, id_kh_v_qg = -1
186 integer :: id_sqg_struct=-1, id_bs_struct=-1, id_khth_struct=-1, id_khtr_struct=-1
187 integer :: id_kdgl90_struct=-1
188 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
189 !! timing of diagnostic output.
190 !>@}
191
192 type(wave_speed_cs) :: wave_speed !< Wave speed control structure
193 type(group_pass_type) :: pass_cg1 !< For group halo pass
194 logical :: debug !< If true, write out checksums of data for debugging
195end type varmix_cs
196
199
200contains
201
202!> Calculates the non-dimensional depth functions.
203subroutine calc_depth_function(G, CS)
204 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
205 type(varmix_cs), intent(inout) :: cs !< Variable mixing control structure
206
207 ! Local variables
208 integer :: is, ie, js, je, isq, ieq, jsq, jeq
209 integer :: i, j
210 real :: h0 ! The depth above which KHTH is linearly scaled away [Z ~> m]
211 real :: h1, h2 ! Temporary total thicknesses [Z ~> m]
212 real :: expo ! exponent used in the depth dependent scaling [nondim]
213 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
214 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
215
216 if (.not. cs%initialized) call mom_error(fatal, "calc_depth_function: "// &
217 "Module must be initialized before it is used.")
218
219 if (.not. cs%calculate_depth_fns) return
220 if (.not. allocated(cs%Depth_fn_u)) call mom_error(fatal, &
221 "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
222 if (.not. allocated(cs%Depth_fn_v)) call mom_error(fatal, &
223 "calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.")
224
225 ! For efficiency, the reciprocal of H0 should be used instead.
226 h0 = cs%depth_scaled_khth_h0
227 expo = cs%depth_scaled_khth_exp
228!$OMP do
229 do j=js,je ; do i=is-1,ieq
230 h1 = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0)
231 h2 = max(g%meanSL(i+1,j) + g%bathyT(i+1,j), 0.0)
232 cs%Depth_fn_u(i,j) = (min(1.0, (0.5 * (h1 + h2)) / h0))**expo
233 enddo ; enddo
234!$OMP do
235 do j=js-1,jeq ; do i=is,ie
236 h1 = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0)
237 h2 = max(g%meanSL(i,j+1) + g%bathyT(i,j+1), 0.0)
238 cs%Depth_fn_v(i,j) = (min(1.0, (0.5 * (h1 + h2)) / h0))**expo
239 enddo ; enddo
240
241end subroutine calc_depth_function
242
243!> Calculates and stores the non-dimensional resolution functions
244subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, OBC, dt)
245 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
246 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
247 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
248 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
249 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
250 type(varmix_cs), intent(inout) :: cs !< Variable mixing control structure
251 type(meke_type), intent(in) :: meke !< MEKE struct
252 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure
253 real, intent(in) :: dt !< Time increment [T ~> s]
254
255 ! Local variables
256 ! Depending on the power-function being used, dimensional rescaling may be limited, so some
257 ! of the following variables have units that depend on that power.
258 real :: cg1_q(szib_(g),szjb_(g)) ! The gravity wave speed interpolated to q points [L T-1 ~> m s-1] or [m s-1].
259 real :: cg1_u(szib_(g),szj_(g)) ! The gravity wave speed interpolated to u points [L T-1 ~> m s-1] or [m s-1].
260 real :: cg1_v(szi_(g),szjb_(g)) ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
261 real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
262 logical :: apply_u_obc, apply_v_obc ! If true, OBCs will be used to set the wave speed at some points on this PE.
263 integer :: power_2
264 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
265 integer :: i, j, k
266 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
267 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
268
269 if (.not. cs%initialized) call mom_error(fatal, "calc_resoln_function: "// &
270 "Module must be initialized before it is used.")
271
272 if (cs%calculate_cg1) then
273 if (.not. allocated(cs%cg1)) call mom_error(fatal, &
274 "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
275 if (cs%khth_use_ebt_struct .or. cs%kdgl90_use_ebt_struct &
276 .or. cs%khtr_use_ebt_struct .or. cs%BS_EBT_power>0.) then
277 if (.not. allocated(cs%ebt_struct)) call mom_error(fatal, &
278 "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
279 if (cs%Resoln_use_ebt) then
280 ! Both resolution fn and vertical structure are using EBT
281 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed, modal_structure=cs%ebt_struct)
282 else
283 ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode
284 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed, modal_structure=cs%ebt_struct, &
285 use_ebt_mode=.true.)
286 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed)
287 endif
288 call pass_var(cs%ebt_struct, g%Domain)
289 else
290 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed)
291 endif
292
293 call create_group_pass(cs%pass_cg1, cs%cg1, g%Domain)
294 call do_group_pass(cs%pass_cg1, g%Domain)
295 endif
296 if (cs%BS_use_sqg_struct .or. cs%khth_use_sqg_struct .or. cs%khtr_use_sqg_struct &
297 .or. cs%kdgl90_use_sqg_struct .or. cs%id_sqg_struct>0) then
298 call calc_sqg_struct(h, tv, g, gv, us, cs, dt, meke, obc)
299 call pass_var(cs%sqg_struct, g%Domain)
300 endif
301
302 if (cs%BS_EBT_power>0.) then
303 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
304 cs%BS_struct(i,j,k) = cs%ebt_struct(i,j,k)**cs%BS_EBT_power
305 enddo ; enddo ; enddo
306 elseif (cs%BS_use_sqg_struct) then
307 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
308 cs%BS_struct(i,j,k) = cs%sqg_struct(i,j,k)
309 enddo ; enddo ; enddo
310 endif
311
312 if (cs%khth_use_ebt_struct) then
313 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
314 cs%khth_struct(i,j,k) = cs%ebt_struct(i,j,k)
315 enddo ; enddo ; enddo
316 elseif (cs%khth_use_sqg_struct) then
317 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
318 cs%khth_struct(i,j,k) = cs%sqg_struct(i,j,k)
319 enddo ; enddo ; enddo
320 endif
321
322 if (cs%khtr_use_ebt_struct) then
323 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
324 cs%khtr_struct(i,j,k) = cs%ebt_struct(i,j,k)
325 enddo ; enddo ; enddo
326 elseif (cs%khtr_use_sqg_struct) then
327 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
328 cs%khtr_struct(i,j,k) = cs%sqg_struct(i,j,k)
329 enddo ; enddo ; enddo
330 endif
331
332 if (cs%kdgl90_use_ebt_struct) then
333 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
334 cs%kdgl90_struct(i,j,k) = cs%ebt_struct(i,j,k)
335 enddo ; enddo ; enddo
336 elseif (cs%kdgl90_use_sqg_struct) then
337 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
338 cs%kdgl90_struct(i,j,k) = cs%sqg_struct(i,j,k)
339 enddo ; enddo ; enddo
340 endif
341
342 ! Calculate and store the ratio between deformation radius and grid-spacing
343 ! at h-points [nondim].
344 if (cs%calculate_rd_dx) then
345 if (.not. allocated(cs%Rd_dx_h)) call mom_error(fatal, &
346 "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
347 !$OMP parallel do default(shared)
348 do j=js-1,je+1 ; do i=is-1,ie+1
349 cs%Rd_dx_h(i,j) = cs%cg1(i,j) / &
350 (sqrt(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))
351 enddo ; enddo
352 if (query_averaging_enabled(cs%diag)) then
353 if (cs%id_Rd_dx > 0) call post_data(cs%id_Rd_dx, cs%Rd_dx_h, cs%diag)
354 endif
355 endif
356
357 if (.not. cs%calculate_res_fns) return
358
359 if (.not. allocated(cs%Res_fn_h)) call mom_error(fatal, &
360 "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
361 if (.not. allocated(cs%Res_fn_q)) call mom_error(fatal, &
362 "calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
363 if (.not. allocated(cs%Res_fn_u)) call mom_error(fatal, &
364 "calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
365 if (.not. allocated(cs%Res_fn_v)) call mom_error(fatal, &
366 "calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
367 if (.not. allocated(cs%f2_dx2_h)) call mom_error(fatal, &
368 "calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
369 if (.not. allocated(cs%f2_dx2_q)) call mom_error(fatal, &
370 "calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
371 if (.not. allocated(cs%f2_dx2_u)) call mom_error(fatal, &
372 "calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
373 if (.not. allocated(cs%f2_dx2_v)) call mom_error(fatal, &
374 "calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
375 if (.not. allocated(cs%beta_dx2_h)) call mom_error(fatal, &
376 "calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
377 if (.not. allocated(cs%beta_dx2_q)) call mom_error(fatal, &
378 "calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
379 if (.not. allocated(cs%beta_dx2_u)) call mom_error(fatal, &
380 "calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
381 if (.not. allocated(cs%beta_dx2_v)) call mom_error(fatal, &
382 "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")
383
384 apply_u_obc = .false. ; apply_v_obc = .false.
385 if (associated(obc) .and. (.not.cs%res_fn_OBC_bug)) then
386 apply_u_obc = obc%u_OBCs_on_PE
387 apply_v_obc = obc%v_OBCs_on_PE
388 endif
389
390 !$OMP parallel default(shared) private(dx_term,power_2)
391
392 if (apply_u_obc .or. apply_v_obc) then
393 !$OMP do
394 do j=js-1,jeq ; do i=is-1,ieq
395 if ((obc%segnum_u(i,j) /= 0) .or. (obc%segnum_u(i,j+1) /= 0) .or. &
396 (obc%segnum_v(i,j) /= 0) .or. (obc%segnum_u(i+1,j) /= 0)) then
397 ! This is an OBC node, so use the fact that G%mask2dT is zero behind OBCs. The nondimensional
398 ! constant 1e-20 in the denominator makes this a de facto implementation of Adcroft's reciprocal
399 ! rule with a value that works for either 64-bit or 32-bit real numbers.
400 cg1_q(i,j) = ((g%mask2dT(i,j) * cs%cg1(i,j) + g%mask2dT(i+1,j+1) * cs%cg1(i+1,j+1)) + &
401 (g%mask2dT(i+1,j) * cs%cg1(i+1,j) + g%mask2dT(i,j+1) * cs%cg1(i,j+1))) / &
402 ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) + 1.0e-20)
403 else
404 cg1_q(i,j) = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
405 endif
406 enddo ; enddo
407 else
408 !$OMP do
409 do j=js-1,jeq ; do i=is-1,ieq
410 cg1_q(i,j) = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
411 enddo ; enddo
412 endif
413
414 ! Do this calculation on the extent used in MOM_hor_visc.F90, and
415 ! MOM_tracer.F90 so that no halo update is needed.
416 if (cs%Res_fn_power_visc >= 100) then
417 !$OMP do
418 do j=js-1,je+1 ; do i=is-1,ie+1
419 dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
420 if ((cs%Res_coef_visc * cs%cg1(i,j))**2 > dx_term) then
421 cs%Res_fn_h(i,j) = 0.0
422 else
423 cs%Res_fn_h(i,j) = 1.0
424 endif
425 enddo ; enddo
426 !$OMP do
427 do j=js-1,jeq ; do i=is-1,ieq
428 dx_term = cs%f2_dx2_q(i,j) + cg1_q(i,j) * cs%beta_dx2_q(i,j)
429 if ((cs%Res_coef_visc * cg1_q(i,j))**2 > dx_term) then
430 cs%Res_fn_q(i,j) = 0.0
431 else
432 cs%Res_fn_q(i,j) = 1.0
433 endif
434 enddo ; enddo
435 elseif (cs%Res_fn_power_visc == 2) then
436 !$OMP do
437 do j=js-1,je+1 ; do i=is-1,ie+1
438 dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
439 cs%Res_fn_h(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**2)
440 enddo ; enddo
441 !$OMP do
442 do j=js-1,jeq ; do i=is-1,ieq
443 dx_term = cs%f2_dx2_q(i,j) + cg1_q(i,j) * cs%beta_dx2_q(i,j)
444 cs%Res_fn_q(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cg1_q(i,j))**2)
445 enddo ; enddo
446 elseif (mod(cs%Res_fn_power_visc, 2) == 0) then
447 power_2 = cs%Res_fn_power_visc / 2
448 !$OMP do
449 do j=js-1,je+1 ; do i=is-1,ie+1
450 dx_term = (us%L_T_to_m_s**2*(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**power_2
451 cs%Res_fn_h(i,j) = dx_term / &
452 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
453 enddo ; enddo
454 !$OMP do
455 do j=js-1,jeq ; do i=is-1,ieq
456 dx_term = (us%L_T_to_m_s**2*(cs%f2_dx2_q(i,j) + cg1_q(i,j) * cs%beta_dx2_q(i,j)))**power_2
457 cs%Res_fn_q(i,j) = dx_term / &
458 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q(i,j))**cs%Res_fn_power_visc)
459 enddo ; enddo
460 else
461 !$OMP do
462 do j=js-1,je+1 ; do i=is-1,ie+1
463 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_h(i,j) + &
464 cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**cs%Res_fn_power_visc
465 cs%Res_fn_h(i,j) = dx_term / &
466 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
467 enddo ; enddo
468 !$OMP do
469 do j=js-1,jeq ; do i=is-1,ieq
470 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_q(i,j) + &
471 cg1_q(i,j) * cs%beta_dx2_q(i,j)))**cs%Res_fn_power_visc
472 cs%Res_fn_q(i,j) = dx_term / &
473 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q(i,j))**cs%Res_fn_power_visc)
474 enddo ; enddo
475 endif
476
477 if (cs%interpolate_Res_fn) then
478 if (apply_u_obc) then
479 do j=js,je ; do i=is-1,ieq
480 cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
481 if (obc%segnum_u(i,j) > 0) cs%Res_fn_u(i,j) = cs%Res_fn_h(i,j) ! Eastern OBC
482 if (obc%segnum_u(i,j) < 0) cs%Res_fn_u(i,j) = cs%Res_fn_h(i+1,j) ! Western OBC
483 enddo ; enddo
484 else
485 do j=js,je ; do i=is-1,ieq
486 cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
487 enddo ; enddo
488 endif
489
490 if (apply_v_obc) then
491 do j=js-1,jeq ; do i=is,ie
492 cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
493 if (obc%segnum_v(i,j) > 0) cs%Res_fn_v(i,j) = cs%Res_fn_h(i,j) ! Northern OBC
494 if (obc%segnum_v(i,j) < 0) cs%Res_fn_v(i,j) = cs%Res_fn_h(i,j+1) ! Southern OBC
495 enddo ; enddo
496 else
497 do j=js-1,jeq ; do i=is,ie
498 cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
499 enddo ; enddo
500 endif
501
502 else ! .not.CS%interpolate_Res_fn
503 if (apply_u_obc) then
504 !$OMP do
505 do j=js,je ; do i=is-1,ieq
506 cg1_u(i,j) = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
507 if (obc%segnum_u(i,j) > 0) cg1_u(i,j) = cs%cg1(i,j) ! Eastern OBC
508 if (obc%segnum_u(i,j) < 0) cg1_u(i,j) = cs%cg1(i+1,j) ! Western OBC
509 enddo ; enddo
510 else
511 !$OMP do
512 do j=js,je ; do i=is-1,ieq
513 cg1_u(i,j) = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
514 enddo ; enddo
515 endif
516
517 if (apply_v_obc) then
518 !$OMP do
519 do j=js-1,jeq ; do i=is,ie
520 cg1_v(i,j) = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
521 if (obc%segnum_v(i,j) > 0) cg1_v(i,j) = cs%cg1(i,j) ! Northern OBC
522 if (obc%segnum_v(i,j) < 0) cg1_v(i,j) = cs%cg1(i,j+1) ! Southern OBC
523 enddo ; enddo
524 else
525 !$OMP do
526 do j=js-1,jeq ; do i=is,ie
527 cg1_v(i,j) = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
528 enddo ; enddo
529 endif
530
531 if (cs%Res_fn_power_khth >= 100) then
532 !$OMP do
533 do j=js,je ; do i=is-1,ieq
534 dx_term = cs%f2_dx2_u(i,j) + cg1_u(i,j) * cs%beta_dx2_u(i,j)
535 if ((cs%Res_coef_khth * cg1_u(i,j))**2 > dx_term) then
536 cs%Res_fn_u(i,j) = 0.0
537 else
538 cs%Res_fn_u(i,j) = 1.0
539 endif
540 enddo ; enddo
541 !$OMP do
542 do j=js-1,jeq ; do i=is,ie
543 dx_term = cs%f2_dx2_v(i,j) + cg1_v(i,j) * cs%beta_dx2_v(i,j)
544 if ((cs%Res_coef_khth * cg1_v(i,j))**2 > dx_term) then
545 cs%Res_fn_v(i,j) = 0.0
546 else
547 cs%Res_fn_v(i,j) = 1.0
548 endif
549 enddo ; enddo
550 elseif (cs%Res_fn_power_khth == 2) then
551 !$OMP do
552 do j=js,je ; do i=is-1,ieq
553 dx_term = cs%f2_dx2_u(i,j) + cg1_u(i,j) * cs%beta_dx2_u(i,j)
554 cs%Res_fn_u(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_u(i,j))**2)
555 enddo ; enddo
556 !$OMP do
557 do j=js-1,jeq ; do i=is,ie
558 dx_term = cs%f2_dx2_v(i,j) + cg1_v(i,j) * cs%beta_dx2_v(i,j)
559 cs%Res_fn_v(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_v(i,j))**2)
560 enddo ; enddo
561 elseif (mod(cs%Res_fn_power_khth, 2) == 0) then
562 power_2 = cs%Res_fn_power_khth / 2
563 !$OMP do
564 do j=js,je ; do i=is-1,ieq
565 dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_u(i,j) + cg1_u(i,j) * cs%beta_dx2_u(i,j)))**power_2
566 cs%Res_fn_u(i,j) = dx_term / &
567 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u(i,j))**cs%Res_fn_power_khth)
568 enddo ; enddo
569 !$OMP do
570 do j=js-1,jeq ; do i=is,ie
571 dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_v(i,j) + cg1_v(i,j) * cs%beta_dx2_v(i,j)))**power_2
572 cs%Res_fn_v(i,j) = dx_term / &
573 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v(i,j))**cs%Res_fn_power_khth)
574 enddo ; enddo
575 else
576 !$OMP do
577 do j=js,je ; do i=is-1,ieq
578 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_u(i,j) + &
579 cg1_u(i,j) * cs%beta_dx2_u(i,j)))**cs%Res_fn_power_khth
580 cs%Res_fn_u(i,j) = dx_term / &
581 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u(i,j))**cs%Res_fn_power_khth)
582 enddo ; enddo
583 !$OMP do
584 do j=js-1,jeq ; do i=is,ie
585 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_v(i,j) + &
586 cg1_v(i,j) * cs%beta_dx2_v(i,j)))**cs%Res_fn_power_khth
587 cs%Res_fn_v(i,j) = dx_term / &
588 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v(i,j))**cs%Res_fn_power_khth)
589 enddo ; enddo
590 endif
591 endif
592 !$OMP end parallel
593
594 if (query_averaging_enabled(cs%diag)) then
595 if (cs%id_Res_fn > 0) call post_data(cs%id_Res_fn, cs%Res_fn_h, cs%diag)
596 if (cs%id_BS_struct > 0) call post_data(cs%id_BS_struct, cs%BS_struct, cs%diag)
597 if (cs%id_khth_struct > 0) call post_data(cs%id_khth_struct, cs%khth_struct, cs%diag)
598 if (cs%id_khtr_struct > 0) call post_data(cs%id_khtr_struct, cs%khtr_struct, cs%diag)
599 if (cs%id_kdgl90_struct > 0) call post_data(cs%id_kdgl90_struct, cs%kdgl90_struct, cs%diag)
600 endif
601
602 if (cs%debug) then
603 call hchksum(cs%cg1, "calc_resoln_fn cg1", g%HI, haloshift=1, unscale=us%L_T_to_m_s)
604 call uvchksum("Res_fn_[uv]", cs%Res_fn_u, cs%Res_fn_v, g%HI, haloshift=0, &
605 unscale=1.0, scalar_pair=.true.)
606 endif
607
608end subroutine calc_resoln_function
609
610!> Calculates and stores functions of SQG mode
611subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE, OBC)
612 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
613 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
614 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
615 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
616 type(thermo_var_ptrs), intent(in) :: tv !<Thermodynamic variables
617 real, intent(in) :: dt !< Time increment [T ~> s]
618 type(varmix_cs), intent(inout) :: cs !< Variable mixing control struct
619 type(meke_type), intent(in) :: meke !< MEKE struct
620 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure
621
622 ! Local variables
623 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
624 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: n2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2]
625 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: n2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2]
626 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
627 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
628 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzsxn ! |Sx| N times dz at u-points [Z T-1 ~> m s-1]
629 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzsyn ! |Sy| N times dz at v-points [Z T-1 ~> m s-1]
630 real, dimension(SZI_(G),SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1]
631 real :: n2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2]
632 real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [Z ~> m]
633 real :: f_subround ! The minimal resolved value of Coriolis parameter to prevent division by zero [T-1 ~> s-1]
634 real, dimension(SZI_(G),SZJ_(G)) :: le ! Eddy length scale [L ~> m]
635
636 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Geometric layer thicknesses in height units [Z ~> m]
637 real :: i_f_le(szi_(g),szj_(g)) ! The inverse of the absolute value of f times the Eddy
638 ! length scale [T L-1 ~> s m-1]
639 real :: p_i(szi_(g),szj_(g)) ! Pressure at the interface [R L2 T-2 ~> Pa]
640 real :: t_i(szi_(g)) ! Temperature at the interface [C ~> degC]
641 real :: s_i(szi_(g)) ! Salinity at the interface [S ~> ppt]
642 real :: drho_ds(szi_(g)) ! Local change in density with salinity using the model EOS and
643 ! state interpolated to an interface [R C-1 ~> kg m-3 ppt-1]
644 real :: drho_dt(szi_(g)) ! Local change in density with salinity using the model EOS and
645 ! state interpolated [R C-1 ~> kg m-3 degC-1]
646 real :: h_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
647 real :: gxspv ! Gravitiational acceleration times the specific volume at an interface
648 ! [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
649 real :: drdk ! Vertical density differences across an interface [R ~> kg m-3]
650 real :: dz_int ! Average of thicknesses around an interface in height units [Z ~> m]
651 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
652 integer :: i, j, k, is, ie, js, je, nz
653
654 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
655 f_subround = 1.0e-40 * us%s_to_T
656
657 if (.not. cs%initialized) call mom_error(fatal, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
658 "Module must be initialized before it is used.")
659
660 if (cs%sqg_expo <= 0.) then
661 cs%sqg_struct(:,:,:) = 1.
662 else
663 if (allocated(meke%Le)) then
664 do j=js,je ; do i=is,ie
665 le(i,j) = meke%Le(i,j)
666 enddo ; enddo
667 else
668 do j=js,je ; do i=is,ie
669 le(i,j) = sqrt(g%areaT(i,j))
670 enddo ; enddo
671 endif
672
673 do j=js,je ; do i=is,ie
674 ! Setting the structure averaged over the top layer to 1 is consistent with it being well mixed.
675 cs%sqg_struct(i,j,1) = 1.0
676 enddo ; enddo
677
678 if (cs%interpolated_sqg_struct) then
679 do j=js,je ; do i=is,ie
680 f(i,j) = max(0.25 * abs((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
681 (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1))), f_subround)
682 enddo ; enddo
683 call find_eta(h, tv, g, gv, us, e, halo_size=2) !### Could be halo_size=1?
684 call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, cs%use_stanley_iso, &
685 cs%slope_x, cs%slope_y, n2_u=n2_u, n2_v=n2_v, dzu=dzu, dzv=dzv, &
686 dzsxn=dzsxn, dzsyn=dzsyn, halo=1, obc=obc, obc_n2=cs%OBC_friendly)
687 do k=2,nz ; do j=js,je ; do i=is,ie
688 n2 = max(0.25 * ((n2_u(i-1,j,k) + n2_u(i,j,k)) + (n2_v(i,j-1,k) + n2_v(i,j,k))), 0.0)
689 dzc = 0.25 * ((dzu(i-1,j,k) + dzu(i,j,k)) + (dzv(i,j-1,k) + dzv(i,j,k)))
690 cs%sqg_struct(i,j,k) = cs%sqg_struct(i,j,k-1) * &
691 exp(-cs%sqg_expo * (dzc * sqrt(n2)/(f(i,j) * le(i,j))))
692 enddo ; enddo ; enddo
693 else
694 do j=js,je ; do i=is,ie
695 i_f_le(i,j) = 1.0 / &
696 (le(i,j) * max(0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
697 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1)))), f_subround))
698 enddo ; enddo
699
700 call thickness_to_dz(h, tv, dz, g, gv, us)
701
702 if (associated(tv%eqn_of_state)) then
703 eosdom(:) = eos_domain(g%HI)
704 h_to_pres = gv%H_to_RZ * gv%g_Earth
705 ! Set the pressure at the topmost interior interface.
706 p_i(:,:) = 0.0
707 if (associated(tv%p_surf)) then
708 do j=js,je ; do i=is,ie ; p_i(i,j) = tv%p_surf(i,j) ; enddo ; enddo
709 endif
710 if (.not.allocated(tv%SpV_avg)) gxspv = gv%g_Earth / gv%Rho0
711 do k=2,nz ; do j=js,je
712 ! Find the derivatives of density with T and S at the interface.
713 do i=is,ie
714 p_i(i,j) = p_i(i,j) + h_to_pres * h(i,j,k-1)
715 t_i(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
716 s_i(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
717 enddo
718 call calculate_density_derivs(t_i, s_i, p_i(:,j), drho_dt, drho_ds, tv%eqn_of_state, eosdom)
719
720 do i=is,ie
721 if (allocated(tv%SpV_avg)) & ! GxSpV is in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
722 gxspv = gv%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
723
724 drdk = max(drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
725 drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1)), 0.0) ! Density difference [R ~> kg m-3]
726 dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
727 cs%sqg_struct(i,j,k) = cs%sqg_struct(i,j,k-1) * &
728 exp(-cs%sqg_expo * (sqrt((gxspv * drdk) * dz_int) * i_f_le(i,j)) )
729 ! To derive the expression above, note that
730 ! N2 = GxSpV * drdk / dzh(i,j,K) ! Square of positive buoyancy freq. [L2 Z-2 T-2 ~> s-2]
731 ! CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
732 ! exp(-CS%sqg_expo * (dz_int(i,j,K) * sqrt(N2) * I_f_Le(i,j)) )
733 enddo
734 enddo ; enddo
735 else ! (GV%Boussinesq .and. .not.use_EOS) then
736 do k=2,nz ; do j=js,je ; do i=is,ie
737 dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
738 cs%sqg_struct(i,j,k) = cs%sqg_struct(i,j,k-1) * &
739 exp(-cs%sqg_expo * (sqrt(gv%g_prime(k) * dz_int) * i_f_le(i,j)) )
740 enddo ; enddo ; enddo
741 endif
742 endif
743 endif
744
745 if (query_averaging_enabled(cs%diag)) then
746 if (cs%id_sqg_struct > 0) call post_data(cs%id_sqg_struct, cs%sqg_struct, cs%diag)
747 if (cs%interpolated_sqg_struct .and. (cs%sqg_expo > 0.)) then
748 if (cs%id_N2_u > 0) call post_data(cs%id_N2_u, n2_u, cs%diag)
749 if (cs%id_N2_v > 0) call post_data(cs%id_N2_v, n2_v, cs%diag)
750 endif
751 endif
752
753end subroutine calc_sqg_struct
754
755!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
756!! style scaling of diffusivity
757subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
758 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
759 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
760 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
761 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
762 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
763 real, intent(in) :: dt !< Time increment [T ~> s]
764 type(varmix_cs), intent(inout) :: cs !< Variable mixing control structure
765 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure
766
767 ! Local variables
768 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
769 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: n2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2]
770 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: n2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2]
771 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
772 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
773 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzsxn ! |Sx| N times dz at u-points [Z T-1 ~> m s-1]
774 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzsyn ! |Sy| N times dz at v-points [Z T-1 ~> m s-1]
775
776 if (.not. cs%initialized) call mom_error(fatal, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
777 "Module must be initialized before it is used.")
778
779 if (cs%calculate_Eady_growth_rate) then
780 call find_eta(h, tv, g, gv, us, e, halo_size=2)
781 if (cs%use_simpler_Eady_growth_rate) then
782 call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, cs%use_stanley_iso, &
783 cs%slope_x, cs%slope_y, n2_u=n2_u, n2_v=n2_v, dzu=dzu, dzv=dzv, &
784 dzsxn=dzsxn, dzsyn=dzsyn, halo=1, obc=obc, obc_n2=cs%OBC_friendly)
785 call calc_eady_growth_rate_2d(cs, g, gv, us, h, e, dzu, dzv, dzsxn, dzsyn, cs%SN_u, cs%SN_v)
786 elseif (cs%use_stored_slopes) then
787 call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, cs%use_stanley_iso, &
788 cs%slope_x, cs%slope_y, n2_u=n2_u, n2_v=n2_v, halo=1, obc=obc, &
789 obc_n2=cs%OBC_friendly)
790 call calc_visbeck_coeffs_old(h, cs%slope_x, cs%slope_y, n2_u, n2_v, g, gv, us, cs, obc)
791 else
792 call calc_slope_functions_using_just_e(h, g, gv, us, cs, e)
793 endif
794 endif
795
796 if (query_averaging_enabled(cs%diag)) then
797 if (cs%id_dzu > 0) call post_data(cs%id_dzu, dzu, cs%diag)
798 if (cs%id_dzv > 0) call post_data(cs%id_dzv, dzv, cs%diag)
799 if (cs%id_dzSxN > 0) call post_data(cs%id_dzSxN, dzsxn, cs%diag)
800 if (cs%id_dzSyN > 0) call post_data(cs%id_dzSyN, dzsyn, cs%diag)
801 if (cs%id_SN_u > 0) call post_data(cs%id_SN_u, cs%SN_u, cs%diag)
802 if (cs%id_SN_v > 0) call post_data(cs%id_SN_v, cs%SN_v, cs%diag)
803 if (cs%id_L2u > 0) call post_data(cs%id_L2u, cs%L2u, cs%diag)
804 if (cs%id_L2v > 0) call post_data(cs%id_L2v, cs%L2v, cs%diag)
805 if (cs%id_N2_u > 0) call post_data(cs%id_N2_u, n2_u, cs%diag)
806 if (cs%id_N2_v > 0) call post_data(cs%id_N2_v, n2_v, cs%diag)
807 endif
808
809end subroutine calc_slope_functions
810
811!> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al., 1997.
812!! This is on older implementation that is susceptible to large values of Eady growth rate
813!! for incropping layers.
814subroutine calc_visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
815 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
816 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
817 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
818 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [Z L-1 ~> nondim]
819 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency
820 !! at u-points [L2 Z-2 T-2 ~> s-2]
821 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope
822 !! [Z L-1 ~> nondim]
823 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency
824 !! at v-points [L2 Z-2 T-2 ~> s-2]
825 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
826 type(varmix_cs), intent(inout) :: CS !< Variable mixing control structure
827 type(ocean_obc_type), pointer :: OBC !< Open boundaries control structure.
828
829 ! Local variables
830 real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
831 real :: N2 ! Positive buoyancy frequency or zero [L2 Z-2 T-2 ~> s-2]
832 real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
833 real :: H_geom ! The geometric mean of Hup and Hdn [H ~> m or kg m-2].
834 real :: S2max ! An upper bound on the squared slopes [Z2 L-2 ~> nondim]
835 real :: wNE, wSE, wSW, wNW ! Weights of adjacent points [nondim]
836 real :: H_u(SZIB_(G)), H_v(SZI_(G)) ! Layer thicknesses at u- and v-points [H ~> m or kg m-2]
837
838 ! Note that at some points in the code S2_u and S2_v hold the running depth
839 ! integrals of the squared slope [H ~> m or kg m-2] before the average is taken.
840 real :: S2_u(SZIB_(G),SZJ_(G)) ! At first the thickness-weighted depth integral of the squared
841 ! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
842 ! squared slope [Z2 L-2 ~> nondim] at u points.
843 real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared
844 ! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
845 ! squared slope [Z2 L-2 ~> nondim] at v points.
846 integer :: OBC_dir_u(SZIB_(G),SZJ_(G)) ! An integer indicating where there are u OBCs: +1 for
847 ! eastern OBCs, -1 for western OBCs and 0 at points with no OBCs.
848 integer :: OBC_dir_v(SZI_(G),SZJB_(G)) ! An integer indicating where there are v OBCs: +1 for
849 ! northern OBCs, -1 for southern OBCs and 0 at points with no OBCs.
850 real :: h4_u(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a u-point
851 ! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
852 real :: h4_v(SZI_(G),SZJB_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a v-point
853 ! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
854 integer :: i, j, k, is, ie, js, je, nz
855
856 if (.not. cs%initialized) call mom_error(fatal, "calc_Visbeck_coeffs_old: "// &
857 "Module must be initialized before it is used.")
858
859 if (.not. cs%calculate_Eady_growth_rate) return
860 if (.not. allocated(cs%SN_u)) call mom_error(fatal, "calc_slope_function: "// &
861 "%SN_u is not associated with use_variable_mixing.")
862 if (.not. allocated(cs%SN_v)) call mom_error(fatal, "calc_slope_function: "// &
863 "%SN_v is not associated with use_variable_mixing.")
864
865 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
866
867 s2max = cs%Visbeck_S_max**2
868
869 cs%SN_u(:,:) = 0.0
870 cs%SN_v(:,:) = 0.0
871
872 ! These settings apply where there are not open boundary conditions.
873 obc_dir_u(:,:) = 0 ; obc_dir_v(:,:) = 0
874
875 if (associated(obc) .and. cs%OBC_friendly) then
876 ! Store the direction of any OBC faces.
877 !$OMP parallel do default(shared)
878 do j=js-1,je+1 ; do i=is-1,ie ; if (obc%segnum_u(i,j) /= 0) then
879 if (obc%segnum_u(i,j) > 0) obc_dir_u(i,j) = 1 ! OBC_DIRECTION_E
880 if (obc%segnum_u(i,j) < 0) obc_dir_u(i,j) = -1 ! OBC_DIRECTION_W
881 endif ; enddo ; enddo
882 !$OMP parallel do default(shared)
883 do j=js-1,je ; do i=is-1,ie+1 ; if (obc%segnum_v(i,j) /= 0) then
884 if (obc%segnum_v(i,j) > 0) obc_dir_v(i,j) = 1 ! OBC_DIRECTION_N
885 if (obc%segnum_v(i,j) < 0) obc_dir_v(i,j) = -1 ! OBC_DIRECTION_S
886 endif ; enddo ; enddo
887
888 ! Use the masked product of the 4 (or 2) thicknesses around a velocity-point interface for weights.
889 !$OMP parallel do default(shared)
890 do k=2,nz
891 do j=js-1,je+1 ; do i=is-1,ie
892 if (obc_dir_u(i,j) == 0) then
893 h4_u(i,j,k) = g%mask2dCu(i,j) * ( (h(i,j,k)*h(i+1,j,k)) * (h(i,j,k-1)*h(i+1,j,k-1)) )
894 elseif (obc_dir_u(i,j) == 1) then ! OBC_DIRECTION_E
895 h4_u(i,j,k) = g%mask2dCu(i,j) * ( (h(i,j,k)**2) * (h(i,j,k-1)**2) )
896 elseif (obc_dir_u(i,j) == -1) then ! OBC_DIRECTION_W
897 h4_u(i,j,k) = g%mask2dCu(i,j) * ( (h(i+1,j,k)**2) * (h(i+1,j,k-1)**2) )
898 endif
899 enddo ; enddo
900 do j=js-1,je ; do i=is-1,ie+1
901 if (obc_dir_v(i,j) == 0) then
902 h4_v(i,j,k) = g%mask2dCv(i,j) * ( (h(i,j,k)*h(i,j+1,k)) * (h(i,j,k-1)*h(i,j+1,k-1)) )
903 elseif (obc_dir_v(i,j) == 1) then ! OBC_DIRECTION_N
904 h4_v(i,j,k) = g%mask2dCv(i,j) * ( (h(i,j,k)**2) * (h(i,j,k-1)**2) )
905 elseif (obc_dir_v(i,j) == -1) then ! OBC_DIRECTION_S
906 h4_v(i,j,k) = g%mask2dCv(i,j) * ( (h(i,j+1,k)**2) * (h(i,j+1,k-1)**2) )
907 endif
908 enddo ; enddo
909 enddo
910 else ! The land mask is sufficient and there are no special considerations taken at OBC points.
911 ! Use the masked product of the 4 thicknesses around a velocity-point interface for weights.
912 !$OMP parallel do default(shared)
913 do k=2,nz
914 do j=js-1,je+1 ; do i=is-1,ie
915 h4_u(i,j,k) = g%mask2dCu(i,j) * ( (h(i,j,k)*h(i+1,j,k)) * (h(i,j,k-1)*h(i+1,j,k-1)) )
916 enddo ; enddo
917 do j=js-1,je ; do i=is-1,ie+1
918 h4_v(i,j,k) = g%mask2dCv(i,j) * ( (h(i,j,k)*h(i,j+1,k)) * (h(i,j,k-1)*h(i,j+1,k-1)) )
919 enddo ; enddo
920 enddo
921 endif
922
923 ! To set the length scale based on the deformation radius, use wave_speed to
924 ! calculate the first-mode gravity wave speed and then blend the equatorial
925 ! and midlatitude deformation radii, using calc_resoln_function as a template.
926
927 !$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
928 do j=js,je
929 do i=is-1,ie
930 cs%SN_u(i,j) = 0. ; h_u(i) = 0. ; s2_u(i,j) = 0.
931 enddo
932 do k=2,nz ; do i=is-1,ie
933 hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
934 hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
935 h_geom = sqrt( hdn * hup )
936 !H_geom = H_geom * sqrt(N2) ! WKB-ish
937 !H_geom = H_geom * N2 ! WKB-ish
938 wse = h4_v(i+1,j-1,k)
939 wnw = h4_v(i,j,k)
940 wne = h4_v(i+1,j,k)
941 wsw = h4_v(i,j-1,k)
942 if (obc_dir_u(i,j) == 1) then ! OBC_DIRECTION_E
943 wse = 0.0 ; wne = 0.0
944 h_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
945 elseif (obc_dir_u(i,j) == -1) then ! OBC_DIRECTION_W
946 wsw = 0.0 ; wnw = 0.0
947 h_geom = sqrt( h(i+1,j,k) * h(i+1,j,k-1) )
948 endif
949 s2 = slope_x(i,j,k)**2 + &
950 (((wnw*slope_y(i,j,k)**2) + (wse*slope_y(i+1,j-1,k)**2)) + &
951 ((wne*slope_y(i+1,j,k)**2) + (wsw*slope_y(i,j-1,k)**2)) ) / &
952 ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
953 if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
954
955 n2 = max(0., n2_u(i,j,k))
956 cs%SN_u(i,j) = cs%SN_u(i,j) + sqrt( s2*n2 )*h_geom
957 s2_u(i,j) = s2_u(i,j) + s2*h_geom
958 h_u(i) = h_u(i) + h_geom
959 enddo ; enddo
960 do i=is-1,ie
961 if (h_u(i)>0.) then
962 cs%SN_u(i,j) = g%OBCmaskCu(i,j) * cs%SN_u(i,j) / h_u(i)
963 s2_u(i,j) = g%OBCmaskCu(i,j) * s2_u(i,j) / h_u(i)
964 else
965 cs%SN_u(i,j) = 0.
966 endif
967 enddo
968 enddo
969
970 !$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
971 do j=js-1,je
972 do i=is,ie
973 cs%SN_v(i,j) = 0. ; h_v(i) = 0. ; s2_v(i,j) = 0.
974 enddo
975 do k=2,nz ; do i=is,ie
976 hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
977 hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
978 h_geom = sqrt( hdn * hup )
979 !H_geom = H_geom * sqrt(N2) ! WKB-ish
980 !H_geom = H_geom * N2 ! WKB-ish
981 wse = h4_u(i,j,k)
982 wnw = h4_u(i-1,j+1,k)
983 wne = h4_u(i,j+1,k)
984 wsw = h4_u(i-1,j,k)
985 if (obc_dir_v(i,j) == 1) then ! OBC_DIRECTION_N
986 wnw = 0.0 ; wne = 0.0
987 h_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
988 elseif (obc_dir_v(i,j) == -1) then ! OBC_DIRECTION_S
989 wsw = 0.0 ; wse = 0.0
990 h_geom = sqrt( h(i,j+1,k) * h(i,j+1,k-1) )
991 endif
992 s2 = slope_y(i,j,k)**2 + &
993 (((wse*slope_x(i,j,k)**2) + (wnw*slope_x(i-1,j+1,k)**2)) + &
994 ((wne*slope_x(i,j+1,k)**2) + (wsw*slope_x(i-1,j,k)**2)) ) / &
995 ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
996 if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
997
998 n2 = max(0., n2_v(i,j,k))
999 cs%SN_v(i,j) = cs%SN_v(i,j) + sqrt( s2*n2 )*h_geom
1000 s2_v(i,j) = s2_v(i,j) + s2*h_geom
1001 h_v(i) = h_v(i) + h_geom
1002 enddo ; enddo
1003 do i=is,ie
1004 if (h_v(i)>0.) then
1005 cs%SN_v(i,j) = g%OBCmaskCv(i,j) * cs%SN_v(i,j) / h_v(i)
1006 s2_v(i,j) = g%OBCmaskCv(i,j) * s2_v(i,j) / h_v(i)
1007 else
1008 cs%SN_v(i,j) = 0.
1009 endif
1010 enddo
1011 enddo
1012
1013 ! Offer diagnostic fields for averaging.
1014 if (query_averaging_enabled(cs%diag)) then
1015 if (cs%id_S2_u > 0) call post_data(cs%id_S2_u, s2_u, cs%diag)
1016 if (cs%id_S2_v > 0) call post_data(cs%id_S2_v, s2_v, cs%diag)
1017 endif
1018
1019 if (cs%debug) then
1020 call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, g%HI, &
1021 unscale=us%Z_to_L, haloshift=1)
1022 ! call uvchksum("calc_Visbeck_coeffs_old S2_[uv]", S2_u, S2_v, G%HI, &
1023 ! unscale=US%Z_to_L**2, scalar_pair=.true.)
1024 call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", n2_u, n2_v, g%HI, &
1025 unscale=us%L_to_Z**2*us%s_to_T**2, scalar_pair=.true.)
1026 call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", cs%SN_u, cs%SN_v, g%HI, &
1027 unscale=us%s_to_T, scalar_pair=.true.)
1028 endif
1029
1030end subroutine calc_visbeck_coeffs_old
1031
1032!> Calculates the Eady growth rate (2D fields) for use in MEKE and the Visbeck schemes
1033subroutine calc_eady_growth_rate_2d(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, SN_u, SN_v)
1034 type(varmix_cs), intent(inout) :: CS !< Variable mixing coefficients
1035 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1036 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1037 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1038 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m]
1039 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m]
1040 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m]
1041 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzv !< dz at v-points [Z ~> m]
1042 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzSxN !< dz Sx N at u-points [Z T-1 ~> m s-1]
1043 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzSyN !< dz Sy N at v-points [Z T-1 ~> m s-1]
1044 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: SN_u !< SN at u-points [T-1 ~> s-1]
1045 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: SN_v !< SN at v-points [T-1 ~> s-1]
1046 ! Local variables
1047 real :: D_scale ! The depth over which to average SN [Z ~> m]
1048 real :: dnew ! Depth of bottom of layer [Z ~> m]
1049 real :: dz ! Limited thickness of this layer [Z ~> m]
1050 real :: weight ! Fraction of this layer that contributes to integral [nondim]
1051 real :: sum_dz(SZI_(G)) ! Cumulative sum of z-thicknesses [Z ~> m]
1052 real :: vint_SN(SZIB_(G)) ! Cumulative integral of SN [Z T-1 ~> m s-1]
1053 real, dimension(SZIB_(G),SZJ_(G)) :: SN_cpy !< SN at u-points [T-1 ~> s-1]
1054 real :: dz_neglect ! A negligibly small distance to avoid division by zero [Z ~> m]
1055 real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~> m-1]
1056 real :: dB, dT ! Elevation variables used when cropping [Z ~> m]
1057 integer :: i, j, k
1058 logical :: crop
1059
1060 dz_neglect = gv%dZ_subroundoff
1061 d_scale = cs%Eady_GR_D_scale
1062 if (d_scale<=0.) d_scale = 64.*gv%max_depth ! 0 means use full depth so choose something big
1063 r_crp_dist = 1. / max( dz_neglect, cs%cropping_distance )
1064 crop = cs%cropping_distance>=0. ! Only filter out in-/out-cropped interface is parameter if non-negative
1065
1066 if (cs%debug) then
1067 call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, g%HI, unscale=us%Z_to_m, scalar_pair=.true.)
1068 call uvchksum("calc_Eady_growth_rate_2D dzS2N2[uv]", dzsxn, dzsyn, g%HI, &
1069 unscale=us%Z_to_m*us%s_to_T, scalar_pair=.true.)
1070 endif
1071
1072 !$OMP parallel do default(shared)
1073 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
1074 cs%SN_u(i,j) = 0.0
1075 cs%SN_v(i,j) = 0.0
1076 enddo ; enddo
1077
1078 !$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1079 do j=g%jsc-1,g%jec+1
1080 do i=g%isc-1,g%iec
1081 vint_sn(i) = 0.
1082 sum_dz(i) = dz_neglect
1083 enddo
1084 if (crop) then
1085 do k=2,gv%ke ; do i=g%isc-1,g%iec
1086 dnew = sum_dz(i) + dzu(i,j,k) ! This is where the bottom of the layer is
1087 dnew = min(dnew, d_scale) ! This limits the depth to D_scale
1088 dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1089 ! When D_scale>dnew, dz=dzu (+roundoff error).
1090 ! When sum_dz<D_scale<dnew, 0<dz<dzu.
1091 ! When D_scale<sum_dz, dz=0.
1092 weight = dz / ( dzu(i,j,k) + dz_neglect ) ! Fraction of this layer to include
1093 dt = min( e(i,j,1), e(i+1,j,1) ) ! Deepest sea surface
1094 db = max( e(i,j,k), e(i+1,j,k) ) ! Shallowest interface
1095 weight = weight * min( max( 0., (dt-db)*r_crp_dist ), 1. )
1096 dt = min( e(i,j,k), e(i+1,j,k) ) ! Deepest interface
1097 db = max( e(i,j,gv%ke+1), e(i+1,j,gv%ke+1) ) ! Shallowest topography
1098 weight = weight * min( max( 0., (dt-db)*r_crp_dist ), 1. )
1099 vint_sn(i) = vint_sn(i) + weight * dzsxn(i,j,k)
1100 sum_dz(i) = sum_dz(i) + weight * dzu(i,j,k)
1101 enddo ; enddo
1102 else
1103 do k=2,gv%ke ; do i=g%isc-1,g%iec
1104 dnew = sum_dz(i) + dzu(i,j,k) ! This is where the bottom of the layer is
1105 dnew = min(dnew, d_scale) ! This limits the depth to D_scale
1106 dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1107 ! When D_scale>dnew, dz=dzu (+roundoff error).
1108 ! When sum_dz<D_scale<dnew, 0<dz<dzu.
1109 ! When D_scale<sum_dz, dz=0.
1110 weight = dz / ( dzu(i,j,k) + dz_neglect ) ! Fraction of this layer to include
1111 vint_sn(i) = vint_sn(i) + weight * dzsxn(i,j,k)
1112 sum_dz(i) = sum_dz(i) + weight * dzu(i,j,k)
1113 enddo ; enddo
1114 endif
1115 do i=g%isc-1,g%iec
1116 cs%SN_u(i,j) = g%OBCmaskCu(i,j) * ( vint_sn(i) / sum_dz(i) )
1117 sn_cpy(i,j) = g%OBCmaskCu(i,j) * ( vint_sn(i) / sum_dz(i) )
1118 enddo
1119 enddo
1120
1121 !$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1122 do j=g%jsc-1,g%jec
1123 do i=g%isc-1,g%iec+1
1124 vint_sn(i) = 0.
1125 sum_dz(i) = dz_neglect
1126 enddo
1127 if (crop) then
1128 do k=2,gv%ke ; do i=g%isc-1,g%iec+1
1129 dnew = sum_dz(i) + dzv(i,j,k) ! This is where the bottom of the layer is
1130 dnew = min(dnew, d_scale) ! This limits the depth to D_scale
1131 dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1132 ! When D_scale>dnew, dz=dzu (+roundoff error).
1133 ! When sum_dz<D_scale<dnew, 0<dz<dzu.
1134 ! When D_scale<sum_dz, dz=0.
1135 weight = dz / ( dzv(i,j,k) + dz_neglect ) ! Fraction of this layer to include
1136 dt = min( e(i,j,1), e(i,j+1,1) ) ! Deepest sea surface
1137 db = max( e(i,j,k), e(i,j+1,k) ) ! Shallowest interface
1138 weight = weight * min( max( 0., (dt-db)*r_crp_dist ), 1. )
1139 dt = min( e(i,j,k), e(i,j+1,k) )! Deepest interface
1140 db = max( e(i,j,gv%ke+1), e(i,j+1,gv%ke+1) ) ! Shallowest topography
1141 weight = weight * min( max( 0., (dt-db)*r_crp_dist ), 1. )
1142 vint_sn(i) = vint_sn(i) + weight**2 * dzsyn(i,j,k)
1143 sum_dz(i) = sum_dz(i) + weight * dzv(i,j,k)
1144 enddo ; enddo
1145 else
1146 do k=2,gv%ke ; do i=g%isc-1,g%iec+1
1147 dnew = sum_dz(i) + dzv(i,j,k) ! This is where the bottom of the layer is
1148 dnew = min(dnew, d_scale) ! This limits the depth to D_scale
1149 dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1150 ! When D_scale>dnew, dz=dzu (+roundoff error).
1151 ! When sum_dz<D_scale<dnew, 0<dz<dzu.
1152 ! When D_scale<sum_dz, dz=0.
1153 weight = dz / ( dzv(i,j,k) + dz_neglect ) ! Fraction of this layer to include
1154 vint_sn(i) = vint_sn(i) + weight**2 * dzsyn(i,j,k)
1155 sum_dz(i) = sum_dz(i) + weight * dzv(i,j,k)
1156 enddo ; enddo
1157 endif
1158 do i=g%isc-1,g%iec+1
1159 cs%SN_v(i,j) = g%OBCmaskCv(i,j) * ( vint_sn(i) / sum_dz(i) )
1160 enddo
1161 enddo
1162
1163 do j=g%jsc,g%jec
1164 do i=g%isc-1,g%iec
1165 cs%SN_u(i,j) = sqrt( sn_cpy(i,j)**2 &
1166 + 0.25*( ((cs%SN_v(i,j)**2) + (cs%SN_v(i+1,j-1)**2)) &
1167 + ((cs%SN_v(i+1,j)**2) + (cs%SN_v(i,j-1)**2)) ) )
1168 enddo
1169 enddo
1170 do j=g%jsc-1,g%jec
1171 do i=g%isc,g%iec
1172 cs%SN_v(i,j) = sqrt( cs%SN_v(i,j)**2 &
1173 + 0.25*( ((sn_cpy(i,j)**2) + (sn_cpy(i-1,j+1)**2)) &
1174 + ((sn_cpy(i,j+1)**2) + (sn_cpy(i-1,j)**2)) ) )
1175 enddo
1176 enddo
1177
1178 if (cs%debug) then
1179 call uvchksum("calc_Eady_growth_rate_2D SN_[uv]", cs%SN_u, cs%SN_v, g%HI, &
1180 unscale=us%s_to_T, scalar_pair=.true.)
1181 endif
1182
1183end subroutine calc_eady_growth_rate_2d
1184
1185!> The original calc_slope_function() that calculated slopes using
1186!! interface positions only, not accounting for density variations.
1187subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
1188 type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1189 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1190 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1191 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1192 type(varmix_cs), intent(inout) :: CS !< Variable mixing control structure
1193 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m]
1194 ! type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1195 ! Local variables
1196 real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics)
1197 real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics)
1198 real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m]
1199 ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m]
1200 real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
1201 real :: dZ_cutoff ! A minimum water column depth for masking [H ~> m or kg m-2]
1202 real :: h1, h2 ! Temporary total thicknesses [Z ~> m]
1203 real :: h_neglect ! A thickness that is so small it is usually lost
1204 ! in roundoff and can be neglected [H ~> m or kg m-2].
1205 real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
1206 real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2]
1207 real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
1208 real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
1209 real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times
1210 ! the buoyancy frequency squared at u-points [Z T-2 ~> m s-2]
1211 real :: S2N2_v_local(SZI_(G),SZJB_(G),SZK_(GV)) ! The depth integral of the slope times
1212 ! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2]
1213 logical :: use_dztot ! If true, use the total water column thickness rather than the
1214 ! bathymetric depth for certain calculations.
1215 integer :: is, ie, js, je, nz
1216 integer :: i, j, k
1217
1218 if (.not. cs%initialized) call mom_error(fatal, "calc_slope_functions_using_just_e: "// &
1219 "Module must be initialized before it is used.")
1220
1221 if (.not. cs%calculate_Eady_growth_rate) return
1222 if (.not. allocated(cs%SN_u)) call mom_error(fatal, "calc_slope_function: "// &
1223 "%SN_u is not associated with use_variable_mixing.")
1224 if (.not. allocated(cs%SN_v)) call mom_error(fatal, "calc_slope_function: "// &
1225 "%SN_v is not associated with use_variable_mixing.")
1226
1227 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1228
1229 h_neglect = gv%H_subroundoff
1230 h_cutoff = real(2*nz) * (gv%Angstrom_H + h_neglect)
1231 dz_cutoff = real(2*nz) * (gv%Angstrom_Z + gv%dz_subroundoff)
1232
1233 use_dztot = cs%full_depth_Eady_growth_rate ! .or. .not.(GV%Boussinesq or GV%semi_Boussinesq)
1234
1235 if (use_dztot) then
1236 !$OMP parallel do default(shared)
1237 do j=js-1,je+1 ; do i=is-1,ie+1
1238 dz_tot(i,j) = e(i,j,1) - e(i,j,nz+1)
1239 enddo ; enddo
1240 ! The following mathematically equivalent expression is more expensive but is less
1241 ! sensitive to roundoff for large Z_ref:
1242 ! call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
1243 ! do j=js-1,je+1
1244 ! do i=is-1,ie+1 ; dz_tot(i,j) = 0.0 ; enddo
1245 ! do k=1,nz ; do i=is-1,ie+1
1246 ! dz_tot(i,j) = dz_tot(i,j) + dz(i,j,k)
1247 ! enddo ; enddo
1248 ! enddo
1249 endif
1250
1251 ! To set the length scale based on the deformation radius, use wave_speed to
1252 ! calculate the first-mode gravity wave speed and then blend the equatorial
1253 ! and midlatitude deformation radii, using calc_resoln_function as a template.
1254
1255 !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
1256 do k=nz,cs%VarMix_Ktop,-1
1257
1258 ! Calculate the interface slopes E_x and E_y and u- and v- points respectively
1259 do j=js-1,je+1 ; do i=is-1,ie
1260 e_x(i,j) = (e(i+1,j,k)-e(i,j,k))*g%IdxCu(i,j)
1261 ! Mask slopes where interface intersects topography
1262 if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
1263 enddo ; enddo
1264 do j=js-1,je ; do i=is-1,ie+1
1265 e_y(i,j) = (e(i,j+1,k)-e(i,j,k))*g%IdyCv(i,j)
1266 ! Mask slopes where interface intersects topography
1267 if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
1268 enddo ; enddo
1269
1270 ! Calculate N*S*h from this layer and add to the sum
1271 do j=js,je ; do i=is-1,ie
1272 s2 = ( e_x(i,j)**2 + 0.25*( &
1273 ((e_y(i,j)**2) + (e_y(i+1,j-1)**2)) + ((e_y(i+1,j)**2) + (e_y(i,j-1)**2)) ) )
1274 if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < h_cutoff) s2 = 0.0
1275
1276 hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1277 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)
1278 h_geom = sqrt(hdn*hup)
1279 ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1280 s2n2_u_local(i,j,k) = (h_geom * s2) * (gv%g_prime(k) / max(hdn, hup, cs%h_min_N2) )
1281 enddo ; enddo
1282 do j=js-1,je ; do i=is,ie
1283 s2 = ( e_y(i,j)**2 + 0.25*( &
1284 ((e_x(i,j)**2) + (e_x(i-1,j+1)**2)) + ((e_x(i,j+1)**2) + (e_x(i-1,j)**2)) ) )
1285 if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < h_cutoff) s2 = 0.0
1286
1287 hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1288 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)
1289 h_geom = sqrt(hdn*hup)
1290 ! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1291 s2n2_v_local(i,j,k) = (h_geom * s2) * (gv%g_prime(k) / (max(hdn, hup, cs%h_min_N2)))
1292 enddo ; enddo
1293
1294 enddo ! k
1295
1296 !$OMP parallel do default(shared)
1297 do j=js,je
1298 do i=is-1,ie ; cs%SN_u(i,j) = 0.0 ; enddo
1299 do k=nz,cs%VarMix_Ktop,-1 ; do i=is-1,ie
1300 cs%SN_u(i,j) = cs%SN_u(i,j) + s2n2_u_local(i,j,k)
1301 enddo ; enddo
1302 ! SN above contains S^2*N^2*H, convert to vertical average of S*N
1303
1304 if (use_dztot) then
1305 do i=is-1,ie
1306 cs%SN_u(i,j) = g%OBCmaskCu(i,j) * sqrt( cs%SN_u(i,j) / &
1307 max(dz_tot(i,j), dz_tot(i+1,j), gv%dz_subroundoff) )
1308 enddo
1309 else
1310 do i=is-1,ie
1311 h1 = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0)
1312 h2 = max(g%meanSL(i+1,j) + g%bathyT(i+1,j), 0.0)
1313 if ( min(h1, h2) > dz_cutoff ) then
1314 cs%SN_u(i,j) = g%OBCmaskCu(i,j) * sqrt( cs%SN_u(i,j) / max(h1, h2) )
1315 else
1316 cs%SN_u(i,j) = 0.0
1317 endif
1318 enddo
1319 endif
1320 enddo
1321 !$OMP parallel do default(shared)
1322 do j=js-1,je
1323 do i=is,ie ; cs%SN_v(i,j) = 0.0 ; enddo
1324 do k=nz,cs%VarMix_Ktop,-1 ; do i=is,ie
1325 cs%SN_v(i,j) = cs%SN_v(i,j) + s2n2_v_local(i,j,k)
1326 enddo ; enddo
1327 if (use_dztot) then
1328 do i=is,ie
1329 cs%SN_v(i,j) = g%OBCmaskCv(i,j) * sqrt( cs%SN_v(i,j) / &
1330 max(dz_tot(i,j), dz_tot(i,j+1), gv%dz_subroundoff) )
1331 enddo
1332 else
1333 do i=is,ie
1334 ! There is a primordial horizontal indexing bug on the following line from the previous
1335 ! versions of the code. This comment should be deleted by the end of 2024.
1336 ! if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
1337 h1 = max(g%meanSL(i,j) + g%bathyT(i,j), 0.0)
1338 h2 = max(g%meanSL(i,j+1) + g%bathyT(i,j+1), 0.0)
1339 if ( min(h1, h2) > dz_cutoff ) then
1340 cs%SN_v(i,j) = g%OBCmaskCv(i,j) * sqrt( cs%SN_v(i,j) / max(h1, h2) )
1341 else
1342 cs%SN_v(i,j) = 0.0
1343 endif
1344 enddo
1345 endif
1346 enddo
1347
1349
1350
1351!> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity.
1352subroutine calc_qg_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC)
1353 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1354 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1355 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1356 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1357 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1358 real, intent(in) :: dt !< Time increment [T ~> s]
1359 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1360 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1361 type(varmix_cs), intent(in) :: cs !< Variable mixing control structure
1362 type(ocean_obc_type), pointer :: obc !< Open boundaries control structure
1363 ! Local variables
1364 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
1365
1366 if (.not. cs%initialized) call mom_error(fatal, "MOM_lateral_mixing_coeffs.F90, calc_QG_slopes: "//&
1367 "Module must be initialized before it is used.")
1368
1369 call find_eta(h, tv, g, gv, us, e, halo_size=3)
1370 call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, cs%use_stanley_iso, &
1371 slope_x, slope_y, halo=2, obc=obc, obc_n2=cs%OBC_friendly)
1372
1373end subroutine calc_qg_slopes
1374
1375!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
1376subroutine calc_qg_leith_viscosity(CS, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, slope_x, slope_y, &
1377 vort_xy_dx, vort_xy_dy)
1378 type(varmix_cs), intent(inout) :: cs !< Variable mixing coefficients
1379 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1380 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1381 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1382 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1383 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer vertical extents [Z ~> m]
1384 integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
1385 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence
1386 !! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1387 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence
1388 !! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1389 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1390 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1391 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
1392 !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1393 real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
1394 !! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1395 ! Local variables
1396 real, dimension(SZI_(G),SZJB_(G)) :: &
1397 dslopey_dz, & ! z-derivative of y-slope at v-points [L-1 ~> m-1]
1398 h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
1399 beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
1400 grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
1401 grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]
1402
1403 real, dimension(SZIB_(G),SZJ_(G)) :: &
1404 dslopex_dz, & ! z-derivative of x-slope at u-points [L-1 ~> m-1]
1405 h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
1406 beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
1407 grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
1408 grad_div_mag_u ! Magnitude of divergence gradient at u-points [T-1 L-1 ~> s-1 m-1]
1409 real :: h_at_slope_above ! The thickness above [H ~> m or kg m-2]
1410 real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2]
1411 real :: ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
1412 real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1]
1413 real :: z_to_h ! A local copy of depth to thickness conversion factors or the inverse of the
1414 ! mass-weighted average specific volumes around an interface [H Z-1 ~> nondim or kg m-3]
1415 real :: inv_pi3 ! The inverse of pi cubed [nondim]
1416 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq, nz
1417
1418 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1419 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1420 nz = gv%ke
1421
1422 inv_pi3 = 1.0 / ((4.0*atan(1.0))**3)
1423 z_to_h = gv%Z_to_H ! This will be replaced with a varying value in non-Boussinesq mode.
1424
1425 if ((k > 1) .and. (k < nz)) then
1426
1427 do j=js-2,je+2 ; do i=is-2,ie+1
1428 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) ) / &
1429 ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
1430 + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + gv%H_subroundoff**3 )
1431 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) ) / &
1432 ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
1433 + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + gv%H_subroundoff**3 )
1434 ih = 1./ ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff )
1435 if (.not.gv%Boussinesq) &
1436 z_to_h = ( (h(i,j,k-1) + h(i+1,j,k-1)) + (h(i,j,k) + h(i+1,j,k)) ) / &
1437 ( (dz(i,j,k-1) + dz(i+1,j,k-1)) + (dz(i,j,k) + dz(i+1,j,k)) + gv%dZ_subroundoff)
1438 dslopex_dz(i,j) = 2. * ( slope_x(i,j,k) - slope_x(i,j,k+1) ) * (z_to_h * ih)
1439 h_at_u(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
1440 enddo ; enddo
1441
1442 do j=js-2,je+1 ; do i=is-2,ie+2
1443 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) ) / &
1444 ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
1445 + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + gv%H_subroundoff**3 )
1446 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) ) / &
1447 ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
1448 + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + gv%H_subroundoff**3 )
1449 ih = 1./ ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff )
1450 if (.not.gv%Boussinesq) &
1451 z_to_h = ( (h(i,j,k-1) + h(i,j+1,k-1)) + (h(i,j,k) + h(i,j+1,k)) ) / &
1452 ( (dz(i,j,k-1) + dz(i,j+1,k-1)) + (dz(i,j,k) + dz(i,j+1,k)) + gv%dZ_subroundoff)
1453 dslopey_dz(i,j) = 2. * ( slope_y(i,j,k) - slope_y(i,j,k+1) ) * (z_to_h * ih)
1454 h_at_v(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
1455 enddo ; enddo
1456
1457 do j=js-2,je+1 ; do i=is-1,ie+1
1458 f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
1459 vort_xy_dx(i,j) = vort_xy_dx(i,j) - f * &
1460 ( ( (h_at_u(i,j) * dslopex_dz(i,j)) + (h_at_u(i-1,j+1) * dslopex_dz(i-1,j+1)) ) &
1461 + ( (h_at_u(i-1,j) * dslopex_dz(i-1,j)) + (h_at_u(i,j+1) * dslopex_dz(i,j+1)) ) ) / &
1462 ( ( 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)
1463 enddo ; enddo
1464
1465 do j=js-1,je+1 ; do i=is-2,ie+1
1466 f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
1467 vort_xy_dy(i,j) = vort_xy_dy(i,j) - f * &
1468 ( ( (h_at_v(i,j) * dslopey_dz(i,j)) + (h_at_v(i+1,j-1) * dslopey_dz(i+1,j-1)) ) &
1469 + ( (h_at_v(i,j-1) * dslopey_dz(i,j-1)) + (h_at_v(i+1,j) * dslopey_dz(i+1,j)) ) ) / &
1470 ( ( 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)
1471 enddo ; enddo
1472 endif ! k > 1
1473
1474 if (cs%use_QG_Leith_GM) then
1475
1476 do j=js,je ; do i=is-1,ieq
1477 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-1)) &
1478 + (vort_xy_dx(i+1,j) + vort_xy_dx(i,j-1))))**2)
1479 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-1)) &
1480 + (div_xx_dy(i+1,j) + div_xx_dy(i,j-1))))**2)
1481 if (cs%use_beta_in_QG_Leith) then
1482 beta_u(i,j) = sqrt((0.5*(g%dF_dx(i,j)+g%dF_dx(i+1,j))**2) + &
1483 (0.5*(g%dF_dy(i,j)+g%dF_dy(i+1,j))**2))
1484 cs%KH_u_QG(i,j,k) = min(grad_vort_mag_u(i,j) + grad_div_mag_u(i,j), 3.0*beta_u(i,j)) * &
1485 cs%Laplac3_const_u(i,j) * inv_pi3
1486 else
1487 cs%KH_u_QG(i,j,k) = (grad_vort_mag_u(i,j) + grad_div_mag_u(i,j)) * &
1488 cs%Laplac3_const_u(i,j) * inv_pi3
1489 endif
1490 enddo ; enddo
1491
1492 do j=js-1,jeq ; do i=is,ie
1493 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+1)) &
1494 + (vort_xy_dy(i,j+1) + vort_xy_dy(i-1,j))))**2)
1495 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+1)) &
1496 + (div_xx_dx(i,j+1) + div_xx_dx(i-1,j))))**2)
1497 if (cs%use_beta_in_QG_Leith) then
1498 beta_v(i,j) = sqrt((0.5*(g%dF_dx(i,j)+g%dF_dx(i,j+1))**2) + &
1499 (0.5*(g%dF_dy(i,j)+g%dF_dy(i,j+1))**2))
1500 cs%KH_v_QG(i,j,k) = min(grad_vort_mag_v(i,j) + grad_div_mag_v(i,j), 3.0*beta_v(i,j)) * &
1501 cs%Laplac3_const_v(i,j) * inv_pi3
1502 else
1503 cs%KH_v_QG(i,j,k) = (grad_vort_mag_v(i,j) + grad_div_mag_v(i,j)) * &
1504 cs%Laplac3_const_v(i,j) * inv_pi3
1505 endif
1506 enddo ; enddo
1507 ! post diagnostics
1508
1509 if (k==nz) then
1510 if (cs%id_KH_v_QG > 0) call post_data(cs%id_KH_v_QG, cs%KH_v_QG, cs%diag)
1511 if (cs%id_KH_u_QG > 0) call post_data(cs%id_KH_u_QG, cs%KH_u_QG, cs%diag)
1512 endif
1513 endif
1514
1515end subroutine calc_qg_leith_viscosity
1516
1517!> Initializes the variables mixing coefficients container
1518subroutine varmix_init(Time, G, GV, US, param_file, diag, CS)
1519 type(time_type), intent(in) :: time !< Current model time
1520 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1521 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1522 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1523 type(param_file_type), intent(in) :: param_file !< Parameter file handles
1524 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
1525 type(varmix_cs), intent(inout) :: cs !< Variable mixing coefficients
1526
1527 ! Local variables
1528 real :: khtr_slope_cff ! The nondimensional coefficient in the Visbeck formula
1529 ! for the epipycnal tracer diffusivity [nondim]
1530 real :: khth_slope_cff ! The nondimensional coefficient in the Visbeck formula
1531 ! for the interface depth diffusivity [nondim]
1532 real :: oneortwo ! A variable that may be 1 or 2, depending on which form
1533 ! of the equatorial deformation radius us used [nondim]
1534 real :: n2_filter_depth ! A depth below which stratification is treated as monotonic when
1535 ! calculating the first-mode wave speed [H ~> m or kg m-2]
1536 real :: khtr_passivity_coeff ! Coefficient setting the ratio between along-isopycnal tracer
1537 ! mixing and interface height mixing [nondim]
1538 real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
1539 ! default value is roughly (pi / (the age of the universe)).
1540 logical :: gill_equatorial_ld, use_fgnv_streamfn, use_meke, in_use
1541 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1542 integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1543 ! for remapping. Values below 20190101 recover the remapping
1544 ! answers from 2018, while higher values use more robust
1545 ! forms of the same remapping expressions.
1546 real :: mle_front_length ! The frontal-length scale used to calculate the upscaling of
1547 ! buoyancy gradients in boundary layer parameterizations [L ~> m]
1548 real :: leith_lap_const ! The non-dimensional coefficient in the Leith viscosity [nondim]
1549 real :: grid_sp_u2, grid_sp_v2 ! Intermediate quantities for Leith metrics [L2 ~> m2]
1550 real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
1551 real :: wave_speed_min ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
1552 real :: wave_speed_tol ! The fractional tolerance for finding the wave speeds [nondim]
1553 logical :: resoln_scaled_meke_visc ! If true, the viscosity contribution from MEKE is
1554 ! scaled by the resolution function.
1555 logical :: better_speed_est ! If true, use a more robust estimate of the first
1556 ! mode wave speed as the starting point for iterations.
1557 real :: stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
1558 ! temperature variance [nondim]
1559 logical :: use_sqg ! This is true if the SQG structure will be used for any parameterizations.
1560 logical :: om4_remap_via_sub_cells ! Use the OM4-era remap_via_sub_cells for calculating the EBT structure
1561 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1562 ! recreate the bugs, or if false bugs are only used if actively selected.
1563 logical :: mixing_coefs_obc_bug ! If false, use only interior data for thickness weighting in
1564 ! lateral mixing coefficient calculations and to calculate stratification
1565 ! and other fields at open boundary condition faces.
1566 ! This include declares and sets the variable "version".
1567# include "version_variable.h"
1568 character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
1569 integer :: number_of_obc_segments
1570 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
1571 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1572 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1573 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1574 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1575 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1576
1577 cs%initialized = .true.
1578 in_use = .false. ! Set to true to avoid deallocating
1579 cs%diag => diag ! Diagnostics pointer
1580 cs%calculate_cg1 = .false.
1581 cs%calculate_Rd_dx = .false.
1582 cs%calculate_res_fns = .false.
1583 cs%use_simpler_Eady_growth_rate = .false.
1584 cs%full_depth_Eady_growth_rate = .false.
1585 cs%calculate_depth_fns = .false.
1586 ! Read all relevant parameters and write them to the model log.
1587 call log_version(param_file, mdl, version, "")
1588 call get_param(param_file, mdl, "USE_VARIABLE_MIXING", cs%use_variable_mixing,&
1589 "If true, the variable mixing code will be called. This "//&
1590 "allows diagnostics to be created even if the scheme is "//&
1591 "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
1592 "this is set to true regardless of what is in the "//&
1593 "parameter file.", default=.false.)
1594 call get_param(param_file, mdl, "USE_VISBECK", cs%use_Visbeck,&
1595 "If true, use the Visbeck et al. (1997) formulation for \n"//&
1596 "thickness diffusivity.", default=.false.)
1597 call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
1598 "If true, the Laplacian lateral viscosity is scaled away "//&
1599 "when the first baroclinic deformation radius is well "//&
1600 "resolved.", default=.false.)
1601 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", cs%Depth_scaled_KhTh, &
1602 "If true, KHTH is scaled away when the depth is shallower "//&
1603 "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//&
1604 "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
1605 "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
1606 default=.false.)
1607 call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
1608 "If true, the interface depth diffusivity is scaled away "//&
1609 "when the first baroclinic deformation radius is well "//&
1610 "resolved.", default=.false.)
1611 call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
1612 "If true, the epipycnal tracer diffusivity is scaled "//&
1613 "away when the first baroclinic deformation radius is "//&
1614 "well resolved.", default=.false.)
1615 call get_param(param_file, mdl, "USE_MEKE", use_meke, &
1616 default=.false., do_not_log=.true.)
1617 call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", resoln_scaled_meke_visc, &
1618 "If true, the viscosity contribution from MEKE is scaled by "//&
1619 "the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
1620 if (.not.use_meke) resoln_scaled_meke_visc = .false.
1621 call get_param(param_file, mdl, "RESOLN_USE_EBT", cs%Resoln_use_ebt, &
1622 "If true, uses the equivalent barotropic wave speed instead "//&
1623 "of first baroclinic wave for calculating the resolution function.",&
1624 default=.false.)
1625 call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", cs%BS_EBT_power, &
1626 "Power to raise EBT vertical structure to when backscatter "// &
1627 "has vertical structure.", units="nondim", default=0.0)
1628 call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", cs%BS_use_sqg_struct, &
1629 "If true, the SQG vertical structure is used for backscatter "//&
1630 "on the condition that BS_EBT_power=0", &
1631 default=.false.)
1632 call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
1633 "If true, uses the equivalent barotropic structure "//&
1634 "as the vertical structure of thickness diffusivity.",&
1635 default=.false.)
1636 call get_param(param_file, mdl, "KHTH_USE_SQG_STRUCT", cs%khth_use_sqg_struct, &
1637 "If true, uses the surface quasigeostrophic structure "//&
1638 "as the vertical structure of thickness diffusivity.",&
1639 default=.false.)
1640 call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", cs%khtr_use_ebt_struct, &
1641 "If true, uses the equivalent barotropic structure "//&
1642 "as the vertical structure of tracer diffusivity.",&
1643 default=.false.)
1644 call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", cs%khtr_use_sqg_struct, &
1645 "If true, uses the surface quasigeostrophic structure "//&
1646 "as the vertical structure of tracer diffusivity.",&
1647 default=.false.)
1648 call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", cs%kdgl90_use_ebt_struct, &
1649 "If true, uses the equivalent barotropic structure "//&
1650 "as the vertical structure of diffusivity in the GL90 scheme.",&
1651 default=.false.)
1652 call get_param(param_file, mdl, "KD_GL90_USE_SQG_STRUCT", cs%kdgl90_use_sqg_struct, &
1653 "If true, uses the equivalent barotropic structure "//&
1654 "as the vertical structure of diffusivity in the GL90 scheme.",&
1655 default=.false.)
1656 call get_param(param_file, mdl, "KHTH_SLOPE_CFF", khth_slope_cff, &
1657 "The nondimensional coefficient in the Visbeck formula "//&
1658 "for the interface depth diffusivity", units="nondim", default=0.0)
1659 call get_param(param_file, mdl, "KHTR_SLOPE_CFF", khtr_slope_cff, &
1660 "The nondimensional coefficient in the Visbeck formula "//&
1661 "for the epipycnal tracer diffusivity", units="nondim", default=0.0)
1662 call get_param(param_file, mdl, "USE_STORED_SLOPES", cs%use_stored_slopes,&
1663 "If true, the isopycnal slopes are calculated once and "//&
1664 "stored for re-use. This uses more memory but avoids calling "//&
1665 "the equation of state more times than should be necessary.", &
1666 default=.false.)
1667 call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1668 "A miniscule frequency that is used to avoid division by 0. The default "//&
1669 "value is roughly (pi / (the age of the universe)).", &
1670 default=1.0e-17, units="s-1", scale=us%T_to_s)
1671 call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
1672 default=.false., do_not_log=.true.)
1673 cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn .or. cs%khth_use_ebt_struct &
1674 .or. cs%kdgl90_use_ebt_struct .or. cs%BS_EBT_power>0.
1675 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
1676 ! Indicate whether to calculate the Eady growth rate
1677 cs%calculate_Eady_growth_rate = use_meke .or. (khtr_slope_cff>0.) .or. (khth_slope_cff>0.)
1678 call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
1679 units="nondim", default=0., do_not_log=.true.)
1680 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
1681 call get_param(param_file, mdl, "MLE_FRONT_LENGTH", mle_front_length, &
1682 units="m", default=0.0, scale=us%m_to_L, do_not_log=.true.)
1683 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
1684
1685 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
1686
1687 call get_param(param_file, mdl, "USE_STANLEY_ISO", cs%use_stanley_iso, &
1688 "If true, turn on Stanley SGS T variance parameterization "// &
1689 "in isopycnal slope code.", default=.false.)
1690 if (cs%use_stanley_iso) then
1691 call get_param(param_file, mdl, "STANLEY_COEFF", stanley_coeff, &
1692 "Coefficient correlating the temperature gradient and SGS T variance.", &
1693 units="nondim", default=-1.0, do_not_log=.true.)
1694 if (stanley_coeff < 0.0) call mom_error(fatal, &
1695 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.")
1696 endif
1697 call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_obc_segments, &
1698 default=0, do_not_log=.true.)
1699 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1700 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1701 call get_param(param_file, mdl, "MIXING_COEFS_OBC_BUG", mixing_coefs_obc_bug, &
1702 "If false, use only interior data for thickness weighting in lateral mixing "//&
1703 "coefficient calculations and to calculate stratification and other fields at "//&
1704 "open boundary condition faces.", &
1705 default=enable_bugs, do_not_log=(number_of_obc_segments<=0))
1706 cs%OBC_friendly = .not. mixing_coefs_obc_bug
1707 call get_param(param_file, mdl, "RESOLN_FUNCTION_OBC_BUG", cs%res_fn_OBC_bug, &
1708 "If false, use only interior data for calculating the resolution functions at "//&
1709 "open boundary condition faces and vertices.", &
1710 default=enable_bugs, do_not_log=(number_of_obc_segments<=0))
1711
1712 if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct .or. cs%kdgl90_use_ebt_struct &
1713 .or. cs%BS_EBT_power>0. .or. cs%khtr_use_ebt_struct) then
1714 in_use = .true.
1715 call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
1716 "The depth below which N2 is monotonized to avoid stratification "//&
1717 "artifacts from altering the equivalent barotropic mode structure. "//&
1718 "This monotonzization is disabled if this parameter is negative.", &
1719 units="m", default=-1.0, scale=gv%m_to_H)
1720 allocate(cs%ebt_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1721 endif
1722
1723 use_sqg = cs%BS_use_sqg_struct .or. cs%khth_use_sqg_struct .or. cs%khtr_use_sqg_struct .or. &
1724 cs%kdgl90_use_sqg_struct
1725 call get_param(param_file, mdl, "SQG_EXPO", cs%sqg_expo, &
1726 "Nondimensional exponent coeffecient of the SQG mode that is used for the "//&
1727 "vertical struture of diffusivities.", &
1728 units="nondim", default=1.0, do_not_log=.not.use_sqg)
1729 call get_param(param_file, mdl, "INTERPOLATED_SQG_STRUCTURE", cs%interpolated_sqg_struct, &
1730 "If true, interpolate properties to velocity points and then interpolate the "//&
1731 "buoyancy frequencies and layer thicknesses back to tracer points when "//&
1732 "calculating the SQG vertical structure.", &
1733 default=.true., do_not_log=.not.use_sqg)
1734 !### Consider changing the default for INTERPOLATED_SQG_STRUCTURE to false.
1735
1736 if ((cs%BS_EBT_power>0.) .and. cs%BS_use_sqg_struct) call mom_error(fatal, &
1737 "calc_resoln_function: BS_EBT_POWER>0. and BS_USE_SQG=True cannot be set together")
1738
1739 if (cs%khth_use_ebt_struct .and. cs%khth_use_sqg_struct) call mom_error(fatal, &
1740 "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT and KHTH_USE_SQG_STRUCT can be true")
1741
1742 if (cs%khtr_use_ebt_struct .and. cs%khtr_use_sqg_struct) call mom_error(fatal, &
1743 "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT and KHTR_USE_SQG_STRUCT can be true")
1744
1745 if (cs%kdgl90_use_ebt_struct .and. cs%kdgl90_use_sqg_struct) call mom_error(fatal, &
1746 "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT and KD_GL90_USE_SQG_STRUCT can be true")
1747
1748 if (cs%BS_EBT_power>0. .or. cs%BS_use_sqg_struct) then
1749 allocate(cs%BS_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1750 endif
1751
1752 if (cs%khth_use_ebt_struct .or. cs%khth_use_sqg_struct) then
1753 allocate(cs%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1754 endif
1755
1756 if (cs%khtr_use_ebt_struct .or. cs%khtr_use_sqg_struct) then
1757 allocate(cs%khtr_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1758 endif
1759
1760 if (cs%kdgl90_use_ebt_struct .or. cs%kdgl90_use_sqg_struct) then
1761 allocate(cs%kdgl90_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1762 endif
1763
1764 if (cs%use_stored_slopes) then
1765 if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1766 call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
1767 "If non-zero, is an upper bound on slopes used in the "//&
1768 "Visbeck formula for diffusivity. This does not affect the "//&
1769 "isopycnal slope calculation used within thickness diffusion.", &
1770 units="nondim", default=0.0, scale=us%L_to_Z)
1771 else
1772 cs%Visbeck_S_max = 0.
1773 endif
1774 endif
1775
1776 if (cs%use_stored_slopes .or. (cs%interpolated_sqg_struct .and. (cs%sqg_expo>0.0))) then
1777 ! CS%calculate_Eady_growth_rate=.true.
1778 in_use = .true.
1779 allocate(cs%slope_x(isdb:iedb,jsd:jed,gv%ke+1), source=0.0)
1780 allocate(cs%slope_y(isd:ied,jsdb:jedb,gv%ke+1), source=0.0)
1781 call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
1782 "A diapycnal diffusivity that is used to interpolate "//&
1783 "more sensible values of T & S into thin layers.", &
1784 units="m2 s-1", default=1.0e-6, scale=gv%m2_s_to_HZ_T)
1785 endif
1786
1787 if (cs%calculate_Eady_growth_rate) then
1788 in_use = .true.
1789 allocate(cs%SN_u(isdb:iedb,jsd:jed), source=0.0)
1790 allocate(cs%SN_v(isd:ied,jsdb:jedb), source=0.0)
1791 cs%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, time, &
1792 'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=us%s_to_T)
1793 cs%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, time, &
1794 'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=us%s_to_T)
1795 call get_param(param_file, mdl, "USE_SIMPLER_EADY_GROWTH_RATE", cs%use_simpler_Eady_growth_rate, &
1796 "If true, use a simpler method to calculate the Eady growth rate "//&
1797 "that avoids division by layer thickness. Recommended.", default=.false.)
1798 if (cs%use_simpler_Eady_growth_rate) then
1799 if (.not. cs%use_stored_slopes) call mom_error(fatal, &
1800 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1801 "When USE_SIMPLER_EADY_GROWTH_RATE=True, USE_STORED_SLOPES must also be True.")
1802 call get_param(param_file, mdl, "EADY_GROWTH_RATE_D_SCALE", cs%Eady_GR_D_scale, &
1803 "The depth from surface over which to average SN when calculating "//&
1804 "a 2D Eady growth rate. Zero mean use full depth.", &
1805 units="m", default=0., scale=us%m_to_Z)
1806 call get_param(param_file, mdl, "EADY_GROWTH_RATE_CROPPING_DISTANCE", cs%cropping_distance, &
1807 "Distance from surface or bottom to filter out outcropped or "//&
1808 "incropped interfaces for the Eady growth rate calc. "//&
1809 "Negative values disables cropping.", units="m", default=0., scale=us%m_to_Z)
1810 else
1811 call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
1812 "The layer number at which to start vertical integration "//&
1813 "of S*N for purposes of finding the Eady growth rate.", &
1814 units="nondim", default=2)
1815 call get_param(param_file, mdl, "MIN_DZ_FOR_SLOPE_N2", cs%h_min_N2, &
1816 "The minimum vertical distance to use in the denominator of the "//&
1817 "bouyancy frequency used in the slope calculation.", &
1818 units="m", default=1.0, scale=gv%m_to_H, do_not_log=cs%use_stored_slopes)
1819
1820 call get_param(param_file, mdl, "FULL_DEPTH_EADY_GROWTH_RATE", cs%full_depth_Eady_growth_rate, &
1821 "If true, calculate the Eady growth rate based on average slope times "//&
1822 "stratification that includes contributions from sea-level changes "//&
1823 "in its denominator, rather than just the nominal depth of the bathymetry. "//&
1824 "This only applies when using the model interface heights as a proxy for "//&
1825 "isopycnal slopes.", default=.not.(gv%Boussinesq.or.gv%semi_Boussinesq), &
1826 do_not_log=cs%use_stored_slopes)
1827 endif
1828 endif
1829
1830 if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1831 in_use = .true.
1832 call get_param(param_file, mdl, "VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1833 "The fixed length scale in the Visbeck formula, or if negative a nondimensional "//&
1834 "scaling factor relating this length scale squared to the cell areas.", &
1835 units="m or nondim", default=0.0, scale=us%m_to_L)
1836 allocate(cs%L2u(isdb:iedb,jsd:jed), source=0.0)
1837 allocate(cs%L2v(isd:ied,jsdb:jedb), source=0.0)
1838 if (cs%Visbeck_L_scale<0) then
1839 ! Undo the rescaling of CS%Visbeck_L_scale.
1840 do j=js,je ; do i=is-1,ieq
1841 cs%L2u(i,j) = (us%L_to_m*cs%Visbeck_L_scale)**2 * g%areaCu(i,j)
1842 enddo ; enddo
1843 do j=js-1,jeq ; do i=is,ie
1844 cs%L2v(i,j) = (us%L_to_m*cs%Visbeck_L_scale)**2 * g%areaCv(i,j)
1845 enddo ; enddo
1846 else
1847 cs%L2u(:,:) = cs%Visbeck_L_scale**2
1848 cs%L2v(:,:) = cs%Visbeck_L_scale**2
1849 endif
1850
1851 cs%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, time, &
1852 'Length scale squared for mixing coefficient, at u-points', &
1853 'm2', conversion=us%L_to_m**2)
1854 cs%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, time, &
1855 'Length scale squared for mixing coefficient, at v-points', &
1856 'm2', conversion=us%L_to_m**2)
1857 endif
1858
1859 cs%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, time, &
1860 'Vertical structure of SQG mode', 'nondim')
1861 if (cs%BS_use_sqg_struct .or. cs%khth_use_sqg_struct .or. cs%khtr_use_sqg_struct &
1862 .or. cs%kdgl90_use_sqg_struct .or. cs%id_sqg_struct>0) then
1863 allocate(cs%sqg_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1864 endif
1865
1866 if (cs%BS_EBT_power>0. .or. cs%BS_use_sqg_struct) then
1867 cs%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, time, &
1868 'Vertical structure of backscatter', 'nondim')
1869 endif
1870 if (cs%khth_use_ebt_struct .or. cs%khth_use_sqg_struct) then
1871 cs%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, time, &
1872 'Vertical structure of thickness diffusivity', 'nondim')
1873 endif
1874 if (cs%khtr_use_ebt_struct .or. cs%khtr_use_sqg_struct) then
1875 cs%id_khtr_struct = register_diag_field('ocean_model', 'khtr_struct', diag%axesTl, time, &
1876 'Vertical structure of tracer diffusivity', 'nondim')
1877 endif
1878 if (cs%kdgl90_use_ebt_struct .or. cs%kdgl90_use_sqg_struct) then
1879 cs%id_kdgl90_struct = register_diag_field('ocean_model', 'kdgl90_struct', diag%axesTl, time, &
1880 'Vertical structure of GL90 diffusivity', 'nondim')
1881 endif
1882
1883 if ((cs%calculate_Eady_growth_rate .and. cs%use_stored_slopes) ) then
1884 cs%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, time, &
1885 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
1886 's-2', conversion=(us%L_to_Z*us%s_to_T)**2)
1887 cs%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, time, &
1888 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', &
1889 's-2', conversion=(us%L_to_Z*us%s_to_T)**2)
1890 endif
1891 if (cs%use_simpler_Eady_growth_rate) then
1892 cs%id_dzu = register_diag_field('ocean_model', 'dzu_Visbeck', diag%axesCui, time, &
1893 'dz at u-points, used in calculating Eady growth rate in Visbeck et al..', &
1894 'm', conversion=us%Z_to_m)
1895 cs%id_dzv = register_diag_field('ocean_model', 'dzv_Visbeck', diag%axesCvi, time, &
1896 'dz at v-points, used in calculating Eady growth rate in Visbeck et al..', &
1897 'm', conversion=us%Z_to_m)
1898 cs%id_dzSxN = register_diag_field('ocean_model', 'dzSxN', diag%axesCui, time, &
1899 'dz * |slope_x| * N, used in calculating Eady growth rate in '//&
1900 'Visbeck et al..', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1901 cs%id_dzSyN = register_diag_field('ocean_model', 'dzSyN', diag%axesCvi, time, &
1902 'dz * |slope_y| * N, used in calculating Eady growth rate in '//&
1903 'Visbeck et al..', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1904 endif
1905 if (cs%use_stored_slopes) then
1906 cs%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, time, &
1907 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', &
1908 'nondim', conversion=us%Z_to_L**2)
1909 cs%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, time, &
1910 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', &
1911 'nondim', conversion=us%Z_to_L**2)
1912 endif
1913
1914 oneortwo = 1.0
1915 cs%Resoln_scaling_used = cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. &
1916 cs%Resoln_scaled_KhTr .or. resoln_scaled_meke_visc
1917 if (cs%Resoln_scaling_used) then
1918 cs%calculate_Rd_dx = .true.
1919 cs%calculate_res_fns = .true.
1920 allocate(cs%Res_fn_h(isd:ied,jsd:jed), source=0.0)
1921 allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb), source=0.0)
1922 allocate(cs%Res_fn_u(isdb:iedb,jsd:jed), source=0.0)
1923 allocate(cs%Res_fn_v(isd:ied,jsdb:jedb), source=0.0)
1924 allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb), source=0.0)
1925 allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed), source=0.0)
1926 allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb), source=0.0)
1927 allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb), source=0.0)
1928 allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed), source=0.0)
1929 allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb), source=0.0)
1930
1931 cs%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, time, &
1932 'Resolution function for scaling diffusivities', 'nondim')
1933
1934 call get_param(param_file, mdl, "KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1935 "A coefficient that determines how KhTh is scaled away if "//&
1936 "RESOLN_SCALED_... is true, as "//&
1937 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1938 units="nondim", default=1.0)
1939 call get_param(param_file, mdl, "KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1940 "The power of dx/Ld in the Kh resolution function. Any "//&
1941 "positive integer may be used, although even integers "//&
1942 "are more efficient to calculate. Setting this greater "//&
1943 "than 100 results in a step-function being used.", &
1944 default=2)
1945 call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1946 "A coefficient that determines how Kh is scaled away if "//&
1947 "RESOLN_SCALED_... is true, as "//&
1948 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1949 "This function affects lateral viscosity, Kh, and not KhTh.", &
1950 units="nondim", default=cs%Res_coef_khth)
1951 call get_param(param_file, mdl, "VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1952 "The power of dx/Ld in the Kh resolution function. Any "//&
1953 "positive integer may be used, although even integers "//&
1954 "are more efficient to calculate. Setting this greater "//&
1955 "than 100 results in a step-function being used. "//&
1956 "This function affects lateral viscosity, Kh, and not KhTh.", &
1957 default=cs%Res_fn_power_khth)
1958 call get_param(param_file, mdl, "INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1959 "If true, interpolate the resolution function to the "//&
1960 "velocity points from the thickness points; otherwise "//&
1961 "interpolate the wave speed and calculate the resolution "//&
1962 "function independently at each point.", default=.false.)
1963 if (cs%interpolate_Res_fn) then
1964 if (cs%Res_coef_visc /= cs%Res_coef_khth) call mom_error(fatal, &
1965 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1966 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1967 if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth) call mom_error(fatal, &
1968 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1969 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1970 endif
1971 call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1972 "If true, uses Gill's definition of the baroclinic "//&
1973 "equatorial deformation radius, otherwise, if false, use "//&
1974 "Pedlosky's definition. These definitions differ by a factor "//&
1975 "of 2 in front of the beta term in the denominator. Gill's "//&
1976 "is the more appropriate definition.", default=.true.)
1977 if (gill_equatorial_ld) then
1978 oneortwo = 2.0
1979 endif
1980
1981 do j=js-1,jeq ; do i=is-1,ieq
1982 cs%f2_dx2_q(i,j) = ((g%dxBu(i,j)**2) + (g%dyBu(i,j)**2)) * &
1983 max(g%Coriolis2Bu(i,j), absurdly_small_freq**2)
1984 cs%beta_dx2_q(i,j) = oneortwo * ((g%dxBu(i,j)**2) + (g%dyBu(i,j)**2)) * (sqrt(0.5 * &
1985 ( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) + &
1986 (((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2)) + &
1987 ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
1988 (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2)) ) ))
1989 enddo ; enddo
1990
1991 do j=js,je ; do i=is-1,ieq
1992 cs%f2_dx2_u(i,j) = ((g%dxCu(i,j)**2) + (g%dyCu(i,j)**2)) * &
1993 max(0.5* (g%Coriolis2Bu(i,j)+g%Coriolis2Bu(i,j-1)), absurdly_small_freq**2)
1994 cs%beta_dx2_u(i,j) = oneortwo * ((g%dxCu(i,j)**2) + (g%dyCu(i,j)**2)) * (sqrt( &
1995 ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1996 0.25*( ((((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1997 (((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2)) + &
1998 ((((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2) + &
1999 (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2)) ) ))
2000 enddo ; enddo
2001
2002 do j=js-1,jeq ; do i=is,ie
2003 cs%f2_dx2_v(i,j) = ((g%dxCv(i,j)**2) + (g%dyCv(i,j)**2)) * &
2004 max(0.5*(g%Coriolis2Bu(i,j)+g%Coriolis2Bu(i-1,j)), absurdly_small_freq**2)
2005 cs%beta_dx2_v(i,j) = oneortwo * ((g%dxCv(i,j)**2) + (g%dyCv(i,j)**2)) * (sqrt( &
2006 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
2007 0.25*( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
2008 (((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2)) + &
2009 ((((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) + &
2010 (((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2)) ) ))
2011 enddo ; enddo
2012
2013 endif
2014
2015 if (cs%Depth_scaled_KhTh) then
2016 cs%calculate_depth_fns = .true.
2017 allocate(cs%Depth_fn_u(isdb:iedb,jsd:jed), source=0.0)
2018 allocate(cs%Depth_fn_v(isd:ied,jsdb:jedb), source=0.0)
2019 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", cs%depth_scaled_khth_h0, &
2020 "The depth above which KHTH is scaled away.", &
2021 units="m", scale=us%m_to_Z, default=1000.)
2022 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", cs%depth_scaled_khth_exp, &
2023 "The exponent used in the depth dependent scaling function for KHTH.", &
2024 units="nondim", default=3.0)
2025 endif
2026
2027 ! Resolution %Rd_dx_h
2028 cs%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, time, &
2029 'Ratio between deformation radius and grid spacing', 'm m-1')
2030 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
2031
2032 if (cs%calculate_Rd_dx) then
2033 cs%calculate_cg1 = .true. ! We will need %cg1
2034 allocate(cs%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
2035 allocate(cs%beta_dx2_h(isd:ied,jsd:jed), source=0.0)
2036 allocate(cs%f2_dx2_h(isd:ied,jsd:jed), source=0.0)
2037 do j=js-1,je+1 ; do i=is-1,ie+1
2038 cs%f2_dx2_h(i,j) = ((g%dxT(i,j)**2) + (g%dyT(i,j)**2)) * &
2039 max(0.25 * ((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
2040 (g%Coriolis2Bu(i-1,j) + g%Coriolis2Bu(i,j-1))), &
2041 absurdly_small_freq**2)
2042 cs%beta_dx2_h(i,j) = oneortwo * ((g%dxT(i,j)**2) + (g%dyT(i,j)**2)) * (sqrt(0.5 * &
2043 ( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) + &
2044 (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2)) + &
2045 ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
2046 (((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2)) ) ))
2047 enddo ; enddo
2048 endif
2049
2050 if (cs%calculate_cg1) then
2051 in_use = .true.
2052 allocate(cs%cg1(isd:ied,jsd:jed), source=0.0)
2053 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2054 "This sets the default value for the various _ANSWER_DATE parameters.", &
2055 default=99991231)
2056 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2057 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
2058 "Values below 20190101 result in the use of older, less accurate expressions "//&
2059 "that were in use at the end of 2018. Higher values result in the use of more "//&
2060 "robust and accurate forms of mathematically equivalent expressions.", &
2061 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
2062 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2063
2064 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
2065 "The fractional tolerance for finding the wave speeds.", &
2066 units="nondim", default=0.001)
2067 !### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
2068 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
2069 "A floor in the first mode speed below which 0 used instead.", &
2070 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
2071 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
2072 "If true, use a more robust estimate of the first mode wave speed as the "//&
2073 "starting point for iterations.", default=.true.)
2074 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2075 do_not_log=.true., default=.true.)
2076 call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2077 "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
2078 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
2079 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2080 call wave_speed_init(cs%wave_speed, gv, use_ebt_mode=cs%Resoln_use_ebt, &
2081 mono_n2_depth=n2_filter_depth, remap_answer_date=remap_answer_date, &
2082 better_speed_est=better_speed_est, min_speed=wave_speed_min, &
2083 om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol)
2084 endif
2085
2086 ! Leith parameters
2087 call get_param(param_file, mdl, "USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
2088 "If true, use the QG Leith viscosity as the GM coefficient.", &
2089 default=.false.)
2090
2091 if (cs%Use_QG_Leith_GM) then
2092 call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
2093 "The nondimensional Laplacian Leith constant, \n"//&
2094 "often set to 1.0", units="nondim", default=0.0)
2095
2096 call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
2097 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
2098 default=.true.)
2099
2100 allocate(cs%Laplac3_const_u(isdb:iedb,jsd:jed), source=0.0)
2101 allocate(cs%Laplac3_const_v(isd:ied,jsdb:jedb), source=0.0)
2102 allocate(cs%KH_u_QG(isdb:iedb,jsd:jed,gv%ke), source=0.0)
2103 allocate(cs%KH_v_QG(isd:ied,jsdb:jedb,gv%ke), source=0.0)
2104
2105 ! register diagnostics
2106 cs%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, time, &
2107 'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2108 cs%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, time, &
2109 'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2110
2111 do j=jsq,jeq+1 ; do i=is-1,ieq
2112 ! Static factors in the Leith schemes
2113 grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
2114 grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
2115 cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
2116 enddo ; enddo
2117 do j=js-1,jeq ; do i=isq,ieq+1
2118 ! Static factors in the Leith schemes
2119 grid_sp_v2 = g%dyCv(i,j)*g%dxCv(i,j)
2120 grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
2121 cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
2122 enddo ; enddo
2123
2124 if (.not. cs%use_stored_slopes) call mom_error(fatal, &
2125 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2126 "USE_STORED_SLOPES must be True when using QG Leith.")
2127 endif
2128
2129 ! Re-enable variable mixing if one of the schemes was enabled
2130 cs%use_variable_mixing = in_use .or. cs%use_variable_mixing
2131end subroutine varmix_init
2132
2133!> Destructor for VarMix control structure
2134subroutine varmix_end(CS)
2135 type(varmix_cs), intent(inout) :: cs
2136
2137 if (allocated(cs%ebt_struct)) deallocate(cs%ebt_struct)
2138 if (allocated(cs%sqg_struct)) deallocate(cs%sqg_struct)
2139 if (allocated(cs%BS_struct)) deallocate(cs%BS_struct)
2140 if (allocated(cs%khth_struct)) deallocate(cs%khth_struct)
2141 if (allocated(cs%khtr_struct)) deallocate(cs%khtr_struct)
2142 if (allocated(cs%kdgl90_struct)) deallocate(cs%kdgl90_struct)
2143
2144 if (allocated(cs%slope_x)) deallocate(cs%slope_x)
2145 if (allocated(cs%slope_y)) deallocate(cs%slope_y)
2146
2147 if (allocated(cs%SN_u)) deallocate(cs%SN_u)
2148 if (allocated(cs%SN_v)) deallocate(cs%SN_v)
2149
2150 if (allocated(cs%L2u)) deallocate(cs%L2u)
2151 if (allocated(cs%L2v)) deallocate(cs%L2v)
2152
2153 if (allocated(cs%Res_fn_h)) deallocate(cs%Res_fn_h)
2154 if (allocated(cs%Res_fn_q)) deallocate(cs%Res_fn_q)
2155 if (allocated(cs%Res_fn_u)) deallocate(cs%Res_fn_u)
2156 if (allocated(cs%Res_fn_v)) deallocate(cs%Res_fn_v)
2157 if (allocated(cs%beta_dx2_q)) deallocate(cs%beta_dx2_q)
2158 if (allocated(cs%beta_dx2_u)) deallocate(cs%beta_dx2_u)
2159 if (allocated(cs%beta_dx2_v)) deallocate(cs%beta_dx2_v)
2160 if (allocated(cs%f2_dx2_q)) deallocate(cs%f2_dx2_q)
2161 if (allocated(cs%f2_dx2_u)) deallocate(cs%f2_dx2_u)
2162 if (allocated(cs%f2_dx2_v)) deallocate(cs%f2_dx2_v)
2163
2164 if (allocated(cs%Depth_fn_u)) deallocate(cs%Depth_fn_u)
2165 if (allocated(cs%Depth_fn_v)) deallocate(cs%Depth_fn_v)
2166
2167 if (allocated(cs%Rd_dx_h)) deallocate(cs%Rd_dx_h)
2168 if (allocated(cs%beta_dx2_h)) deallocate(cs%beta_dx2_h)
2169 if (allocated(cs%f2_dx2_h)) deallocate(cs%f2_dx2_h)
2170
2171 if (allocated(cs%cg1)) deallocate(cs%cg1)
2172
2173 if (allocated(cs%Laplac3_const_u)) deallocate(cs%Laplac3_const_u)
2174 if (allocated(cs%Laplac3_const_v)) deallocate(cs%Laplac3_const_v)
2175 if (allocated(cs%KH_u_QG)) deallocate(cs%KH_u_QG)
2176 if (allocated(cs%KH_v_QG)) deallocate(cs%KH_v_QG)
2177
2178end subroutine varmix_end
2179
2180!> \namespace mom_lateral_mixing_coeffs
2181!!
2182!! This module provides a container for various factors used in prescribing diffusivities, that are
2183!! a function of the state (in particular the stratification and isoneutral slopes).
2184!!
2185!! \section section_Resolution_Function The resolution function
2186!!
2187!! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
2188!! The square of the resolution parameter is
2189!!
2190!! \f[
2191!! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
2192!! \f]
2193!!
2194!! where the grid spacing is calculated as
2195!!
2196!! \f[
2197!! \Delta^2 = \Delta x^2 + \Delta y^2 .
2198!! \f]
2199!!
2200!! \todo Check this reference to Bob on/off paper.
2201!! The resolution function used in scaling diffusivities (\cite hallberg2013) is
2202!!
2203!! \f[
2204!! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
2205!! \f]
2206!!
2207!! The resolution function can be applied independently to thickness diffusion \(module mom_thickness_diffuse\),
2208!! tracer diffusion \(mom_tracer_hordiff\) lateral viscosity \(mom_hor_visc\).
2209!!
2210!! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
2211!! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007
2212!!
2213!! | Symbol | Module parameter |
2214!! | ------ | --------------- |
2215!! | - | <code>USE_VARIABLE_MIXING</code> |
2216!! | - | <code>RESOLN_SCALED_KH</code> |
2217!! | - | <code>RESOLN_SCALED_KHTH</code> |
2218!! | - | <code>RESOLN_SCALED_KHTR</code> |
2219!! | \f$ \alpha \f$ | <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) |
2220!! | \f$ p \f$ | <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) |
2221!! | \f$ \alpha \f$ | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) |
2222!! | \f$ p \f$ | <code>VISC_RES_FN_POWER</code> (for lateral viscosity) |
2223!! | - | <code>GILL_EQUATORIAL_LD</code> |
2224!!
2225!!
2226!!
2227!! \section section_Vicbeck Visbeck diffusivity
2228!!
2229!! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
2230!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse but calculated in this module.
2231!!
2232!! \f[
2233!! \kappa_h = \alpha_s L_s^2 S N
2234!! \f]
2235!!
2236!! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
2237!!
2238!! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
2239!! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
2240!!
2241!! | Symbol | Module parameter |
2242!! | ------ | --------------- |
2243!! | - | <code>USE_VARIABLE_MIXING</code> |
2244!! | \f$ \alpha_s \f$ | <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)|
2245!! | \f$ \alpha_s \f$ | <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)|
2246!! | \f$ L_{s} \f$ | <code>VISBECK_L_SCALE</code> |
2247!! | \f$ S_{max} \f$ | <code>VISBECK_MAX_SLOPE</code> |
2248!!
2249!!
2250!! \section section_vertical_structure_khth Vertical structure function for KhTh
2251!!
2252!! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
2253!! velocity mode. The structure function is stored in the control structure for this module (varmix_cs) but is
2254!! calculated using subroutines in mom_wave_speed.
2255!!
2256!! | Symbol | Module parameter |
2257!! | ------ | --------------- |
2258!! | - | <code>KHTH_USE_EBT_STRUCT</code> |
2259