17 use cvmix_shear,
only : cvmix_init_shear, cvmix_coeffs_shear
19 implicit none ;
private
21 #include <MOM_memory.h>
39 real,
allocatable,
dimension(:,:,:) :: n2
40 real,
allocatable,
dimension(:,:,:) :: s2
41 real,
allocatable,
dimension(:,:,:) :: ri_grad
42 real,
allocatable,
dimension(:,:,:) :: ri_grad_smooth
44 character(10) :: mix_scheme
48 integer :: id_n2 = -1, id_s2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1
49 integer :: id_ri_grad_smooth = -1
54 character(len=40) ::
mdl =
"MOM_CVMix_shear"
63 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_h
64 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: v_h
66 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
68 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd
70 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kv
75 integer :: i, j, k, kk, km1
77 real :: pref, du, dv, drho, dz, n2, s2, dummy
78 real,
dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d
79 real,
dimension(G%ke+1) :: ri_grad
80 real,
dimension(G%ke+1) :: kvisc
81 real,
dimension(G%ke+1) :: kdiff
82 real,
parameter :: epsln = 1.e-10
85 gorho = gv%mks_g_Earth / (us%R_to_kg_m3*gv%Rho0)
91 if (g%mask2dT(i,j)==0.) cycle
104 temp_1d(kk+1) = tv%T(i,j,k)
105 temp_1d(kk+2) = tv%T(i,j,km1)
106 salt_1d(kk+1) = tv%S(i,j,k)
107 salt_1d(kk+2) = tv%S(i,j,km1)
111 pref = pref + gv%H_to_Pa * h(i,j,k)
116 call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 2*g%ke, tv%EQN_OF_STATE)
122 du = us%L_T_to_m_s*(u_h(i,j,k) - u_h(i,j,km1))
123 dv = us%L_T_to_m_s*(v_h(i,j,k) - v_h(i,j,km1))
124 drho = (gorho * (rho_1d(kk+1) - rho_1d(kk+2)) )
125 dz = ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
127 s2 = (du*du+dv*dv)/(dz*dz)
128 ri_grad(k) = max(0.,n2)/max(s2,1.e-10)
131 if (cs%id_N2 > 0) cs%N2(i,j,k) = n2
132 if (cs%id_S2 > 0) cs%S2(i,j,k) = s2
136 ri_grad(g%ke+1) = ri_grad(g%ke)
138 if (cs%id_ri_grad > 0) cs%ri_grad(i,j,:) = ri_grad(:)
140 if (cs%smooth_ri)
then
143 if (h(i,j,k) .le. epsln) ri_grad(k) = ri_grad(k-1)
146 ri_grad(g%ke+1) = ri_grad(g%ke)
149 dummy = 0.25 * ri_grad(2)
150 ri_grad(g%ke+1) = ri_grad(g%ke)
152 ri_grad(k) = dummy + 0.5 * ri_grad(k) + 0.25 * ri_grad(k+1)
153 dummy = 0.25 * ri_grad(k)
156 if (cs%id_ri_grad_smooth > 0) cs%ri_grad_smooth(i,j,:) = ri_grad(:)
160 kvisc(k) = us%Z2_T_to_m2_s * kv(i,j,k)
161 kdiff(k) = us%Z2_T_to_m2_s * kd(i,j,k)
165 call cvmix_coeffs_shear(mdiff_out=kvisc(:), &
166 tdiff_out=kdiff(:), &
171 kv(i,j,k) = us%m2_s_to_Z2_T * kvisc(k)
172 kd(i,j,k) = us%m2_s_to_Z2_T * kdiff(k)
178 if (cs%id_kd > 0)
call post_data(cs%id_kd, kd, cs%diag)
179 if (cs%id_kv > 0)
call post_data(cs%id_kv, kv, cs%diag)
180 if (cs%id_N2 > 0)
call post_data(cs%id_N2, cs%N2, cs%diag)
181 if (cs%id_S2 > 0)
call post_data(cs%id_S2, cs%S2, cs%diag)
182 if (cs%id_ri_grad > 0)
call post_data(cs%id_ri_grad, cs%ri_grad, cs%diag)
183 if (cs%id_ri_grad_smooth > 0)
call post_data(cs%id_ri_grad_smooth ,cs%ri_grad_smooth, cs%diag)
193 type(time_type),
intent(in) :: time
198 type(
diag_ctrl),
target,
intent(inout) :: diag
201 integer :: numbertrue=0
204 #include "version_variable.h"
206 if (
associated(cs))
then
207 call mom_error(warning,
"CVMix_shear_init called with an associated "// &
208 "control structure.")
215 "Parameterization of shear-driven turbulence via CVMix (various options)")
216 call get_param(param_file,
mdl,
"USE_LMD94", cs%use_LMD94, &
217 "If true, use the Large-McWilliams-Doney (JGR 1994) "//&
218 "shear mixing parameterization.", default=.false.)
219 if (cs%use_LMD94)
then
220 numbertrue=numbertrue + 1
223 call get_param(param_file,
mdl,
"USE_PP81", cs%use_PP81, &
224 "If true, use the Pacanowski and Philander (JPO 1981) "//&
225 "shear mixing parameterization.", default=.false.)
226 if (cs%use_PP81)
then
227 numbertrue = numbertrue + 1
231 if (use_jhl) numbertrue = numbertrue + 1
234 if ((numbertrue) > 1)
then
235 call mom_error(fatal,
'MOM_CVMix_shear_init: '// &
236 'Multiple shear driven internal mixing schemes selected,'//&
237 ' please disable all but one scheme to proceed.')
243 call get_param(param_file,
mdl,
"NU_ZERO", cs%Nu_Zero, &
244 "Leading coefficient in KPP shear mixing.", &
245 units=
"nondim", default=5.e-3)
246 call get_param(param_file,
mdl,
"RI_ZERO", cs%Ri_Zero, &
247 "Critical Richardson for KPP shear mixing, "// &
248 "NOTE this the internal mixing and this is "// &
249 "not for setting the boundary layer depth." &
250 ,units=
"nondim", default=0.8)
251 call get_param(param_file,
mdl,
"KPP_EXP", cs%KPP_exp, &
252 "Exponent of unitless factor of diffusivities, "// &
253 "for KPP internal shear mixing scheme." &
254 ,units=
"nondim", default=3.0)
255 call get_param(param_file,
mdl,
"SMOOTH_RI", cs%smooth_ri, &
256 "If true, vertically smooth the Richardson "// &
257 "number by applying a 1-2-1 filter once.", &
259 call cvmix_init_shear(mix_scheme=cs%Mix_Scheme, &
260 kpp_nu_zero=cs%Nu_Zero, &
261 kpp_ri_zero=cs%Ri_zero, &
267 cs%id_N2 = register_diag_field(
'ocean_model',
'N2_shear', diag%axesTi, time, &
268 'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module',
'1/s2')
269 if (cs%id_N2 > 0)
then
270 allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%N2(:,:,:) = 0.
273 cs%id_S2 = register_diag_field(
'ocean_model',
'S2_shear', diag%axesTi, time, &
274 'Square of vertical shear used by MOM_CVMix_shear module',
'1/s2')
275 if (cs%id_S2 > 0)
then
276 allocate( cs%S2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%S2(:,:,:) = 0.
279 cs%id_ri_grad = register_diag_field(
'ocean_model',
'ri_grad_shear', diag%axesTi, time, &
280 'Gradient Richarson number used by MOM_CVMix_shear module',
'nondim')
281 if (cs%id_ri_grad > 0)
then
282 allocate( cs%ri_grad( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad(:,:,:) = 1.e8
285 cs%id_ri_grad_smooth = register_diag_field(
'ocean_model',
'ri_grad_shear_smooth', &
287 'Smoothed gradient Richarson number used by MOM_CVMix_shear module',
'nondim')
288 if (cs%id_ri_grad_smooth > 0)
then
289 allocate( cs%ri_grad_smooth( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad_smooth(:,:,:) = 1.e8
292 cs%id_kd = register_diag_field(
'ocean_model',
'kd_shear_CVMix', diag%axesTi, time, &
293 'Vertical diffusivity added by MOM_CVMix_shear module',
'm2/s', conversion=us%Z2_T_to_m2_s)
294 cs%id_kv = register_diag_field(
'ocean_model',
'kv_shear_CVMix', diag%axesTi, time, &
295 'Vertical viscosity added by MOM_CVMix_shear module',
'm2/s', conversion=us%Z2_T_to_m2_s)
305 logical :: lmd94, pp81
307 default=.false., do_not_log = .true.)
309 default=.false., do_not_log = .true.)
318 if (.not.
associated(cs))
return
320 if (cs%id_N2 > 0)
deallocate(cs%N2)
321 if (cs%id_S2 > 0)
deallocate(cs%S2)
322 if (cs%id_ri_grad > 0)
deallocate(cs%ri_grad)