21 implicit none ;
private
23 #include <MOM_memory.h>
27 logical :: use_variable_mixing
28 logical :: resoln_scaled_kh
30 logical :: resoln_scaled_khth
32 logical :: depth_scaled_khth
34 logical :: resoln_scaled_khtr
36 logical :: interpolate_res_fn
41 logical :: use_stored_slopes
42 logical :: resoln_use_ebt
44 logical :: khth_use_ebt_struct
46 logical :: calculate_cg1
49 logical :: calculate_rd_dx
51 logical :: calculate_res_fns
53 logical :: calculate_depth_fns
55 logical :: calculate_eady_growth_rate
57 real,
dimension(:,:),
pointer :: &
71 depth_fn_u => null(), &
73 depth_fn_v => null(), &
75 beta_dx2_h => null(), &
77 beta_dx2_q => null(), &
79 beta_dx2_u => null(), &
81 beta_dx2_v => null(), &
93 real,
dimension(:,:,:),
pointer :: &
100 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: &
103 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: &
106 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
109 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
113 logical :: use_visbeck
114 integer :: varmix_ktop
115 real :: visbeck_l_scale
116 real :: res_coef_khth
119 real :: res_coef_visc
122 real :: depth_scaled_khth_h0
123 real :: depth_scaled_khth_exp
125 integer :: res_fn_power_khth
128 integer :: res_fn_power_visc
131 real :: visbeck_s_max
134 logical :: use_qg_leith_gm
135 logical :: use_beta_in_qg_leith
140 integer :: id_sn_u=-1, id_sn_v=-1, id_l2u=-1, id_l2v=-1, id_res_fn = -1
141 integer :: id_n2_u=-1, id_n2_v=-1, id_s2_u=-1, id_s2_v=-1
142 integer :: id_rd_dx=-1, id_kh_u_qg = -1, id_kh_v_qg = -1
148 type(group_pass_type) :: pass_cg1
163 integer :: is, ie, js, je, isq, ieq, jsq, jeq
167 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
168 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
170 if (.not.
associated(cs))
call mom_error(fatal,
"calc_depth_function:"// &
171 "Module must be initialized before it is used.")
172 if (.not. cs%calculate_depth_fns)
return
173 if (.not.
associated(cs%Depth_fn_u))
call mom_error(fatal, &
174 "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
175 if (.not.
associated(cs%Depth_fn_v))
call mom_error(fatal, &
176 "calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.")
178 h0 = cs%depth_scaled_khth_h0
179 expo = cs%depth_scaled_khth_exp
181 do j=js,je ;
do i=is-1,ieq
182 cs%Depth_fn_u(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))/h0))**expo
185 do j=js-1,jeq ;
do i=is,ie
186 cs%Depth_fn_v(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))/h0))**expo
194 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
208 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
210 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
211 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
213 if (.not.
associated(cs))
call mom_error(fatal,
"calc_resoln_function:"// &
214 "Module must be initialized before it is used.")
215 if (cs%calculate_cg1)
then
216 if (.not.
associated(cs%cg1))
call mom_error(fatal, &
217 "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
218 if (cs%khth_use_ebt_struct)
then
219 if (.not.
associated(cs%ebt_struct))
call mom_error(fatal, &
220 "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
221 if (cs%Resoln_use_ebt)
then
223 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct)
226 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct, &
228 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
230 call pass_var(cs%ebt_struct, g%Domain)
232 call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
236 call do_group_pass(cs%pass_cg1, g%Domain)
241 if (cs%calculate_rd_dx)
then
242 if (.not.
associated(cs%Rd_dx_h))
call mom_error(fatal, &
243 "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
245 do j=js-1,je+1 ;
do i=is-1,ie+1
246 cs%Rd_dx_h(i,j) = cs%cg1(i,j) / &
247 (sqrt(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))
250 if (cs%id_Rd_dx > 0)
call post_data(cs%id_Rd_dx, cs%Rd_dx_h, cs%diag)
254 if (.not. cs%calculate_res_fns)
return
256 if (.not.
associated(cs%Res_fn_h))
call mom_error(fatal, &
257 "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
258 if (.not.
associated(cs%Res_fn_q))
call mom_error(fatal, &
259 "calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
260 if (.not.
associated(cs%Res_fn_u))
call mom_error(fatal, &
261 "calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
262 if (.not.
associated(cs%Res_fn_v))
call mom_error(fatal, &
263 "calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
264 if (.not.
associated(cs%f2_dx2_h))
call mom_error(fatal, &
265 "calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
266 if (.not.
associated(cs%f2_dx2_q))
call mom_error(fatal, &
267 "calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
268 if (.not.
associated(cs%f2_dx2_u))
call mom_error(fatal, &
269 "calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
270 if (.not.
associated(cs%f2_dx2_v))
call mom_error(fatal, &
271 "calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
272 if (.not.
associated(cs%beta_dx2_h))
call mom_error(fatal, &
273 "calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
274 if (.not.
associated(cs%beta_dx2_q))
call mom_error(fatal, &
275 "calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
276 if (.not.
associated(cs%beta_dx2_u))
call mom_error(fatal, &
277 "calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
278 if (.not.
associated(cs%beta_dx2_v))
call mom_error(fatal, &
279 "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")
286 if (cs%Res_fn_power_visc >= 100)
then
288 do j=js-1,je+1 ;
do i=is-1,ie+1
289 dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
290 if ((cs%Res_coef_visc * cs%cg1(i,j))**2 > dx_term)
then
291 cs%Res_fn_h(i,j) = 0.0
293 cs%Res_fn_h(i,j) = 1.0
297 do j=js-1,jeq ;
do i=is-1,ieq
298 cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
299 dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
300 if ((cs%Res_coef_visc * cg1_q)**2 > dx_term)
then
301 cs%Res_fn_q(i,j) = 0.0
303 cs%Res_fn_q(i,j) = 1.0
306 elseif (cs%Res_fn_power_visc == 2)
then
308 do j=js-1,je+1 ;
do i=is-1,ie+1
309 dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
310 cs%Res_fn_h(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**2)
313 do j=js-1,jeq ;
do i=is-1,ieq
314 cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
315 dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
316 cs%Res_fn_q(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cg1_q)**2)
318 elseif (mod(cs%Res_fn_power_visc, 2) == 0)
then
319 power_2 = cs%Res_fn_power_visc / 2
321 do j=js-1,je+1 ;
do i=is-1,ie+1
322 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
323 cs%Res_fn_h(i,j) = dx_term / &
324 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
327 do j=js-1,jeq ;
do i=is-1,ieq
328 cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
329 dx_term = (us%L_T_to_m_s**2*(cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)))**power_2
330 cs%Res_fn_q(i,j) = dx_term / &
331 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
335 do j=js-1,je+1 ;
do i=is-1,ie+1
336 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_h(i,j) + &
337 cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**cs%Res_fn_power_visc
338 cs%Res_fn_h(i,j) = dx_term / &
339 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
342 do j=js-1,jeq ;
do i=is-1,ieq
343 cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
344 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_q(i,j) + &
345 cg1_q * cs%beta_dx2_q(i,j)))**cs%Res_fn_power_visc
346 cs%Res_fn_q(i,j) = dx_term / &
347 (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
351 if (cs%interpolate_Res_fn)
then
352 do j=js,je ;
do i=is-1,ieq
353 cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
355 do j=js-1,jeq ;
do i=is,ie
356 cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
359 if (cs%Res_fn_power_khth >= 100)
then
361 do j=js,je ;
do i=is-1,ieq
362 cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
363 dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
364 if ((cs%Res_coef_khth * cg1_u)**2 > dx_term)
then
365 cs%Res_fn_u(i,j) = 0.0
367 cs%Res_fn_u(i,j) = 1.0
371 do j=js-1,jeq ;
do i=is,ie
372 cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
373 dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
374 if ((cs%Res_coef_khth * cg1_v)**2 > dx_term)
then
375 cs%Res_fn_v(i,j) = 0.0
377 cs%Res_fn_v(i,j) = 1.0
380 elseif (cs%Res_fn_power_khth == 2)
then
382 do j=js,je ;
do i=is-1,ieq
383 cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
384 dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
385 cs%Res_fn_u(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_u)**2)
388 do j=js-1,jeq ;
do i=is,ie
389 cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
390 dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
391 cs%Res_fn_v(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_v)**2)
393 elseif (mod(cs%Res_fn_power_khth, 2) == 0)
then
394 power_2 = cs%Res_fn_power_khth / 2
396 do j=js,je ;
do i=is-1,ieq
397 cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
398 dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)))**power_2
399 cs%Res_fn_u(i,j) = dx_term / &
400 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
403 do j=js-1,jeq ;
do i=is,ie
404 cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
405 dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)))**power_2
406 cs%Res_fn_v(i,j) = dx_term / &
407 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
411 do j=js,je ;
do i=is-1,ieq
412 cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
413 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_u(i,j) + &
414 cg1_u * cs%beta_dx2_u(i,j)))**cs%Res_fn_power_khth
415 cs%Res_fn_u(i,j) = dx_term / &
416 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
419 do j=js-1,jeq ;
do i=is,ie
420 cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
421 dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_v(i,j) + &
422 cg1_v * cs%beta_dx2_v(i,j)))**cs%Res_fn_power_khth
423 cs%Res_fn_v(i,j) = dx_term / &
424 (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
431 if (cs%id_Res_fn > 0)
call post_data(cs%id_Res_fn, cs%Res_fn_h, cs%diag)
442 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
444 real,
intent(in) :: dt
447 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
449 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: n2_u
450 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: n2_v
452 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//&
453 "Module must be initialized before it is used.")
455 if (cs%calculate_Eady_growth_rate)
then
456 call find_eta(h, tv, g, gv, us, e, halo_size=2)
457 if (cs%use_stored_slopes)
then
459 cs%slope_x, cs%slope_y, n2_u, n2_v, 1)
469 if (cs%id_SN_u > 0)
call post_data(cs%id_SN_u, cs%SN_u, cs%diag)
470 if (cs%id_SN_v > 0)
call post_data(cs%id_SN_v, cs%SN_v, cs%diag)
471 if (cs%id_L2u > 0)
call post_data(cs%id_L2u, cs%L2u, cs%diag)
472 if (cs%id_L2v > 0)
call post_data(cs%id_L2v, cs%L2v, cs%diag)
475 if (cs%id_N2_u > 0)
call post_data(cs%id_N2_u, cs%N2_u, cs%diag)
476 if (cs%id_N2_v > 0)
call post_data(cs%id_N2_v, cs%N2_v, cs%diag)
485 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
486 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: slope_x
487 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: N2_u
489 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(in) :: slope_y
490 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(in) :: N2_v
500 integer :: is, ie, js, je, nz
501 integer :: i, j, k, kb_max
502 real :: S2max, wNE, wSE, wSW, wNW
503 real :: H_u(SZIB_(G)), H_v(SZI_(G))
504 real :: S2_u(SZIB_(G), SZJ_(G))
505 real :: S2_v(SZI_(G), SZJB_(G))
507 if (.not.
associated(cs))
call mom_error(fatal,
"calc_slope_function:"// &
508 "Module must be initialized before it is used.")
509 if (.not. cs%calculate_Eady_growth_rate)
return
510 if (.not.
associated(cs%SN_u))
call mom_error(fatal,
"calc_slope_function:"// &
511 "%SN_u is not associated with use_variable_mixing.")
512 if (.not.
associated(cs%SN_v))
call mom_error(fatal,
"calc_slope_function:"// &
513 "%SN_v is not associated with use_variable_mixing.")
515 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
517 s2max = cs%Visbeck_S_max**2
520 do j=js-1,je+1 ;
do i=is-1,ie+1
532 cs%SN_u(i,j) = 0. ; h_u(i) = 0. ; s2_u(i,j) = 0.
534 do k=2,nz ;
do i=is-1,ie
535 hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
536 hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
537 h_geom = sqrt( hdn * hup )
540 wse = g%mask2dCv(i+1,j-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) )
541 wnw = g%mask2dCv(i ,j ) * ( (h(i ,j,k)*h(i ,j+1,k)) * (h(i ,j,k-1)*h(i ,j+1,k-1)) )
542 wne = g%mask2dCv(i+1,j ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) )
543 wsw = g%mask2dCv(i ,j-1) * ( (h(i ,j,k)*h(i ,j-1,k)) * (h(i ,j,k-1)*h(i ,j-1,k-1)) )
544 s2 = slope_x(i,j,k)**2 + &
545 ((wnw*slope_y(i,j,k)**2 + wse*slope_y(i+1,j-1,k)**2) + &
546 (wne*slope_y(i+1,j,k)**2 + wsw*slope_y(i,j-1,k)**2) ) / &
547 ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
548 if (s2max>0.) s2 = s2 * s2max / (s2 + s2max)
550 n2 = max(0., n2_u(i,j,k))
551 cs%SN_u(i,j) = cs%SN_u(i,j) + sqrt( s2*n2 )*h_geom
552 s2_u(i,j) = s2_u(i,j) + s2*h_geom
553 h_u(i) = h_u(i) + h_geom
557 cs%SN_u(i,j) = g%mask2dCu(i,j) * cs%SN_u(i,j) / h_u(i)
558 s2_u(i,j) = g%mask2dCu(i,j) * s2_u(i,j) / h_u(i)
568 cs%SN_v(i,j) = 0. ; h_v(i) = 0. ; s2_v(i,j) = 0.
570 do k=2,nz ;
do i=is,ie
571 hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
572 hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
573 h_geom = sqrt( hdn * hup )
576 wse = g%mask2dCu(i,j) * ( (h(i,j ,k)*h(i+1,j ,k)) * (h(i,j ,k-1)*h(i+1,j ,k-1)) )
577 wnw = g%mask2dCu(i-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) )
578 wne = g%mask2dCu(i,j+1) * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) )
579 wsw = g%mask2dCu(i-1,j) * ( (h(i,j ,k)*h(i-1,j ,k)) * (h(i,j ,k-1)*h(i-1,j ,k-1)) )
580 s2 = slope_y(i,j,k)**2 + &
581 ((wse*slope_x(i,j,k)**2 + wnw*slope_x(i-1,j+1,k)**2) + &
582 (wne*slope_x(i,j+1,k)**2 + wsw*slope_x(i-1,j,k)**2) ) / &
583 ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
584 if (s2max>0.) s2 = s2 * s2max / (s2 + s2max)
586 n2 = max(0., n2_v(i,j,k))
587 cs%SN_v(i,j) = cs%SN_v(i,j) + sqrt( s2*n2 )*h_geom
588 s2_v(i,j) = s2_v(i,j) + s2*h_geom
589 h_v(i) = h_v(i) + h_geom
593 cs%SN_v(i,j) = g%mask2dCv(i,j) * cs%SN_v(i,j) / h_v(i)
594 s2_v(i,j) = g%mask2dCv(i,j) * s2_v(i,j) / h_v(i)
603 if (cs%id_S2_u > 0)
call post_data(cs%id_S2_u, s2_u, cs%diag)
604 if (cs%id_S2_v > 0)
call post_data(cs%id_S2_v, s2_v, cs%diag)
608 call uvchksum(
"calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, g%HI, haloshift=1)
609 call uvchksum(
"calc_Visbeck_coeffs N2_u, N2_v", n2_u, n2_v, g%HI, scale=us%s_to_T**2)
610 call uvchksum(
"calc_Visbeck_coeffs SN_[uv]", cs%SN_u, cs%SN_v, g%HI, scale=us%s_to_T)
619 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
623 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
624 logical,
intent(in) :: calculate_slopes
627 real :: E_x(SZIB_(G), SZJ_(G))
628 real :: E_y(SZI_(G), SZJB_(G))
637 integer :: is, ie, js, je, nz
638 integer :: i, j, k, kb_max
639 real :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(G))
640 real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(G))
642 if (.not.
associated(cs))
call mom_error(fatal,
"calc_slope_function:"// &
643 "Module must be initialized before it is used.")
644 if (.not. cs%calculate_Eady_growth_rate)
return
645 if (.not.
associated(cs%SN_u))
call mom_error(fatal,
"calc_slope_function:"// &
646 "%SN_u is not associated with use_variable_mixing.")
647 if (.not.
associated(cs%SN_v))
call mom_error(fatal,
"calc_slope_function:"// &
648 "%SN_v is not associated with use_variable_mixing.")
650 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
652 one_meter = 1.0 * gv%m_to_H
653 h_neglect = gv%H_subroundoff
654 h_cutoff = real(2*nz) * (gv%Angstrom_H + h_neglect)
661 do k=nz,cs%VarMix_Ktop,-1
663 if (calculate_slopes)
then
665 do j=js-1,je+1 ;
do i=is-1,ie
666 e_x(i,j) = us%Z_to_L*(e(i+1,j,k)-e(i,j,k))*g%IdxCu(i,j)
668 if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
670 do j=js-1,je ;
do i=is-1,ie+1
671 e_y(i,j) = us%Z_to_L*(e(i,j+1,k)-e(i,j,k))*g%IdyCv(i,j)
673 if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
676 do j=js-1,je+1 ;
do i=is-1,ie
677 e_x(i,j) = cs%slope_x(i,j,k)
678 if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
680 do j=js-1,je ;
do i=is-1,ie+1
681 e_y(i,j) = cs%slope_y(i,j,k)
682 if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
687 do j=js,je ;
do i=is-1,ie
688 s2 = ( e_x(i,j)**2 + 0.25*( &
689 (e_y(i,j)**2+e_y(i+1,j-1)**2)+(e_y(i+1,j)**2+e_y(i,j-1)**2) ) )
690 hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
691 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)
692 h_geom = sqrt(hdn*hup)
693 n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
694 if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < h_cutoff) &
696 s2n2_u_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
698 do j=js-1,je ;
do i=is,ie
699 s2 = ( e_y(i,j)**2 + 0.25*( &
700 (e_x(i,j)**2+e_x(i-1,j+1)**2)+(e_x(i,j+1)**2+e_x(i-1,j)**2) ) )
701 hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
702 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)
703 h_geom = sqrt(hdn*hup)
704 n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
705 if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < h_cutoff) &
707 s2n2_v_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
713 do i=is-1,ie ; cs%SN_u(i,j) = 0.0 ;
enddo
714 do k=nz,cs%VarMix_Ktop,-1 ;
do i=is-1,ie
715 cs%SN_u(i,j) = cs%SN_u(i,j) + s2n2_u_local(i,j,k)
721 if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z )
then
722 cs%SN_u(i,j) = g%mask2dCu(i,j) * sqrt( cs%SN_u(i,j) / &
723 (max(g%bathyT(i,j), g%bathyT(i+1,j))) )
731 do i=is,ie ; cs%SN_v(i,j) = 0.0 ;
enddo
732 do k=nz,cs%VarMix_Ktop,-1 ;
do i=is,ie
733 cs%SN_v(i,j) = cs%SN_v(i,j) + s2n2_v_local(i,j,k)
738 if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z )
then
739 cs%SN_v(i,j) = g%mask2dCv(i,j) * sqrt( cs%SN_v(i,j) / &
740 (max(g%bathyT(i,j), g%bathyT(i,j+1))) )
757 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
758 integer,
intent(in) :: k
759 real,
dimension(SZIB_(G),SZJ_(G)),
intent(in) :: div_xx_dx
761 real,
dimension(SZI_(G),SZJB_(G)),
intent(in) :: div_xx_dy
763 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vort_xy_dx
765 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: vort_xy_dy
777 real,
dimension(SZI_(G),SZJB_(G)) :: &
784 real,
dimension(SZIB_(G),SZJ_(G)) :: &
790 real :: h_at_slope_above
791 real :: h_at_slope_below
794 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq,nz
797 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
798 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
801 inv_pi3 = 1.0/((4.0*atan(1.0))**3)
806 if ((k > 1) .and. (k < nz))
then
812 do j=js-2,jeq+2 ;
do i=is-2,ieq+1
813 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) ) / &
814 ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
815 + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + gv%H_subroundoff )
816 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) ) / &
817 ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
818 + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + gv%H_subroundoff )
819 ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
820 dslopex_dz(i,j) = 2. * ( cs%slope_x(i,j,k) - cs%slope_x(i,j,k+1) ) * ih
821 h_at_u(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
825 do j=js-2,jeq+1 ;
do i=is-2,ieq+2
826 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) ) / &
827 ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
828 + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + gv%H_subroundoff )
829 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) ) / &
830 ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
831 + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + gv%H_subroundoff )
832 ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
833 dslopey_dz(i,j) = 2. * ( cs%slope_y(i,j,k) - cs%slope_y(i,j,k+1) ) * ih
834 h_at_v(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
838 do j=js-2,jeq+1 ;
do i=is-1,ieq+1
839 f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
840 vort_xy_dx(i,j) = vort_xy_dx(i,j) - f * us%L_to_Z * &
841 ( ( h_at_u(i,j) * dslopex_dz(i,j) + h_at_u(i-1,j+1) * dslopex_dz(i-1,j+1) ) &
842 + ( h_at_u(i-1,j) * dslopex_dz(i-1,j) + h_at_u(i,j+1) * dslopex_dz(i,j+1) ) ) / &
843 ( ( 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)
847 do j=js-1,jeq+1 ;
do i=is-2,ieq+1
848 f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
850 vort_xy_dy(i,j) = vort_xy_dx(i,j) - f * us%L_to_Z * &
851 ( ( h_at_v(i,j) * dslopey_dz(i,j) + h_at_v(i+1,j-1) * dslopey_dz(i+1,j-1) ) &
852 + ( h_at_v(i,j-1) * dslopey_dz(i,j-1) + h_at_v(i+1,j) * dslopey_dz(i+1,j) ) ) / &
853 ( ( 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)
860 if (cs%use_QG_Leith_GM)
then
862 do j=js,je ;
do i=is-1,ieq
866 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) &
867 + vort_xy_dx(i,j-1) + vort_xy_dx(i+1,j-1)))**2)
868 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) &
869 + div_xx_dy(i,j-1) + div_xx_dy(i+1,j-1)))**2)
870 if (cs%use_beta_in_QG_Leith)
then
871 beta_u(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i+1,j))**2) + &
872 (0.5*(g%dF_dy(i,j)+g%dF_dy(i+1,j))**2) )
873 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)) * &
874 cs%Laplac3_const_u(i,j) * inv_pi3
876 cs%KH_u_QG(i,j,k) = (grad_vort_mag_u(i,j) + grad_div_mag_u(i,j)) * &
877 cs%Laplac3_const_u(i,j) * inv_pi3
881 do j=js-1,jeq ;
do i=is,ie
883 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) &
884 + vort_xy_dy(i,j+1) + vort_xy_dy(i-1,j+1)))**2)
885 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) &
886 + div_xx_dx(i,j+1) + div_xx_dx(i-1,j+1)))**2)
887 if (cs%use_beta_in_QG_Leith)
then
888 beta_v(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i,j+1))**2) + &
889 (0.5*(g%dF_dy(i,j)+g%dF_dy(i,j+1))**2) )
890 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)) * &
891 cs%Laplac3_const_v(i,j) * inv_pi3
893 cs%KH_v_QG(i,j,k) = (grad_vort_mag_v(i,j) + grad_div_mag_v(i,j)) * &
894 cs%Laplac3_const_v(i,j) * inv_pi3
900 if (cs%id_KH_v_QG > 0)
call post_data(cs%id_KH_v_QG, cs%KH_v_QG, cs%diag)
901 if (cs%id_KH_u_QG > 0)
call post_data(cs%id_KH_u_QG, cs%KH_u_QG, cs%diag)
908 subroutine varmix_init(Time, G, GV, US, param_file, diag, CS)
909 type(time_type),
intent(in) :: time
914 type(
diag_ctrl),
target,
intent(inout) :: diag
917 real :: khtr_slope_cff, khth_slope_cff, oneortwo
918 real :: n2_filter_depth
920 real :: khtr_passivity_coeff
921 real :: absurdly_small_freq
923 logical :: gill_equatorial_ld, use_fgnv_streamfn, use_meke, in_use
924 real :: mle_front_length
925 real :: leith_lap_const
926 real :: grid_sp_u2, grid_sp_v2
927 real :: grid_sp_u3, grid_sp_v3
929 #include "version_variable.h"
930 character(len=40) :: mdl =
"MOM_lateral_mixing_coeffs"
931 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
932 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
933 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
934 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
935 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
936 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
938 if (
associated(cs))
then
939 call mom_error(warning,
"VarMix_init called with an associated "// &
940 "control structure.")
947 cs%calculate_cg1 = .false.
948 cs%calculate_Rd_dx = .false.
949 cs%calculate_res_fns = .false.
950 cs%calculate_Eady_growth_rate = .false.
951 cs%calculate_depth_fns = .false.
954 call get_param(param_file, mdl,
"USE_VARIABLE_MIXING", cs%use_variable_mixing,&
955 "If true, the variable mixing code will be called. This "//&
956 "allows diagnostics to be created even if the scheme is "//&
957 "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
958 "this is set to true regardless of what is in the "//&
959 "parameter file.", default=.false.)
960 call get_param(param_file, mdl,
"USE_VISBECK", cs%use_Visbeck,&
961 "If true, use the Visbeck et al. (1997) formulation for \n"//&
962 "thickness diffusivity.", default=.false.)
963 call get_param(param_file, mdl,
"RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
964 "If true, the Laplacian lateral viscosity is scaled away "//&
965 "when the first baroclinic deformation radius is well "//&
966 "resolved.", default=.false.)
967 call get_param(param_file, mdl,
"DEPTH_SCALED_KHTH", cs%Depth_scaled_KhTh, &
968 "If true, KHTH is scaled away when the depth is shallower"//&
969 "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//&
970 "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
971 "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
973 call get_param(param_file, mdl,
"RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
974 "If true, the interface depth diffusivity is scaled away "//&
975 "when the first baroclinic deformation radius is well "//&
976 "resolved.", default=.false.)
977 call get_param(param_file, mdl,
"RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
978 "If true, the epipycnal tracer diffusivity is scaled "//&
979 "away when the first baroclinic deformation radius is "//&
980 "well resolved.", default=.false.)
981 call get_param(param_file, mdl,
"RESOLN_USE_EBT", cs%Resoln_use_ebt, &
982 "If true, uses the equivalent barotropic wave speed instead "//&
983 "of first baroclinic wave for calculating the resolution fn.",&
985 call get_param(param_file, mdl,
"KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
986 "If true, uses the equivalent barotropic structure "//&
987 "as the vertical structure of thickness diffusivity.",&
989 call get_param(param_file, mdl,
"KHTH_SLOPE_CFF", khth_slope_cff, &
990 "The nondimensional coefficient in the Visbeck formula "//&
991 "for the interface depth diffusivity", units=
"nondim", &
993 call get_param(param_file, mdl,
"KHTR_SLOPE_CFF", khtr_slope_cff, &
994 "The nondimensional coefficient in the Visbeck formula "//&
995 "for the epipycnal tracer diffusivity", units=
"nondim", &
997 call get_param(param_file, mdl,
"USE_STORED_SLOPES", cs%use_stored_slopes,&
998 "If true, the isopycnal slopes are calculated once and "//&
999 "stored for re-use. This uses more memory but avoids calling "//&
1000 "the equation of state more times than should be necessary.", &
1002 call get_param(param_file, mdl,
"VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1003 "A miniscule frequency that is used to avoid division by 0. The default "//&
1004 "value is roughly (pi / (the age of the universe)).", &
1005 default=1.0e-17, units=
"s-1", scale=us%T_to_s)
1006 call get_param(param_file, mdl,
"KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
1007 default=.false., do_not_log=.true.)
1008 cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn
1009 call get_param(param_file, mdl,
"USE_MEKE", use_meke, &
1010 default=.false., do_not_log=.true.)
1011 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
1012 cs%calculate_Eady_growth_rate = cs%calculate_Eady_growth_rate .or. use_meke
1013 call get_param(param_file, mdl,
"KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
1014 default=0., do_not_log=.true.)
1015 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
1016 call get_param(param_file, mdl,
"MLE_FRONT_LENGTH", mle_front_length, &
1017 default=0., do_not_log=.true.)
1018 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
1020 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false., do_not_log=.true.)
1023 if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct)
then
1025 call get_param(param_file, mdl,
"RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
1026 "The depth below which N2 is monotonized to avoid stratification "//&
1027 "artifacts from altering the equivalent barotropic mode structure.",&
1028 units=
"m", default=2000., scale=us%m_to_Z)
1029 allocate(cs%ebt_struct(isd:ied,jsd:jed,g%ke)) ; cs%ebt_struct(:,:,:) = 0.0
1032 if (khtr_slope_cff>0. .or. khth_slope_cff>0.)
then
1033 cs%calculate_Eady_growth_rate = .true.
1034 call get_param(param_file, mdl,
"VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
1035 "If non-zero, is an upper bound on slopes used in the "//&
1036 "Visbeck formula for diffusivity. This does not affect the "//&
1037 "isopycnal slope calculation used within thickness diffusion.", &
1038 units=
"nondim", default=0.0)
1041 if (cs%use_stored_slopes)
then
1043 allocate(cs%slope_x(isdb:iedb,jsd:jed,g%ke+1)) ; cs%slope_x(:,:,:) = 0.0
1044 allocate(cs%slope_y(isd:ied,jsdb:jedb,g%ke+1)) ; cs%slope_y(:,:,:) = 0.0
1045 allocate(cs%N2_u(isdb:iedb,jsd:jed,g%ke+1)) ; cs%N2_u(:,:,:) = 0.0
1046 allocate(cs%N2_v(isd:ied,jsdb:jedb,g%ke+1)) ; cs%N2_v(:,:,:) = 0.0
1047 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_smooth, &
1048 "A diapycnal diffusivity that is used to interpolate "//&
1049 "more sensible values of T & S into thin layers.", &
1050 units=
"m2 s-1", default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
1053 if (cs%calculate_Eady_growth_rate)
then
1055 allocate(cs%SN_u(isdb:iedb,jsd:jed)) ; cs%SN_u(:,:) = 0.0
1056 allocate(cs%SN_v(isd:ied,jsdb:jedb)) ; cs%SN_v(:,:) = 0.0
1057 cs%id_SN_u = register_diag_field(
'ocean_model',
'SN_u', diag%axesCu1, time, &
1058 'Inverse eddy time-scale, S*N, at u-points',
's-1', conversion=us%s_to_T)
1059 cs%id_SN_v = register_diag_field(
'ocean_model',
'SN_v', diag%axesCv1, time, &
1060 'Inverse eddy time-scale, S*N, at v-points',
's-1', conversion=us%s_to_T)
1061 call get_param(param_file, mdl,
"VARMIX_KTOP", cs%VarMix_Ktop, &
1062 "The layer number at which to start vertical integration "//&
1063 "of S*N for purposes of finding the Eady growth rate.", &
1064 units=
"nondim", default=2)
1067 if (khtr_slope_cff>0. .or. khth_slope_cff>0.)
then
1069 call get_param(param_file, mdl,
"VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1070 "The fixed length scale in the Visbeck formula.", units=
"m", &
1072 allocate(cs%L2u(isdb:iedb,jsd:jed)) ; cs%L2u(:,:) = 0.0
1073 allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = 0.0
1074 if (cs%Visbeck_L_scale<0)
then
1075 do j=js,je ;
do i=is-1,ieq
1076 cs%L2u(i,j) = cs%Visbeck_L_scale**2 * g%areaCu(i,j)
1078 do j=js-1,jeq ;
do i=is,ie
1079 cs%L2v(i,j) = cs%Visbeck_L_scale**2 * g%areaCv(i,j)
1082 cs%L2u(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
1083 cs%L2v(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
1086 cs%id_L2u = register_diag_field(
'ocean_model',
'L2u', diag%axesCu1, time, &
1087 'Length scale squared for mixing coefficient, at u-points', &
1088 'm2', conversion=us%L_to_m**2)
1089 cs%id_L2v = register_diag_field(
'ocean_model',
'L2v', diag%axesCv1, time, &
1090 'Length scale squared for mixing coefficient, at v-points', &
1091 'm2', conversion=us%L_to_m**2)
1094 if (cs%use_stored_slopes)
then
1095 cs%id_N2_u = register_diag_field(
'ocean_model',
'N2_u', diag%axesCui, time, &
1096 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.',
's-2')
1097 cs%id_N2_v = register_diag_field(
'ocean_model',
'N2_v', diag%axesCvi, time, &
1098 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.',
's-2')
1100 cs%id_S2_u = register_diag_field(
'ocean_model',
'S2_u', diag%axesCu1, time, &
1101 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.',
's-2')
1102 cs%id_S2_v = register_diag_field(
'ocean_model',
'S2_v', diag%axesCv1, time, &
1103 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.',
's-2')
1107 if (cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. cs%Resoln_scaled_KhTr)
then
1108 cs%calculate_Rd_dx = .true.
1109 cs%calculate_res_fns = .true.
1110 allocate(cs%Res_fn_h(isd:ied,jsd:jed)) ; cs%Res_fn_h(:,:) = 0.0
1111 allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb)) ; cs%Res_fn_q(:,:) = 0.0
1112 allocate(cs%Res_fn_u(isdb:iedb,jsd:jed)) ; cs%Res_fn_u(:,:) = 0.0
1113 allocate(cs%Res_fn_v(isd:ied,jsdb:jedb)) ; cs%Res_fn_v(:,:) = 0.0
1114 allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%beta_dx2_q(:,:) = 0.0
1115 allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed)) ; cs%beta_dx2_u(:,:) = 0.0
1116 allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb)) ; cs%beta_dx2_v(:,:) = 0.0
1117 allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%f2_dx2_q(:,:) = 0.0
1118 allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed)) ; cs%f2_dx2_u(:,:) = 0.0
1119 allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb)) ; cs%f2_dx2_v(:,:) = 0.0
1121 cs%id_Res_fn = register_diag_field(
'ocean_model',
'Res_fn', diag%axesT1, time, &
1122 'Resolution function for scaling diffusivities',
'nondim')
1124 call get_param(param_file, mdl,
"KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1125 "A coefficient that determines how KhTh is scaled away if "//&
1126 "RESOLN_SCALED_... is true, as "//&
1127 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1128 units=
"nondim", default=1.0)
1129 call get_param(param_file, mdl,
"KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1130 "The power of dx/Ld in the Kh resolution function. Any "//&
1131 "positive integer may be used, although even integers "//&
1132 "are more efficient to calculate. Setting this greater "//&
1133 "than 100 results in a step-function being used.", &
1134 units=
"nondim", default=2)
1135 call get_param(param_file, mdl,
"VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1136 "A coefficient that determines how Kh is scaled away if "//&
1137 "RESOLN_SCALED_... is true, as "//&
1138 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1139 "This function affects lateral viscosity, Kh, and not KhTh.", &
1140 units=
"nondim", default=cs%Res_coef_khth)
1141 call get_param(param_file, mdl,
"VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1142 "The power of dx/Ld in the Kh resolution function. Any "//&
1143 "positive integer may be used, although even integers "//&
1144 "are more efficient to calculate. Setting this greater "//&
1145 "than 100 results in a step-function being used. "//&
1146 "This function affects lateral viscosity, Kh, and not KhTh.", &
1147 units=
"nondim", default=cs%Res_fn_power_khth)
1148 call get_param(param_file, mdl,
"INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1149 "If true, interpolate the resolution function to the "//&
1150 "velocity points from the thickness points; otherwise "//&
1151 "interpolate the wave speed and calculate the resolution "//&
1152 "function independently at each point.", default=.true.)
1153 if (cs%interpolate_Res_fn)
then
1154 if (cs%Res_coef_visc /= cs%Res_coef_khth)
call mom_error(fatal, &
1155 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1156 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1157 if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth)
call mom_error(fatal, &
1158 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1159 "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1162 call get_param(param_file, mdl,
"GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1163 "If true, uses Gill's definition of the baroclinic "//&
1164 "equatorial deformation radius, otherwise, if false, use "//&
1165 "Pedlosky's definition. These definitions differ by a factor "//&
1166 "of 2 in front of the beta term in the denominator. Gill's "//&
1167 "is the more appropriate definition.", default=.false.)
1168 if (gill_equatorial_ld)
then
1172 do j=js-1,jeq ;
do i=is-1,ieq
1173 cs%f2_dx2_q(i,j) = (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * &
1174 max(g%CoriolisBu(i,j)**2, absurdly_small_freq**2)
1175 cs%beta_dx2_q(i,j) = oneortwo * ((g%dxBu(i,j))**2 + (g%dyBu(i,j))**2) * (sqrt(0.5 * &
1176 ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1177 ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1178 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1179 ((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) ) ))
1182 do j=js,je ;
do i=is-1,ieq
1183 cs%f2_dx2_u(i,j) = (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * &
1184 max(0.5* (g%CoriolisBu(i,j)**2+g%CoriolisBu(i,j-1)**2), absurdly_small_freq**2)
1185 cs%beta_dx2_u(i,j) = oneortwo * ((g%dxCu(i,j))**2 + (g%dyCu(i,j))**2) * (sqrt( &
1186 0.25*( (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2 + &
1187 ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1188 (((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2 + &
1189 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) ) + &
1190 ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 ))
1193 do j=js-1,jeq ;
do i=is,ie
1194 cs%f2_dx2_v(i,j) = ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * &
1195 max(0.5*(g%CoriolisBu(i,j)**2+g%CoriolisBu(i-1,j)**2), absurdly_small_freq**2)
1196 cs%beta_dx2_v(i,j) = oneortwo * ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * (sqrt( &
1197 ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1198 0.25*( (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1199 ((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2) + &
1200 (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2 + &
1201 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1206 if (cs%Depth_scaled_KhTh)
then
1207 cs%calculate_depth_fns = .true.
1208 allocate(cs%Depth_fn_u(isdb:iedb,jsd:jed)) ; cs%Depth_fn_u(:,:) = 0.0
1209 allocate(cs%Depth_fn_v(isd:ied,jsdb:jedb)) ; cs%Depth_fn_v(:,:) = 0.0
1210 call get_param(param_file, mdl,
"DEPTH_SCALED_KHTH_H0", cs%depth_scaled_khth_h0, &
1211 "The depth above which KHTH is scaled away.",&
1212 units=
"m", default=1000.)
1213 call get_param(param_file, mdl,
"DEPTH_SCALED_KHTH_EXP", cs%depth_scaled_khth_exp, &
1214 "The exponent used in the depth dependent scaling function for KHTH.",&
1215 units=
"nondim", default=3.0)
1219 cs%id_Rd_dx = register_diag_field(
'ocean_model',
'Rd_dx', diag%axesT1, time, &
1220 'Ratio between deformation radius and grid spacing',
'm m-1')
1221 cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
1223 if (cs%calculate_Rd_dx)
then
1224 cs%calculate_cg1 = .true.
1225 allocate(cs%Rd_dx_h(isd:ied,jsd:jed)) ; cs%Rd_dx_h(:,:) = 0.0
1226 allocate(cs%beta_dx2_h(isd:ied,jsd:jed)); cs%beta_dx2_h(:,:) = 0.0
1227 allocate(cs%f2_dx2_h(isd:ied,jsd:jed)) ; cs%f2_dx2_h(:,:) = 0.0
1228 do j=js-1,je+1 ;
do i=is-1,ie+1
1229 cs%f2_dx2_h(i,j) = (g%dxT(i,j)**2 + g%dyT(i,j)**2) * &
1230 max(0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1231 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)), &
1232 absurdly_small_freq**2)
1233 cs%beta_dx2_h(i,j) = oneortwo * ((g%dxT(i,j))**2 + (g%dyT(i,j))**2) * (sqrt(0.5 * &
1234 ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1235 ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1236 (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1237 ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1241 if (cs%calculate_cg1)
then
1243 allocate(cs%cg1(isd:ied,jsd:jed)) ; cs%cg1(:,:) = 0.0
1244 call wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, mono_n2_depth=n2_filter_depth)
1248 call get_param(param_file, mdl,
"USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
1249 "If true, use the QG Leith viscosity as the GM coefficient.", &
1252 if (cs%Use_QG_Leith_GM)
then
1253 call get_param(param_file, mdl,
"LEITH_LAP_CONST", leith_lap_const, &
1254 "The nondimensional Laplacian Leith constant, \n"//&
1255 "often set to 1.0", units=
"nondim", default=0.0)
1257 call get_param(param_file, mdl,
"USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
1258 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1261 alloc_(cs%Laplac3_const_u(isdb:iedb,jsd:jed)) ; cs%Laplac3_const_u(:,:) = 0.0
1262 alloc_(cs%Laplac3_const_v(isd:ied,jsdb:jedb)) ; cs%Laplac3_const_v(:,:) = 0.0
1263 alloc_(cs%KH_u_QG(isdb:iedb,jsd:jed,g%ke)) ; cs%KH_u_QG(:,:,:) = 0.0
1264 alloc_(cs%KH_v_QG(isd:ied,jsdb:jedb,g%ke)) ; cs%KH_v_QG(:,:,:) = 0.0
1267 cs%id_KH_u_QG = register_diag_field(
'ocean_model',
'KH_u_QG', diag%axesCuL, time, &
1268 'Horizontal viscosity from Leith QG, at u-points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1269 cs%id_KH_v_QG = register_diag_field(
'ocean_model',
'KH_v_QG', diag%axesCvL, time, &
1270 'Horizontal viscosity from Leith QG, at v-points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1272 do j=jsq,jeq+1 ;
do i=is-1,ieq
1274 grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
1275 grid_sp_u3 = sqrt(grid_sp_u2)
1276 cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
1278 do j=js-1,jeq ;
do i=isq,ieq+1
1281 grid_sp_v2 = g%dyCv(i,j)*g%dxCu(i,j)
1282 grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
1283 cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
1286 if (.not. cs%use_stored_slopes)
call mom_error(fatal, &
1287 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1288 "USE_STORED_SLOPES must be True when using QG Leith.")
1293 cs%use_variable_mixing = .true.