MOM_CVMix_shear.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!> Interface to CVMix interior shear schemes
7
8!> \author Brandon Reichl
9
10use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
11use mom_diag_mediator, only : diag_ctrl, time_type
12use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
13use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_grid, only : ocean_grid_type
15use mom_interface_heights, only : thickness_to_dz
19use mom_eos, only : calculate_density
20use cvmix_shear, only : cvmix_init_shear, cvmix_coeffs_shear
22implicit none ; private
23
24#include <MOM_memory.h>
25
27
28! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
29! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
30! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
31! vary with the Boussinesq approximation, the Boussinesq variant is given first.
32
33!> Control structure including parameters for CVMix interior shear schemes.
34type, public :: cvmix_shear_cs ; private
35 logical :: use_lmd94 !< Flags to use the LMD94 scheme
36 logical :: use_pp81 !< Flags to use Pacanowski and Philander (JPO 1981)
37 integer :: n_smooth_ri !< Number of times to smooth Ri using a 1-2-1 filter
38 real :: ri_zero !< LMD94 critical Richardson number [nondim]
39 real :: nu_zero !< LMD94 maximum interior diffusivity [Z2 T-1 ~> m2 s-1]
40 real :: kpp_exp !< Exponent of unitless factor of diffusivities
41 !! for KPP internal shear mixing scheme [nondim]
42 real, allocatable, dimension(:,:,:) :: n2 !< Squared Brunt-Vaisala frequency [T-2 ~> s-2]
43 real, allocatable, dimension(:,:,:) :: s2 !< Squared shear frequency [T-2 ~> s-2]
44 real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number [nondim]
45 real, allocatable, dimension(:,:,:) :: ri_grad_orig !< Gradient Richardson number
46 !! after smoothing [nondim]
47 character(10) :: mix_scheme !< Mixing scheme name (string)
48
49 type(diag_ctrl), pointer :: diag => null() !< Pointer to the diagnostics control structure
50 !>@{ Diagnostic handles
51 integer :: id_n2 = -1, id_s2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1
52 integer :: id_ri_grad_orig = -1
53 !>@}
54
55end type cvmix_shear_cs
56
57character(len=40) :: mdl = "MOM_CVMix_shear" !< This module's name.
58
59contains
60
61!> Subroutine for calculating (internal) vertical diffusivities/viscosities
62subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
63 type(ocean_grid_type), intent(in) :: g !< Grid structure.
64 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
65 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
66 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_h !< Initial zonal velocity on T points [L T-1 ~> m s-1]
67 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: v_h !< Initial meridional velocity on T
68 !! points [L T-1 ~> m s-1]
69 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
70 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
71 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: kd !< The vertical diffusivity at each interface
72 !! (not layer!) [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
73 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: kv !< The vertical viscosity at each interface
74 !! (not layer!) [H Z T-1 ~> m2 s-1 or Pa s]
75 type(cvmix_shear_cs), pointer :: cs !< The control structure returned by a previous
76 !! call to CVMix_shear_init.
77 ! Local variables
78 integer :: i, j, k, kk, km1, s
79 real :: gorho ! Gravitational acceleration divided by density [Z T-2 R-1 ~> m4 s-2 kg-1]
80 real :: pref ! Interface pressures [R L2 T-2 ~> Pa]
81 real :: du, dv ! Velocity differences [L T-1 ~> m s-1]
82 real :: dz_int ! Grid spacing around an interface [Z ~> m]
83 real :: n2 ! Buoyancy frequency at an interface [T-2 ~> s-2]
84 real :: s2 ! Shear squared at an interface [T-2 ~> s-2]
85 real :: dummy ! A dummy variable [nondim]
86 real :: drho ! Buoyancy differences [Z T-2 ~> m s-2]
87 real, dimension(SZI_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]
88 real, dimension(2*(GV%ke)) :: pres_1d ! A column of interface pressures [R L2 T-2 ~> Pa]
89 real, dimension(2*(GV%ke)) :: temp_1d ! A column of temperatures [C ~> degC]
90 real, dimension(2*(GV%ke)) :: salt_1d ! A column of salinities [S ~> ppt]
91 real, dimension(2*(GV%ke)) :: rho_1d ! A column of densities at interface pressures [R ~> kg m-3]
92 real, dimension(GV%ke+1) :: ri_grad !< Gradient Richardson number [nondim]
93 real, dimension(GV%ke+1) :: ri_grad_prev !< Gradient Richardson number before s.th smoothing iteration [nondim]
94 real, dimension(GV%ke+1) :: kvisc !< Vertical viscosity at interfaces [m2 s-1]
95 real, dimension(GV%ke+1) :: kdiff !< Diapycnal diffusivity at interfaces [m2 s-1]
96 real :: epsln !< Threshold to identify vanished layers [H ~> m or kg m-2]
97
98 ! some constants
99 gorho = gv%g_Earth_Z_T2 / gv%Rho0
100 epsln = 1.e-10 * gv%m_to_H
101
102 do j = g%jsc, g%jec
103
104 ! Find the vertical distances across layers.
105 call thickness_to_dz(h, tv, dz, j, g, gv)
106
107 do i = g%isc, g%iec
108
109 ! skip calling for land points
110 if (g%mask2dT(i,j)==0.) cycle
111
112 ! Richardson number computed for each cell in a column.
113 pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
114 ri_grad(:)=1.e8 !Initialize w/ large Richardson value
115 do k=1,gv%ke
116 ! pressure, temp, and saln for EOS
117 ! kk+1 = k fields
118 ! kk+2 = km1 fields
119 km1 = max(1, k-1)
120 kk = 2*(k-1)
121 pres_1d(kk+1) = pref
122 pres_1d(kk+2) = pref
123 temp_1d(kk+1) = tv%T(i,j,k)
124 temp_1d(kk+2) = tv%T(i,j,km1)
125 salt_1d(kk+1) = tv%S(i,j,k)
126 salt_1d(kk+2) = tv%S(i,j,km1)
127
128 ! pRef is pressure at interface between k and km1.
129 ! iterate pRef for next pass through k-loop.
130 pref = pref + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k)
131
132 enddo ! k-loop finishes
133
134 ! compute in-situ density [R ~> kg m-3]
135 call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, tv%eqn_of_state)
136
137 ! N2 (can be negative) on interface
138 do k = 1, gv%ke
139 km1 = max(1, k-1)
140 kk = 2*(k-1)
141 du = u_h(i,j,k) - u_h(i,j,km1)
142 dv = v_h(i,j,k) - v_h(i,j,km1)
143 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
144 drho = gorho * (rho_1d(kk+1) - rho_1d(kk+2))
145 else
146 drho = gv%g_Earth_Z_T2 * (rho_1d(kk+1) - rho_1d(kk+2)) / (0.5*(rho_1d(kk+1) + rho_1d(kk+2)))
147 endif
148 dz_int = 0.5*(dz(i,km1) + dz(i,k)) + gv%dZ_subroundoff
149 n2 = drho / dz_int
150 s2 = us%L_to_Z**2*((du*du) + (dv*dv)) / (dz_int*dz_int)
151 ri_grad(k) = max(0., n2) / max(s2, 1.e-10*us%T_to_s**2)
152
153 ! fill 3d arrays, if user asks for diagnostics
154 if (cs%id_N2 > 0) cs%N2(i,j,k) = n2
155 if (cs%id_S2 > 0) cs%S2(i,j,k) = s2
156
157 enddo
158
159 ri_grad(gv%ke+1) = ri_grad(gv%ke)
160
161 if (cs%n_smooth_ri > 0) then
162
163 if (cs%id_ri_grad_orig > 0) cs%ri_grad_orig(i,j,:) = ri_grad(:)
164
165 ! 1) fill Ri_grad in vanished layers with adjacent value
166 do k = 2, gv%ke
167 if (h(i,j,k) <= epsln) ri_grad(k) = ri_grad(k-1)
168 enddo
169
170 ri_grad(gv%ke+1) = ri_grad(gv%ke)
171
172 do s=1,cs%n_smooth_ri
173
174 ri_grad_prev(:) = ri_grad(:)
175
176 ! 2) vertically smooth Ri with 1-2-1 filter
177 dummy = 0.25 * ri_grad_prev(2)
178 do k = 3, gv%ke
179 ri_grad(k) = dummy + 0.5 * ri_grad_prev(k) + 0.25 * ri_grad_prev(k+1)
180 dummy = 0.25 * ri_grad(k)
181 enddo
182 enddo
183
184 ri_grad(gv%ke+1) = ri_grad(gv%ke)
185
186 endif
187
188 if (cs%id_ri_grad > 0) cs%ri_grad(i,j,:) = ri_grad(:)
189
190 do k=1,gv%ke+1
191 kvisc(k) = gv%HZ_T_to_m2_s * kv(i,j,k)
192 kdiff(k) = gv%HZ_T_to_m2_s * kd(i,j,k)
193 enddo
194
195 ! Call to CVMix wrapper for computing interior mixing coefficients.
196 call cvmix_coeffs_shear(mdiff_out=kvisc(:), &
197 tdiff_out=kdiff(:), &
198 rich=ri_grad(:), &
199 nlev=gv%ke, &
200 max_nlev=gv%ke)
201 do k=1,gv%ke+1
202 kv(i,j,k) = gv%m2_s_to_HZ_T * kvisc(k)
203 kd(i,j,k) = gv%m2_s_to_HZ_T * kdiff(k)
204 enddo
205 enddo
206 enddo
207
208 ! write diagnostics
209 if (cs%id_kd > 0) call post_data(cs%id_kd, kd, cs%diag)
210 if (cs%id_kv > 0) call post_data(cs%id_kv, kv, cs%diag)
211 if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
212 if (cs%id_S2 > 0) call post_data(cs%id_S2, cs%S2, cs%diag)
213 if (cs%id_ri_grad > 0) call post_data(cs%id_ri_grad, cs%ri_grad, cs%diag)
214 if (cs%id_ri_grad_orig > 0) call post_data(cs%id_ri_grad_orig ,cs%ri_grad_orig, cs%diag)
215
216end subroutine calculate_cvmix_shear
217
218
219!> Initialized the CVMix internal shear mixing routine.
220!! \todo Does this note require emphasis?
221!! \note *This is where we test to make sure multiple internal shear
222!! mixing routines (including JHL) are not enabled at the same time.*
223!! (returns) CVMix_shear_init - True if module is to be used, False otherwise
224logical function cvmix_shear_init(Time, G, GV, US, param_file, diag, CS)
225 type(time_type), intent(in) :: time !< The current time.
226 type(ocean_grid_type), intent(in) :: g !< Grid structure.
227 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
228 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
229 type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
230 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
231 type(cvmix_shear_cs), pointer :: cs !< This module's control structure.
232 ! Local variables
233 integer :: numbertrue=0
234 logical :: use_jhl
235 logical :: use_lmd94
236 logical :: use_pp81
237
238! This include declares and sets the variable "version".
239#include "version_variable.h"
240
241 if (associated(cs)) then
242 call mom_error(warning, "CVMix_shear_init called with an associated "// &
243 "control structure.")
244 return
245 endif
246
247! Set default, read and log parameters
248 call get_param(param_file, mdl, "USE_LMD94", use_lmd94, default=.false., do_not_log=.true.)
249 call get_param(param_file, mdl, "USE_PP81", use_pp81, default=.false., do_not_log=.true.)
250 call log_version(param_file, mdl, version, &
251 "Parameterization of shear-driven turbulence via CVMix (various options)", &
252 all_default=.not.(use_pp81.or.use_lmd94))
253 call get_param(param_file, mdl, "USE_LMD94", use_lmd94, &
254 "If true, use the Large-McWilliams-Doney (JGR 1994) "//&
255 "shear mixing parameterization.", default=.false.)
256 if (use_lmd94) &
257 numbertrue=numbertrue + 1
258 call get_param(param_file, mdl, "USE_PP81", use_pp81, &
259 "If true, use the Pacanowski and Philander (JPO 1981) "//&
260 "shear mixing parameterization.", default=.false.)
261 if (use_pp81) &
262 numbertrue = numbertrue + 1
263 use_jhl=kappa_shear_is_used(param_file)
264 if (use_jhl) numbertrue = numbertrue + 1
265 ! After testing for interior schemes, make sure only 0 or 1 are enabled.
266 ! Otherwise, warn user and kill job.
267 if ((numbertrue) > 1) then
268 call mom_error(fatal, 'MOM_CVMix_shear_init: '// &
269 'Multiple shear driven internal mixing schemes selected, '//&
270 'please disable all but one scheme to proceed.')
271 endif
272
273 cvmix_shear_init = use_pp81 .or. use_lmd94
274
275 ! Forego remainder of initialization if not using this scheme
276 if (.not. cvmix_shear_init) return
277
278 allocate(cs)
279 cs%use_LMD94 = use_lmd94
280 cs%use_PP81 = use_pp81
281 if (use_lmd94) &
282 cs%Mix_Scheme = 'KPP'
283 if (use_pp81) &
284 cs%Mix_Scheme = 'PP'
285
286 call get_param(param_file, mdl, "NU_ZERO", cs%Nu_Zero, &
287 "Leading coefficient in KPP shear mixing.", &
288 units="m2 s-1", default=5.e-3, scale=us%m2_s_to_Z2_T)
289 call get_param(param_file, mdl, "RI_ZERO", cs%Ri_Zero, &
290 "Critical Richardson for KPP shear mixing, "// &
291 "NOTE this the internal mixing and this is "// &
292 "not for setting the boundary layer depth.", &
293 units="nondim", default=0.8)
294 call get_param(param_file, mdl, "KPP_EXP", cs%KPP_exp, &
295 "Exponent of unitless factor of diffusivities, "// &
296 "for KPP internal shear mixing scheme.", &
297 units="nondim", default=3.0)
298 call get_param(param_file, mdl, "N_SMOOTH_RI", cs%n_smooth_ri, &
299 "If > 0, vertically smooth the Richardson "// &
300 "number by applying a 1-2-1 filter N_SMOOTH_RI times.", &
301 default=0)
302 call cvmix_init_shear(mix_scheme=cs%Mix_Scheme, &
303 kpp_nu_zero=us%Z2_T_to_m2_s*cs%Nu_Zero, &
304 kpp_ri_zero=cs%Ri_zero, &
305 kpp_exp=cs%KPP_exp)
306
307 ! Register diagnostics; allocation and initialization
308 cs%diag => diag
309
310 cs%id_N2 = register_diag_field('ocean_model', 'N2_shear', diag%axesTi, time, &
311 'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module', '1/s2', conversion=us%s_to_T**2)
312 if (cs%id_N2 > 0) then
313 allocate( cs%N2( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
314 endif
315
316 cs%id_S2 = register_diag_field('ocean_model', 'S2_shear', diag%axesTi, time, &
317 'Square of vertical shear used by MOM_CVMix_shear module','1/s2', conversion=us%s_to_T**2)
318 if (cs%id_S2 > 0) then
319 allocate( cs%S2( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
320 endif
321
322 cs%id_ri_grad = register_diag_field('ocean_model', 'ri_grad_shear', diag%axesTi, time, &
323 'Gradient Richarson number used by MOM_CVMix_shear module','nondim')
324 if (cs%id_ri_grad > 0) then !Initialize w/ large Richardson value
325 allocate( cs%ri_grad( szi_(g), szj_(g), szk_(gv)+1 ), source=1.e8 )
326 endif
327
328 if (cs%n_smooth_ri > 0) then
329 cs%id_ri_grad_orig = register_diag_field('ocean_model', 'ri_grad_shear_orig', &
330 diag%axesTi, time, &
331 'Original gradient Richarson number, before smoothing was applied. This is '//&
332 'part of the MOM_CVMix_shear module and only available when N_SMOOTH_RI > 0','nondim')
333 endif
334 if (cs%id_ri_grad_orig > 0 .or. cs%n_smooth_ri > 0) then !Initialize w/ large Richardson value
335 allocate( cs%ri_grad_orig( szi_(g), szj_(g), szk_(gv)+1 ), source=1.e8 )
336 endif
337
338 cs%id_kd = register_diag_field('ocean_model', 'kd_shear_CVMix', diag%axesTi, time, &
339 'Vertical diffusivity added by MOM_CVMix_shear module', 'm2/s', conversion=gv%HZ_T_to_m2_s)
340 cs%id_kv = register_diag_field('ocean_model', 'kv_shear_CVMix', diag%axesTi, time, &
341 'Vertical viscosity added by MOM_CVMix_shear module', 'm2/s', conversion=gv%HZ_T_to_m2_s)
342
343end function cvmix_shear_init
344
345!> Reads the parameters "USE_LMD94" and "USE_PP81" and returns true if either is true.
346!! This function allows other modules to know whether this parameterization will
347!! be used without needing to duplicate the log entry.
348logical function cvmix_shear_is_used(param_file)
349 type(param_file_type), intent(in) :: param_file !< Run-time parameter files handle.
350 ! Local variables
351 logical :: lmd94, pp81
352 call get_param(param_file, mdl, "USE_LMD94", lmd94, &
353 default=.false., do_not_log=.true.)
354 call get_param(param_file, mdl, "USE_PP81", pp81, &
355 default=.false., do_not_log=.true.)
356 cvmix_shear_is_used = (lmd94 .or. pp81)
357end function cvmix_shear_is_used
358
359!> Clear pointers and deallocate memory
360subroutine cvmix_shear_end(CS)
361 type(cvmix_shear_cs), intent(inout) :: cs !< Control structure for this module that
362 !! will be deallocated in this subroutine
363 if (cs%id_N2 > 0) deallocate(cs%N2)
364 if (cs%id_S2 > 0) deallocate(cs%S2)
365 if (cs%id_ri_grad > 0) deallocate(cs%ri_grad)
366end subroutine cvmix_shear_end
367
368end module mom_cvmix_shear