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 logical :: stoch_eos ! Can't use Stanley param here unless stoch_eos is true
1567 ! This include declares and sets the variable "version".
1568# include "version_variable.h"
1569 character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
1570 integer :: number_of_obc_segments
1571 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
1572 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1573 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1574 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1575 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1576 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1577
1578 cs%initialized = .true.
1579 in_use = .false. ! Set to true to avoid deallocating
1580 cs%diag => diag ! Diagnostics pointer
1581 cs%calculate_cg1 = .false.
1582 cs%calculate_Rd_dx = .false.
1583 cs%calculate_res_fns = .false.
1584 cs%use_simpler_Eady_growth_rate = .false.
1585 cs%full_depth_Eady_growth_rate = .false.
1586 cs%calculate_depth_fns = .false.
1587 ! Read all relevant parameters and write them to the model log.
1588 call log_version(param_file, mdl, version, "")
1589 call get_param(param_file, mdl, "USE_VARIABLE_MIXING", cs%use_variable_mixing,&
1590 "If true, the variable mixing code will be called. This "//&
1591 "allows diagnostics to be created even if the scheme is "//&
1592 "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
1593 "this is set to true regardless of what is in the "//&
1594 "parameter file.", default=.false.)
1595 call get_param(param_file, mdl, "USE_VISBECK", cs%use_Visbeck,&
1596 "If true, use the Visbeck et al. (1997) formulation for \n"//&
1597 "thickness diffusivity.", default=.false.)
1598 call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
1599 "If true, the Laplacian lateral viscosity is scaled away "//&
1600 "when the first baroclinic deformation radius is well "//&
1601 "resolved.", default=.false.)
1602 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", cs%Depth_scaled_KhTh, &
1603 "If true, KHTH is scaled away when the depth is shallower "//&
1604 "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//&
1605 "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
1606 "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
1607 default=.false.)
1608 call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
1609 "If true, the interface depth diffusivity is scaled away "//&
1610 "when the first baroclinic deformation radius is well "//&
1611 "resolved.", default=.false.)
1612 call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
1613 "If true, the epipycnal tracer diffusivity is scaled "//&
1614 "away when the first baroclinic deformation radius is "//&
1615 "well resolved.", default=.false.)
1616 call get_param(param_file, mdl, "USE_MEKE", use_meke, &
1617 default=.false., do_not_log=.true.)
1618 call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", resoln_scaled_meke_visc, &
1619 "If true, the viscosity contribution from MEKE is scaled by "//&
1620 "the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
1621 if (.not.use_meke) resoln_scaled_meke_visc = .false.
1622 call get_param(param_file, mdl, "RESOLN_USE_EBT", cs%Resoln_use_ebt, &
1623 "If true, uses the equivalent barotropic wave speed instead "//&
1624 "of first baroclinic wave for calculating the resolution function.",&
1625 default=.false.)
1626 call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", cs%BS_EBT_power, &
1627 "Power to raise EBT vertical structure to when backscatter "// &
1628 "has vertical structure.", units="nondim", default=0.0)
1629 call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", cs%BS_use_sqg_struct, &
1630 "If true, the SQG vertical structure is used for backscatter "//&
1631 "on the condition that BS_EBT_power=0", &
1632 default=.false.)
1633 call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
1634 "If true, uses the equivalent barotropic structure "//&
1635 "as the vertical structure of thickness diffusivity.",&
1636 default=.false.)
1637 call get_param(param_file, mdl, "KHTH_USE_SQG_STRUCT", cs%khth_use_sqg_struct, &
1638 "If true, uses the surface quasigeostrophic structure "//&
1639 "as the vertical structure of thickness diffusivity.",&
1640 default=.false.)
1641 call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", cs%khtr_use_ebt_struct, &
1642 "If true, uses the equivalent barotropic structure "//&
1643 "as the vertical structure of tracer diffusivity.",&
1644 default=.false.)
1645 call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", cs%khtr_use_sqg_struct, &
1646 "If true, uses the surface quasigeostrophic structure "//&
1647 "as the vertical structure of tracer diffusivity.",&
1648 default=.false.)
1649 call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", cs%kdgl90_use_ebt_struct, &
1650 "If true, uses the equivalent barotropic structure "//&
1651 "as the vertical structure of diffusivity in the GL90 scheme.",&
1652 default=.false.)
1653 call get_param(param_file, mdl, "KD_GL90_USE_SQG_STRUCT", cs%kdgl90_use_sqg_struct, &
1654 "If true, uses the equivalent barotropic structure "//&
1655 "as the vertical structure of diffusivity in the GL90 scheme.",&
1656 default=.false.)
1657 call get_param(param_file, mdl, "KHTH_SLOPE_CFF", khth_slope_cff, &
1658 "The nondimensional coefficient in the Visbeck formula "//&
1659 "for the interface depth diffusivity", units="nondim", default=0.0)
1660 call get_param(param_file, mdl, "KHTR_SLOPE_CFF", khtr_slope_cff, &
1661 "The nondimensional coefficient in the Visbeck formula "//&
1662 "for the epipycnal tracer diffusivity", units="nondim", default=0.0)
1663 call get_param(param_file, mdl, "USE_STORED_SLOPES", cs%use_stored_slopes,&
1664 "If true, the isopycnal slopes are calculated once and "//&
1665 "stored for re-use. This uses more memory but avoids calling "//&
1666 "the equation of state more times than should be necessary.", &
1667 default=.false.)
1668 call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1669 "A miniscule frequency that is used to avoid division by 0. The default "//&
1670 "value is roughly (pi / (the age of the universe)).", &
1671 default=1.0e-17, units="s-1", scale=us%T_to_s)
1672 call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
1673 default=.false., do_not_log=.true.)
1674 cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn .or. cs%khth_use_ebt_struct &
1675 .or. cs%kdgl90_use_ebt_struct .or. cs%BS_EBT_power>0.
1676 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
1677 ! Indicate whether to calculate the Eady growth rate
1678 cs%calculate_Eady_growth_rate = use_meke .or. (khtr_slope_cff>0.) .or. (khth_slope_cff>0.)
1679 call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
1680 units="nondim", default=0., do_not_log=.true.)
1681 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
1682 call get_param(param_file, mdl, "MLE_FRONT_LENGTH", mle_front_length, &
1683 units="m", default=0.0, scale=us%m_to_L, do_not_log=.true.)
1684 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
1685
1686 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
1687
1688 call get_param(param_file, mdl, "STOCH_EOS", stoch_eos, &
1689 default=.false., do_not_log=.true.)
1690 call get_param(param_file, mdl, "USE_STANLEY_ISO", cs%use_stanley_iso, &
1691 "If true, turn on Stanley SGS T variance parameterization "// &
1692 "in isopycnal slope code.", default=.false.)
1693 if (cs%use_Stanley_ISO .and. .not.stoch_eos) then
1694 call mom_error(fatal, "VarMix_init: USE_STANLEY_ISO requires STOCH_EOS")
1695 endif
1696 call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_obc_segments, &
1697 default=0, do_not_log=.true.)
1698 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1699 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1700 call get_param(param_file, mdl, "MIXING_COEFS_OBC_BUG", mixing_coefs_obc_bug, &
1701 "If false, use only interior data for thickness weighting in lateral mixing "//&
1702 "coefficient calculations and to calculate stratification and other fields at "//&
1703 "open boundary condition faces.", &
1704 default=enable_bugs, do_not_log=(number_of_obc_segments<=0))
1705 cs%OBC_friendly = .not. mixing_coefs_obc_bug
1706 call get_param(param_file, mdl, "RESOLN_FUNCTION_OBC_BUG", cs%res_fn_OBC_bug, &
1707 "If false, use only interior data for calculating the resolution functions at "//&
1708 "open boundary condition faces and vertices.", &
1709 default=enable_bugs, do_not_log=(number_of_obc_segments<=0))
1710
1711 if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct .or. cs%kdgl90_use_ebt_struct &
1712 .or. cs%BS_EBT_power>0. .or. cs%khtr_use_ebt_struct) then
1713 in_use = .true.
1714 call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
1715 "The depth below which N2 is monotonized to avoid stratification "//&
1716 "artifacts from altering the equivalent barotropic mode structure. "//&
1717 "This monotonzization is disabled if this parameter is negative.", &
1718 units="m", default=-1.0, scale=gv%m_to_H)
1719 allocate(cs%ebt_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1720 endif
1721
1722 use_sqg = cs%BS_use_sqg_struct .or. cs%khth_use_sqg_struct .or. cs%khtr_use_sqg_struct .or. &
1723 cs%kdgl90_use_sqg_struct
1724 call get_param(param_file, mdl, "SQG_EXPO", cs%sqg_expo, &
1725 "Nondimensional exponent coeffecient of the SQG mode that is used for the "//&
1726 "vertical struture of diffusivities.", &
1727 units="nondim", default=1.0, do_not_log=.not.use_sqg)
1728 call get_param(param_file, mdl, "INTERPOLATED_SQG_STRUCTURE", cs%interpolated_sqg_struct, &
1729 "If true, interpolate properties to velocity points and then interpolate the "//&
1730 "buoyancy frequencies and layer thicknesses back to tracer points when "//&
1731 "calculating the SQG vertical structure.", &
1732 default=.true., do_not_log=.not.use_sqg)
1733 !### Consider changing the default for INTERPOLATED_SQG_STRUCTURE to false.
1734
1735 if ((cs%BS_EBT_power>0.) .and. cs%BS_use_sqg_struct) call mom_error(fatal, &
1736 "calc_resoln_function: BS_EBT_POWER>0. and BS_USE_SQG=True cannot be set together")
1737
1738 if (cs%khth_use_ebt_struct .and. cs%khth_use_sqg_struct) call mom_error(fatal, &
1739 "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT and KHTH_USE_SQG_STRUCT can be true")
1740
1741 if (cs%khtr_use_ebt_struct .and. cs%khtr_use_sqg_struct) call mom_error(fatal, &
1742 "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT and KHTR_USE_SQG_STRUCT can be true")
1743
1744 if (cs%kdgl90_use_ebt_struct .and. cs%kdgl90_use_sqg_struct) call mom_error(fatal, &
1745 "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT and KD_GL90_USE_SQG_STRUCT can be true")
1746
1747 if (cs%BS_EBT_power>0. .or. cs%BS_use_sqg_struct) then
1748 allocate(cs%BS_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1749 endif
1750
1751 if (cs%khth_use_ebt_struct .or. cs%khth_use_sqg_struct) then
1752 allocate(cs%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1753 endif
1754
1755 if (cs%khtr_use_ebt_struct .or. cs%khtr_use_sqg_struct) then
1756 allocate(cs%khtr_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1757 endif
1758
1759 if (cs%kdgl90_use_ebt_struct .or. cs%kdgl90_use_sqg_struct) then
1760 allocate(cs%kdgl90_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1761 endif
1762
1763 if (cs%use_stored_slopes) then
1764 if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1765 call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
1766 "If non-zero, is an upper bound on slopes used in the "//&
1767 "Visbeck formula for diffusivity. This does not affect the "//&
1768 "isopycnal slope calculation used within thickness diffusion.", &
1769 units="nondim", default=0.0, scale=us%L_to_Z)
1770 else
1771 cs%Visbeck_S_max = 0.
1772 endif
1773 endif
1774
1775 if (cs%use_stored_slopes .or. (cs%interpolated_sqg_struct .and. (cs%sqg_expo>0.0))) then
1776 ! CS%calculate_Eady_growth_rate=.true.
1777 in_use = .true.
1778 allocate(cs%slope_x(isdb:iedb,jsd:jed,gv%ke+1), source=0.0)
1779 allocate(cs%slope_y(isd:ied,jsdb:jedb,gv%ke+1), source=0.0)
1780 call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
1781 "A diapycnal diffusivity that is used to interpolate "//&
1782 "more sensible values of T & S into thin layers.", &
1783 units="m2 s-1", default=1.0e-6, scale=gv%m2_s_to_HZ_T)
1784 endif
1785
1786 if (cs%calculate_Eady_growth_rate) then
1787 in_use = .true.
1788 allocate(cs%SN_u(isdb:iedb,jsd:jed), source=0.0)
1789 allocate(cs%SN_v(isd:ied,jsdb:jedb), source=0.0)
1790 cs%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, time, &
1791 'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=us%s_to_T)
1792 cs%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, time, &
1793 'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=us%s_to_T)
1794 call get_param(param_file, mdl, "USE_SIMPLER_EADY_GROWTH_RATE", cs%use_simpler_Eady_growth_rate, &
1795 "If true, use a simpler method to calculate the Eady growth rate "//&
1796 "that avoids division by layer thickness. Recommended.", default=.false.)
1797 if (cs%use_simpler_Eady_growth_rate) then
1798 if (.not. cs%use_stored_slopes) call mom_error(fatal, &
1799 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1800 "When USE_SIMPLER_EADY_GROWTH_RATE=True, USE_STORED_SLOPES must also be True.")
1801 call get_param(param_file, mdl, "EADY_GROWTH_RATE_D_SCALE", cs%Eady_GR_D_scale, &
1802 "The depth from surface over which to average SN when calculating "//&
1803 "a 2D Eady growth rate. Zero mean use full depth.", &
1804 units="m", default=0., scale=us%m_to_Z)
1805 call get_param(param_file, mdl, "EADY_GROWTH_RATE_CROPPING_DISTANCE", cs%cropping_distance, &
1806 "Distance from surface or bottom to filter out outcropped or "//&
1807 "incropped interfaces for the Eady growth rate calc. "//&
1808 "Negative values disables cropping.", units="m", default=0., scale=us%m_to_Z)
1809 else
1810 call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
1811 "The layer number at which to start vertical integration "//&
1812 "of S*N for purposes of finding the Eady growth rate.", &
1813 units="nondim", default=2)
1814 call get_param(param_file, mdl, "MIN_DZ_FOR_SLOPE_N2", cs%h_min_N2, &
1815 "The minimum vertical distance to use in the denominator of the "//&
1816 "bouyancy frequency used in the slope calculation.", &
1817 units="m", default=1.0, scale=gv%m_to_H, do_not_log=cs%use_stored_slopes)
1818
1819 call get_param(param_file, mdl, "FULL_DEPTH_EADY_GROWTH_RATE", cs%full_depth_Eady_growth_rate, &
1820 "If true, calculate the Eady growth rate based on average slope times "//&
1821 "stratification that includes contributions from sea-level changes "//&
1822 "in its denominator, rather than just the nominal depth of the bathymetry. "//&
1823 "This only applies when using the model interface heights as a proxy for "//&
1824 "isopycnal slopes.", default=.not.(gv%Boussinesq.or.gv%semi_Boussinesq), &
1825 do_not_log=cs%use_stored_slopes)
1826 endif
1827 endif
1828
1829 if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1830 in_use = .true.
1831 call get_param(param_file, mdl, "VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1832 "The fixed length scale in the Visbeck formula, or if negative a nondimensional "//&
1833 "scaling factor relating this length scale squared to the cell areas.", &
1834 units="m or nondim", default=0.0, scale=us%m_to_L)
1835 allocate(cs%L2u(isdb:iedb,jsd:jed), source=0.0)
1836 allocate(cs%L2v(isd:ied,jsdb:jedb), source=0.0)
1837 if (cs%Visbeck_L_scale<0) then
1838 ! Undo the rescaling of CS%Visbeck_L_scale.
1839 do j=js,je ; do i=is-1,ieq
1840 cs%L2u(i,j) = (us%L_to_m*cs%Visbeck_L_scale)**2 * g%areaCu(i,j)
1841 enddo ; enddo
1842 do j=js-1,jeq ; do i=is,ie
1843 cs%L2v(i,j) = (us%L_to_m*cs%Visbeck_L_scale)**2 * g%areaCv(i,j)
1844 enddo ; enddo
1845 else
1846 cs%L2u(:,:) = cs%Visbeck_L_scale**2
1847 cs%L2v(:,:) = cs%Visbeck_L_scale**2
1848 endif
1849
1850 cs%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, time, &
1851 'Length scale squared for mixing coefficient, at u-points', &
1852 'm2', conversion=us%L_to_m**2)
1853 cs%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, time, &
1854 'Length scale squared for mixing coefficient, at v-points', &
1855 'm2', conversion=us%L_to_m**2)
1856 endif
1857
1858 cs%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, time, &
1859 'Vertical structure of SQG mode', 'nondim')
1860 if (cs%BS_use_sqg_struct .or. cs%khth_use_sqg_struct .or. cs%khtr_use_sqg_struct &
1861 .or. cs%kdgl90_use_sqg_struct .or. cs%id_sqg_struct>0) then
1862 allocate(cs%sqg_struct(isd:ied,jsd:jed,gv%ke), source=0.0)
1863 endif
1864
1865 if (cs%BS_EBT_power>0. .or. cs%BS_use_sqg_struct) then
1866 cs%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, time, &
1867 'Vertical structure of backscatter', 'nondim')
1868 endif
1869 if (cs%khth_use_ebt_struct .or. cs%khth_use_sqg_struct) then
1870 cs%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, time, &
1871 'Vertical structure of thickness diffusivity', 'nondim')
1872 endif
1873 if (cs%khtr_use_ebt_struct .or. cs%khtr_use_sqg_struct) then
1874 cs%id_khtr_struct = register_diag_field('ocean_model', 'khtr_struct', diag%axesTl, time, &
1875 'Vertical structure of tracer diffusivity', 'nondim')
1876 endif
1877 if (cs%kdgl90_use_ebt_struct .or. cs%kdgl90_use_sqg_struct) then
1878 cs%id_kdgl90_struct = register_diag_field('ocean_model', 'kdgl90_struct', diag%axesTl, time, &
1879 'Vertical structure of GL90 diffusivity', 'nondim')
1880 endif
1881
1882 if ((cs%calculate_Eady_growth_rate .and. cs%use_stored_slopes) ) then
1883 cs%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, time, &
1884 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
1885 's-2', conversion=(us%L_to_Z*us%s_to_T)**2)
1886 cs%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, time, &
1887 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', &
1888 's-2', conversion=(us%L_to_Z*us%s_to_T)**2)
1889 endif
1890 if (cs%use_simpler_Eady_growth_rate) then
1891 cs%id_dzu = register_diag_field('ocean_model', 'dzu_Visbeck', diag%axesCui, time, &
1892 'dz at u-points, used in calculating Eady growth rate in Visbeck et al..', &
1893 'm', conversion=us%Z_to_m)
1894 cs%id_dzv = register_diag_field('ocean_model', 'dzv_Visbeck', diag%axesCvi, time, &
1895 'dz at v-points, used in calculating Eady growth rate in Visbeck et al..', &
1896 'm', conversion=us%Z_to_m)
1897 cs%id_dzSxN = register_diag_field('ocean_model', 'dzSxN', diag%axesCui, time, &
1898 'dz * |slope_x| * N, used in calculating Eady growth rate in '//&
1899 'Visbeck et al..', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1900 cs%id_dzSyN = register_diag_field('ocean_model', 'dzSyN', diag%axesCvi, time, &
1901 'dz * |slope_y| * N, used in calculating Eady growth rate in '//&
1902 'Visbeck et al..', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1903 endif
1904 if (cs%use_stored_slopes) then
1905 cs%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, time, &
1906 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', &
1907 'nondim', conversion=us%Z_to_L**2)
1908 cs%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, time, &
1909 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', &
1910 'nondim', conversion=us%Z_to_L**2)
1911 endif
1912
1913 oneortwo = 1.0
1914 cs%Resoln_scaling_used = cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. &
1915 cs%Resoln_scaled_KhTr .or. resoln_scaled_meke_visc
1916 if (cs%Resoln_scaling_used) then
1917 cs%calculate_Rd_dx = .true.
1918 cs%calculate_res_fns = .true.
1919 allocate(cs%Res_fn_h(isd:ied,jsd:jed), source=0.0)
1920 allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb), source=0.0)
1921 allocate(cs%Res_fn_u(isdb:iedb,jsd:jed), source=0.0)
1922 allocate(cs%Res_fn_v(isd:ied,jsdb:jedb), source=0.0)
1923 allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb), source=0.0)
1924 allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed), source=0.0)
1925 allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb), source=0.0)
1926 allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb), source=0.0)
1927 allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed), source=0.0)
1928 allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb), source=0.0)
1929
1930 cs%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, time, &
1931 'Resolution function for scaling diffusivities', 'nondim')
1932
1933 call get_param(param_file, mdl, "KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1934 "A coefficient that determines how KhTh is scaled away if "//&
1935 "RESOLN_SCALED_... is true, as "//&
1936 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1937 units="nondim", default=1.0)
1938 call get_param(param_file, mdl, "KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1939 "The power of dx/Ld in the Kh resolution function. Any "//&
1940 "positive integer may be used, although even integers "//&
1941 "are more efficient to calculate. Setting this greater "//&
1942 "than 100 results in a step-function being used.", &
1943 default=2)
1944 call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1945 "A coefficient that determines how Kh is scaled away if "//&
1946 "RESOLN_SCALED_... is true, as "//&
1947 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1948 "This function affects lateral viscosity, Kh, and not KhTh.", &
1949 units="nondim", default=cs%Res_coef_khth)
1950 call get_param(param_file, mdl, "VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1951 "The power of dx/Ld in the Kh resolution function. Any "//&
1952 "positive integer may be used, although even integers "//&
1953 "are more efficient to calculate. Setting this greater "//&
1954 "than 100 results in a step-function being used. "//&
1955 "This function affects lateral viscosity, Kh, and not KhTh.", &
1956 default=cs%Res_fn_power_khth)
1957 call get_param(param_file, mdl, "INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1958 "If true, interpolate the resolution function to the "//&
1959 "velocity points from the thickness points; otherwise "//&
1960 "interpolate the wave speed and calculate the resolution "//&
1961 "function independently at each point.", default=.false.)
1962 if (cs%interpolate_Res_fn) then
1963 if (cs%Res_coef_visc /= cs%Res_coef_khth) call mom_error(fatal, &
1964 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1965 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1966 if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth) call mom_error(fatal, &
1967 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1968 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1969 endif
1970 call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1971 "If true, uses Gill's definition of the baroclinic "//&
1972 "equatorial deformation radius, otherwise, if false, use "//&
1973 "Pedlosky's definition. These definitions differ by a factor "//&
1974 "of 2 in front of the beta term in the denominator. Gill's "//&
1975 "is the more appropriate definition.", default=.true.)
1976 if (gill_equatorial_ld) then
1977 oneortwo = 2.0
1978 endif
1979
1980 do j=js-1,jeq ; do i=is-1,ieq
1981 cs%f2_dx2_q(i,j) = ((g%dxBu(i,j)**2) + (g%dyBu(i,j)**2)) * &
1982 max(g%Coriolis2Bu(i,j), absurdly_small_freq**2)
1983 cs%beta_dx2_q(i,j) = oneortwo * ((g%dxBu(i,j)**2) + (g%dyBu(i,j)**2)) * (sqrt(0.5 * &
1984 ( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) + &
1985 (((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2)) + &
1986 ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
1987 (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2)) ) ))
1988 enddo ; enddo
1989
1990 do j=js,je ; do i=is-1,ieq
1991 cs%f2_dx2_u(i,j) = ((g%dxCu(i,j)**2) + (g%dyCu(i,j)**2)) * &
1992 max(0.5* (g%Coriolis2Bu(i,j)+g%Coriolis2Bu(i,j-1)), absurdly_small_freq**2)
1993 cs%beta_dx2_u(i,j) = oneortwo * ((g%dxCu(i,j)**2) + (g%dyCu(i,j)**2)) * (sqrt( &
1994 ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1995 0.25*( ((((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1996 (((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2)) + &
1997 ((((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2) + &
1998 (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2)) ) ))
1999 enddo ; enddo
2000
2001 do j=js-1,jeq ; do i=is,ie
2002 cs%f2_dx2_v(i,j) = ((g%dxCv(i,j)**2) + (g%dyCv(i,j)**2)) * &
2003 max(0.5*(g%Coriolis2Bu(i,j)+g%Coriolis2Bu(i-1,j)), absurdly_small_freq**2)
2004 cs%beta_dx2_v(i,j) = oneortwo * ((g%dxCv(i,j)**2) + (g%dyCv(i,j)**2)) * (sqrt( &
2005 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
2006 0.25*( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
2007 (((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2)) + &
2008 ((((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) + &
2009 (((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2)) ) ))
2010 enddo ; enddo
2011
2012 endif
2013
2014 if (cs%Depth_scaled_KhTh) then
2015 cs%calculate_depth_fns = .true.
2016 allocate(cs%Depth_fn_u(isdb:iedb,jsd:jed), source=0.0)
2017 allocate(cs%Depth_fn_v(isd:ied,jsdb:jedb), source=0.0)
2018 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", cs%depth_scaled_khth_h0, &
2019 "The depth above which KHTH is scaled away.", &
2020 units="m", scale=us%m_to_Z, default=1000.)
2021 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", cs%depth_scaled_khth_exp, &
2022 "The exponent used in the depth dependent scaling function for KHTH.", &
2023 units="nondim", default=3.0)
2024 endif
2025
2026 ! Resolution %Rd_dx_h
2027 cs%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, time, &
2028 'Ratio between deformation radius and grid spacing', 'm m-1')
2029 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
2030
2031 if (cs%calculate_Rd_dx) then
2032 cs%calculate_cg1 = .true. ! We will need %cg1
2033 allocate(cs%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
2034 allocate(cs%beta_dx2_h(isd:ied,jsd:jed), source=0.0)
2035 allocate(cs%f2_dx2_h(isd:ied,jsd:jed), source=0.0)
2036 do j=js-1,je+1 ; do i=is-1,ie+1
2037 cs%f2_dx2_h(i,j) = ((g%dxT(i,j)**2) + (g%dyT(i,j)**2)) * &
2038 max(0.25 * ((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
2039 (g%Coriolis2Bu(i-1,j) + g%Coriolis2Bu(i,j-1))), &
2040 absurdly_small_freq**2)
2041 cs%beta_dx2_h(i,j) = oneortwo * ((g%dxT(i,j)**2) + (g%dyT(i,j)**2)) * (sqrt(0.5 * &
2042 ( ((((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) + &
2043 (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2)) + &
2044 ((((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2) + &
2045 (((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2)) ) ))
2046 enddo ; enddo
2047 endif
2048
2049 if (cs%calculate_cg1) then
2050 in_use = .true.
2051 allocate(cs%cg1(isd:ied,jsd:jed), source=0.0)
2052 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2053 "This sets the default value for the various _ANSWER_DATE parameters.", &
2054 default=99991231)
2055 call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2056 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
2057 "Values below 20190101 result in the use of older, less accurate expressions "//&
2058 "that were in use at the end of 2018. Higher values result in the use of more "//&
2059 "robust and accurate forms of mathematically equivalent expressions.", &
2060 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
2061 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2062
2063 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
2064 "The fractional tolerance for finding the wave speeds.", &
2065 units="nondim", default=0.001)
2066 !### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
2067 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
2068 "A floor in the first mode speed below which 0 used instead.", &
2069 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
2070 call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
2071 "If true, use a more robust estimate of the first mode wave speed as the "//&
2072 "starting point for iterations.", default=.true.)
2073 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2074 do_not_log=.true., default=.true.)
2075 call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2076 "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
2077 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
2078 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2079 call wave_speed_init(cs%wave_speed, gv, use_ebt_mode=cs%Resoln_use_ebt, &
2080 mono_n2_depth=n2_filter_depth, remap_answer_date=remap_answer_date, &
2081 better_speed_est=better_speed_est, min_speed=wave_speed_min, &
2082 om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol)
2083 endif
2084
2085 ! Leith parameters
2086 call get_param(param_file, mdl, "USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
2087 "If true, use the QG Leith viscosity as the GM coefficient.", &
2088 default=.false.)
2089
2090 if (cs%Use_QG_Leith_GM) then
2091 call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
2092 "The nondimensional Laplacian Leith constant, \n"//&
2093 "often set to 1.0", units="nondim", default=0.0)
2094
2095 call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
2096 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
2097 default=.true.)
2098
2099 allocate(cs%Laplac3_const_u(isdb:iedb,jsd:jed), source=0.0)
2100 allocate(cs%Laplac3_const_v(isd:ied,jsdb:jedb), source=0.0)
2101 allocate(cs%KH_u_QG(isdb:iedb,jsd:jed,gv%ke), source=0.0)
2102 allocate(cs%KH_v_QG(isd:ied,jsdb:jedb,gv%ke), source=0.0)
2103
2104 ! register diagnostics
2105 cs%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, time, &
2106 'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2107 cs%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, time, &
2108 'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2109
2110 do j=jsq,jeq+1 ; do i=is-1,ieq
2111 ! Static factors in the Leith schemes
2112 grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
2113 grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
2114 cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
2115 enddo ; enddo
2116 do j=js-1,jeq ; do i=isq,ieq+1
2117 ! Static factors in the Leith schemes
2118 grid_sp_v2 = g%dyCv(i,j)*g%dxCv(i,j)
2119 grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
2120 cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
2121 enddo ; enddo
2122
2123 if (.not. cs%use_stored_slopes) call mom_error(fatal, &
2124 "MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2125 "USE_STORED_SLOPES must be True when using QG Leith.")
2126 endif
2127
2128 ! Re-enable variable mixing if one of the schemes was enabled
2129 cs%use_variable_mixing = in_use .or. cs%use_variable_mixing
2130end subroutine varmix_init
2131
2132!> Destructor for VarMix control structure
2133subroutine varmix_end(CS)
2134 type(varmix_cs), intent(inout) :: cs
2135
2136 if (allocated(cs%ebt_struct)) deallocate(cs%ebt_struct)
2137 if (allocated(cs%sqg_struct)) deallocate(cs%sqg_struct)
2138 if (allocated(cs%BS_struct)) deallocate(cs%BS_struct)
2139 if (allocated(cs%khth_struct)) deallocate(cs%khth_struct)
2140 if (allocated(cs%khtr_struct)) deallocate(cs%khtr_struct)
2141 if (allocated(cs%kdgl90_struct)) deallocate(cs%kdgl90_struct)
2142
2143 if (allocated(cs%slope_x)) deallocate(cs%slope_x)
2144 if (allocated(cs%slope_y)) deallocate(cs%slope_y)
2145
2146 if (allocated(cs%SN_u)) deallocate(cs%SN_u)
2147 if (allocated(cs%SN_v)) deallocate(cs%SN_v)
2148
2149 if (allocated(cs%L2u)) deallocate(cs%L2u)
2150 if (allocated(cs%L2v)) deallocate(cs%L2v)
2151
2152 if (allocated(cs%Res_fn_h)) deallocate(cs%Res_fn_h)
2153 if (allocated(cs%Res_fn_q)) deallocate(cs%Res_fn_q)
2154 if (allocated(cs%Res_fn_u)) deallocate(cs%Res_fn_u)
2155 if (allocated(cs%Res_fn_v)) deallocate(cs%Res_fn_v)
2156 if (allocated(cs%beta_dx2_q)) deallocate(cs%beta_dx2_q)
2157 if (allocated(cs%beta_dx2_u)) deallocate(cs%beta_dx2_u)
2158 if (allocated(cs%beta_dx2_v)) deallocate(cs%beta_dx2_v)
2159 if (allocated(cs%f2_dx2_q)) deallocate(cs%f2_dx2_q)
2160 if (allocated(cs%f2_dx2_u)) deallocate(cs%f2_dx2_u)
2161 if (allocated(cs%f2_dx2_v)) deallocate(cs%f2_dx2_v)
2162
2163 if (allocated(cs%Depth_fn_u)) deallocate(cs%Depth_fn_u)
2164 if (allocated(cs%Depth_fn_v)) deallocate(cs%Depth_fn_v)
2165
2166 if (allocated(cs%Rd_dx_h)) deallocate(cs%Rd_dx_h)
2167 if (allocated(cs%beta_dx2_h)) deallocate(cs%beta_dx2_h)
2168 if (allocated(cs%f2_dx2_h)) deallocate(cs%f2_dx2_h)
2169
2170 if (allocated(cs%cg1)) deallocate(cs%cg1)
2171
2172 if (allocated(cs%Laplac3_const_u)) deallocate(cs%Laplac3_const_u)
2173 if (allocated(cs%Laplac3_const_v)) deallocate(cs%Laplac3_const_v)
2174 if (allocated(cs%KH_u_QG)) deallocate(cs%KH_u_QG)
2175 if (allocated(cs%KH_v_QG)) deallocate(cs%KH_v_QG)
2176
2177end subroutine varmix_end
2178
2179!> \namespace mom_lateral_mixing_coeffs
2180!!
2181!! This module provides a container for various factors used in prescribing diffusivities, that are
2182!! a function of the state (in particular the stratification and isoneutral slopes).
2183!!
2184!! \section section_Resolution_Function The resolution function
2185!!
2186!! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
2187!! The square of the resolution parameter is
2188!!
2189!! \f[
2190!! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
2191!! \f]
2192!!
2193!! where the grid spacing is calculated as
2194!!
2195!! \f[
2196!! \Delta^2 = \Delta x^2 + \Delta y^2 .
2197!! \f]
2198!!
2199!! \todo Check this reference to Bob on/off paper.
2200!! The resolution function used in scaling diffusivities (\cite hallberg2013) is
2201!!
2202!! \f[
2203!! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
2204!! \f]
2205!!
2206!! The resolution function can be applied independently to thickness diffusion \(module mom_thickness_diffuse\),
2207!! tracer diffusion \(mom_tracer_hordiff\) lateral viscosity \(mom_hor_visc\).
2208!!
2209!! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
2210!! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007
2211!!
2212!! | Symbol | Module parameter |
2213!! | ------ | --------------- |
2214!! | - | <code>USE_VARIABLE_MIXING</code> |
2215!! | - | <code>RESOLN_SCALED_KH</code> |
2216!! | - | <code>RESOLN_SCALED_KHTH</code> |
2217!! | - | <code>RESOLN_SCALED_KHTR</code> |
2218!! | \f$ \alpha \f$ | <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) |
2219!! | \f$ p \f$ | <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) |
2220!! | \f$ \alpha \f$ | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) |
2221!! | \f$ p \f$ | <code>VISC_RES_FN_POWER</code> (for lateral viscosity) |
2222!! | - | <code>GILL_EQUATORIAL_LD</code> |
2223!!
2224!!
2225!!
2226!! \section section_Vicbeck Visbeck diffusivity
2227!!
2228!! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
2229!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse but calculated in this module.
2230!!
2231!! \f[
2232!! \kappa_h = \alpha_s L_s^2 S N
2233!! \f]
2234!!
2235!! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
2236!!
2237!! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
2238!! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
2239!!
2240!! | Symbol | Module parameter |
2241!! | ------ | --------------- |
2242!! | - | <code>USE_VARIABLE_MIXING</code> |
2243!! | \f$ \alpha_s \f$ | <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)|
2244!! | \f$ \alpha_s \f$ | <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)|
2245!! | \f$ L_{s} \f$ | <code>VISBECK_L_SCALE</code> |
2246!! | \f$ S_{max} \f$ | <code>VISBECK_MAX_SLOPE</code> |
2247!!
2248!!
2249!! \section section_vertical_structure_khth Vertical structure function for KhTh
2250!!
2251!! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
2252!! velocity mode. The structure function is stored in the control structure for this module (varmix_cs) but is
2253!! calculated using subroutines in mom_wave_speed.
2254!!
2255!! | Symbol | Module parameter |
2256!! | ------ | --------------- |
2257!! | - | <code>KHTH_USE_EBT_STRUCT</code> |
2258