17 use cvmix_shear,
only : cvmix_init_shear, cvmix_coeffs_shear
19 implicit none ;
private
21 #include <MOM_memory.h>
23 public calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used, cvmix_shear_end
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"
59 subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
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
65 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
67 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kd
69 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(out) :: kv
74 integer :: i, j, k, kk, km1
76 real :: pref, du, dv, drho, dz, n2, s2, dummy
77 real,
dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d
78 real,
dimension(G%ke+1) :: ri_grad
79 real,
dimension(G%ke+1) :: kvisc
80 real,
dimension(G%ke+1) :: kdiff
81 real,
parameter :: epsln = 1.e-10
84 gorho = gv%mks_g_Earth / gv%Rho0
90 if (g%mask2dT(i,j)==0.) cycle
103 temp_1d(kk+1) = tv%T(i,j,k)
104 temp_1d(kk+2) = tv%T(i,j,km1)
105 salt_1d(kk+1) = tv%S(i,j,k)
106 salt_1d(kk+2) = tv%S(i,j,km1)
110 pref = pref + gv%H_to_Pa * h(i,j,k)
115 call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 2*g%ke, tv%EQN_OF_STATE)
121 du = (u_h(i,j,k))-(u_h(i,j,km1))
122 dv = (v_h(i,j,k))-(v_h(i,j,km1))
123 drho = (gorho * (rho_1d(kk+1) - rho_1d(kk+2)) )
124 dz = ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
126 s2 = (du*du+dv*dv)/(dz*dz)
127 ri_grad(k) = max(0.,n2)/max(s2,1.e-10)
130 if (cs%id_N2 > 0) cs%N2(i,j,k) = n2
131 if (cs%id_S2 > 0) cs%S2(i,j,k) = s2
135 ri_grad(g%ke+1) = ri_grad(g%ke)
137 if (cs%id_ri_grad > 0) cs%ri_grad(i,j,:) = ri_grad(:)
139 if (cs%smooth_ri)
then
142 if (h(i,j,k) .le. epsln) ri_grad(k) = ri_grad(k-1)
145 ri_grad(g%ke+1) = ri_grad(g%ke)
148 dummy = 0.25 * ri_grad(2)
149 ri_grad(g%ke+1) = ri_grad(g%ke)
151 ri_grad(k) = dummy + 0.5 * ri_grad(k) + 0.25 * ri_grad(k+1)
152 dummy = 0.25 * ri_grad(k)
155 if (cs%id_ri_grad_smooth > 0) cs%ri_grad_smooth(i,j,:) = ri_grad(:)
159 kvisc(k) = us%Z2_T_to_m2_s * kv(i,j,k)
160 kdiff(k) = us%Z2_T_to_m2_s * kd(i,j,k)
164 call cvmix_coeffs_shear(mdiff_out=kvisc(:), &
165 tdiff_out=kdiff(:), &
170 kv(i,j,k) = us%m2_s_to_Z2_T * kvisc(k)
171 kd(i,j,k) = us%m2_s_to_Z2_T * kdiff(k)
177 if (cs%id_kd > 0)
call post_data(cs%id_kd, kd, cs%diag)
178 if (cs%id_kv > 0)
call post_data(cs%id_kv, kv, cs%diag)
179 if (cs%id_N2 > 0)
call post_data(cs%id_N2, cs%N2, cs%diag)
180 if (cs%id_S2 > 0)
call post_data(cs%id_S2, cs%S2, cs%diag)
181 if (cs%id_ri_grad > 0)
call post_data(cs%id_ri_grad, cs%ri_grad, cs%diag)
182 if (cs%id_ri_grad_smooth > 0)
call post_data(cs%id_ri_grad_smooth ,cs%ri_grad_smooth, cs%diag)
184 end subroutine calculate_cvmix_shear
191 logical function cvmix_shear_init(Time, G, GV, US, param_file, diag, CS)
192 type(time_type),
intent(in) :: time
197 type(
diag_ctrl),
target,
intent(inout) :: diag
200 integer :: numbertrue=0
203 #include "version_variable.h"
205 if (
associated(cs))
then
206 call mom_error(warning,
"CVMix_shear_init called with an associated "// &
207 "control structure.")
214 "Parameterization of shear-driven turbulence via CVMix (various options)")
215 call get_param(param_file, mdl,
"USE_LMD94", cs%use_LMD94, &
216 "If true, use the Large-McWilliams-Doney (JGR 1994) "//&
217 "shear mixing parameterization.", default=.false.)
218 if (cs%use_LMD94)
then
219 numbertrue=numbertrue + 1
222 call get_param(param_file, mdl,
"USE_PP81", cs%use_PP81, &
223 "If true, use the Pacanowski and Philander (JPO 1981) "//&
224 "shear mixing parameterization.", default=.false.)
225 if (cs%use_PP81)
then
226 numbertrue = numbertrue + 1
229 use_jhl=kappa_shear_is_used(param_file)
230 if (use_jhl) numbertrue = numbertrue + 1
233 if ((numbertrue) > 1)
then
234 call mom_error(fatal,
'MOM_CVMix_shear_init: '// &
235 'Multiple shear driven internal mixing schemes selected,'//&
236 ' please disable all but one scheme to proceed.')
238 cvmix_shear_init=(cs%use_PP81.or.cs%use_LMD94)
241 if (.not. cvmix_shear_init)
return
242 call get_param(param_file, mdl,
"NU_ZERO", cs%Nu_Zero, &
243 "Leading coefficient in KPP shear mixing.", &
244 units=
"nondim", default=5.e-3)
245 call get_param(param_file, mdl,
"RI_ZERO", cs%Ri_Zero, &
246 "Critical Richardson for KPP shear mixing, "// &
247 "NOTE this the internal mixing and this is "// &
248 "not for setting the boundary layer depth." &
249 ,units=
"nondim", default=0.8)
250 call get_param(param_file, mdl,
"KPP_EXP", cs%KPP_exp, &
251 "Exponent of unitless factor of diffusivities, "// &
252 "for KPP internal shear mixing scheme." &
253 ,units=
"nondim", default=3.0)
254 call get_param(param_file, mdl,
"SMOOTH_RI", cs%smooth_ri, &
255 "If true, vertically smooth the Richardson "// &
256 "number by applying a 1-2-1 filter once.", &
258 call cvmix_init_shear(mix_scheme=cs%Mix_Scheme, &
259 kpp_nu_zero=cs%Nu_Zero, &
260 kpp_ri_zero=cs%Ri_zero, &
266 cs%id_N2 = register_diag_field(
'ocean_model',
'N2_shear', diag%axesTi, time, &
267 'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module',
'1/s2')
268 if (cs%id_N2 > 0)
then
269 allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%N2(:,:,:) = 0.
272 cs%id_S2 = register_diag_field(
'ocean_model',
'S2_shear', diag%axesTi, time, &
273 'Square of vertical shear used by MOM_CVMix_shear module',
'1/s2')
274 if (cs%id_S2 > 0)
then
275 allocate( cs%S2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%S2(:,:,:) = 0.
278 cs%id_ri_grad = register_diag_field(
'ocean_model',
'ri_grad_shear', diag%axesTi, time, &
279 'Gradient Richarson number used by MOM_CVMix_shear module',
'nondim')
280 if (cs%id_ri_grad > 0)
then
281 allocate( cs%ri_grad( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad(:,:,:) = 1.e8
284 cs%id_ri_grad_smooth = register_diag_field(
'ocean_model',
'ri_grad_shear_smooth', &
286 'Smoothed gradient Richarson number used by MOM_CVMix_shear module',
'nondim')
287 if (cs%id_ri_grad_smooth > 0)
then
288 allocate( cs%ri_grad_smooth( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad_smooth(:,:,:) = 1.e8
291 cs%id_kd = register_diag_field(
'ocean_model',
'kd_shear_CVMix', diag%axesTi, time, &
292 'Vertical diffusivity added by MOM_CVMix_shear module',
'm2/s', conversion=us%Z2_T_to_m2_s)
293 cs%id_kv = register_diag_field(
'ocean_model',
'kv_shear_CVMix', diag%axesTi, time, &
294 'Vertical viscosity added by MOM_CVMix_shear module',
'm2/s', conversion=us%Z2_T_to_m2_s)
296 end function cvmix_shear_init
301 logical function cvmix_shear_is_used(param_file)
304 logical :: lmd94, pp81
305 call get_param(param_file, mdl,
"USE_LMD94", lmd94, &
306 default=.false., do_not_log = .true.)
307 call get_param(param_file, mdl,
"Use_PP81", pp81, &
308 default=.false., do_not_log = .true.)
309 cvmix_shear_is_used = (lmd94 .or. pp81)
310 end function cvmix_shear_is_used
313 subroutine cvmix_shear_end(CS)
317 if (.not.
associated(cs))
return
319 if (cs%id_N2 > 0)
deallocate(cs%N2)
320 if (cs%id_S2 > 0)
deallocate(cs%S2)
321 if (cs%id_ri_grad > 0)
deallocate(cs%ri_grad)
324 end subroutine cvmix_shear_end