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