7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
8 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
41 implicit none ;
private
43 #include <MOM_memory.h>
59 logical :: bulkmixedlayer
65 logical :: bottomdraglaw
67 logical :: bbl_mixing_as_max
71 logical :: use_lotw_bbl_diffusivity
72 logical :: lotw_bbl_use_omega
89 logical :: limit_dissipation
98 logical :: ml_radiation
111 real :: ml_rad_kd_max
113 real :: ml_rad_efold_coeff
117 logical :: ml_rad_bug
120 logical :: ml_rad_tke_decay
129 logical :: ml_use_omega
132 real :: ml_omega_frac
135 logical :: user_change_diff
136 logical :: usekappashear
138 logical :: vertex_shear
140 logical :: use_cvmix_shear
143 logical :: use_cvmix_ddiff
144 logical :: simple_tke_to_kd
146 real :: max_rrho_salt_fingers
147 real :: max_salt_diff_salt_fingers
150 logical :: answers_2018
154 character(len=200) :: inputdir
164 integer :: id_maxtke = -1, id_tke_to_kd = -1, id_kd_user = -1
165 integer :: id_kd_layer = -1, id_kd_bbl = -1, id_n2 = -1
166 integer :: id_kd_work = -1, id_kt_extra = -1, id_ks_extra = -1
173 real,
pointer,
dimension(:,:,:) :: &
179 kt_extra => null(), &
181 real,
pointer,
dimension(:,:,:) :: tke_to_kd => null()
204 subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, &
205 G, GV, US, CS, Kd_lay, Kd_int)
209 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
211 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
213 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
215 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
217 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
221 type(
forcing),
intent(in) :: fluxes
222 type(optics_type),
pointer :: optics
226 real,
intent(in) :: dt
228 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
229 intent(out) :: kd_lay
230 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
231 optional,
intent(out) :: kd_int
234 real,
dimension(SZI_(G)) :: &
239 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
242 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
246 real,
dimension(SZI_(G),SZK_(G)) :: &
247 n2_lay, & !< squared buoyancy frequency associated with layers [T-2 ~> s-2]
248 maxtke, & !< energy required to entrain to h_max [Z3 T-3 ~> m3 s-3]
253 real,
dimension(SZI_(G),SZK_(G)+1) :: &
254 n2_int, & !< squared buoyancy frequency associated at interfaces [T-2 ~> s-2]
255 drho_int, & !< locally ref potential density difference across interfaces [R ~> kg m-3]
256 kt_extra, & !< double difusion diffusivity of temperature [Z2 T-1 ~> m2 s-1]
263 integer :: kb(szi_(g))
265 logical :: showcalltree
267 integer :: i, j, k, is, ie, js, je, nz
268 integer :: isd, ied, jsd, jed
270 real :: kappa_dt_fill
272 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
273 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
274 showcalltree = calltree_showquery()
275 if (showcalltree)
call calltree_enter(
"set_diffusivity(), MOM_set_diffusivity.F90")
277 if (.not.
associated(cs))
call mom_error(fatal,
"set_diffusivity: "//&
278 "Module must be initialized before it is used.")
280 if (cs%answers_2018)
then
282 kappa_dt_fill = us%m_to_Z**2 * 1.e-3 * 7200.
284 kappa_dt_fill = cs%Kd_smooth * dt
286 omega2 = cs%omega * cs%omega
288 use_eos =
associated(tv%eqn_of_state)
290 if ((cs%use_CVMix_ddiff .or. cs%double_diffusion) .and. .not. &
291 (
associated(visc%Kd_extra_T) .and.
associated(visc%Kd_extra_S))) &
292 call mom_error(fatal,
"set_diffusivity: both visc%Kd_extra_T and "//&
293 "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are true.")
297 kd_lay(:,:,:) = cs%Kd
298 kd_int(:,:,:) = cs%Kd
299 if (
associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = cs%Kv
303 if (cs%id_N2 > 0)
then
304 allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1)) ; dd%N2_3d(:,:,:) = 0.0
306 if (cs%id_Kd_user > 0)
then
307 allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1)) ; dd%Kd_user(:,:,:) = 0.0
309 if (cs%id_Kd_work > 0)
then
310 allocate(dd%Kd_work(isd:ied,jsd:jed,nz)) ; dd%Kd_work(:,:,:) = 0.0
312 if (cs%id_maxTKE > 0)
then
313 allocate(dd%maxTKE(isd:ied,jsd:jed,nz)) ; dd%maxTKE(:,:,:) = 0.0
315 if (cs%id_TKE_to_Kd > 0)
then
316 allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0
318 if (cs%id_KT_extra > 0)
then
319 allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0
321 if (cs%id_KS_extra > 0)
then
322 allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0
324 if (cs%id_Kd_BBL > 0)
then
325 allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0
334 call hchksum(tv%T,
"before vert_fill_TS tv%T",g%HI)
335 call hchksum(tv%S,
"before vert_fill_TS tv%S",g%HI)
336 call hchksum(h,
"before vert_fill_TS h",g%HI, scale=gv%H_to_m)
338 call vert_fill_ts(h, tv%T, tv%S, kappa_dt_fill, t_f, s_f, g, gv, larger_h_denom=.true.)
340 call hchksum(tv%T,
"after vert_fill_TS tv%T",g%HI)
341 call hchksum(tv%S,
"after vert_fill_TS tv%S",g%HI)
342 call hchksum(h,
"after vert_fill_TS h",g%HI, scale=gv%H_to_m)
346 if (cs%useKappaShear)
then
348 call hchksum_pair(
"before calc_KS [uv]_h", u_h, v_h, g%HI, scale=us%L_T_to_m_s)
351 if (cs%Vertex_shear)
then
353 (gv%Z_to_H**2)*kappa_dt_fill, halo=1)
356 visc%TKE_turb, visc%Kv_shear_Bu, dt, g, gv, us, cs%kappaShear_CSp)
357 if (
associated(visc%Kv_shear)) visc%Kv_shear(:,:,:) = 0.0
359 call hchksum(visc%Kd_shear,
"after calc_KS_vert visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
360 call bchksum(visc%Kv_shear_Bu,
"after calc_KS_vert visc%Kv_shear_Bu", g%HI, scale=us%Z2_T_to_m2_s)
361 call bchksum(visc%TKE_turb,
"after calc_KS_vert visc%TKE_turb", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
365 call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, &
366 visc%Kv_shear, dt, g, gv, us, cs%kappaShear_CSp)
368 call hchksum(visc%Kd_shear,
"after calc_KS visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
369 call hchksum(visc%Kv_shear,
"after calc_KS visc%Kv_shear", g%HI, scale=us%Z2_T_to_m2_s)
370 call hchksum(visc%TKE_turb,
"after calc_KS visc%TKE_turb", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
374 if (showcalltree)
call calltree_waypoint(
"done with calculate_kappa_shear (set_diffusivity)")
375 elseif (cs%use_CVMix_shear)
then
377 call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear, g, gv, us, cs%CVMix_shear_CSp)
379 call hchksum(visc%Kd_shear,
"after CVMix_shear visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
380 call hchksum(visc%Kv_shear,
"after CVMix_shear visc%Kv_shear", g%HI, scale=us%Z2_T_to_m2_s)
382 elseif (
associated(visc%Kv_shear))
then
383 visc%Kv_shear(:,:,:) = 0.0
398 call find_n2(h, tv, t_f, s_f, fluxes, j, g, gv, us, cs, drho_int, n2_lay, n2_int, n2_bot)
400 if (
associated(dd%N2_3d))
then
401 do k=1,nz+1 ;
do i=is,ie ; dd%N2_3d(i,j,k) = n2_int(i,k) ;
enddo ;
enddo
405 call calculate_bkgnd_mixing(h, tv, n2_lay, kd_lay, visc%Kv_slow, j, g, gv, us, cs%bkgnd_mixing_csp)
408 if (cs%double_diffusion)
then
409 call double_diffusion(tv, h, t_f, s_f, j, g, gv, us, cs, kt_extra, ks_extra)
410 do k=2,nz ;
do i=is,ie
411 if (ks_extra(i,k) > kt_extra(i,k))
then
412 kd_lay(i,j,k-1) = kd_lay(i,j,k-1) + 0.5 * kt_extra(i,k)
413 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * kt_extra(i,k)
414 visc%Kd_extra_S(i,j,k) = (ks_extra(i,k) - kt_extra(i,k))
415 visc%Kd_extra_T(i,j,k) = 0.0
416 elseif (kt_extra(i,k) > 0.0)
then
417 kd_lay(i,j,k-1) = kd_lay(i,j,k-1) + 0.5 * ks_extra(i,k)
418 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * ks_extra(i,k)
419 visc%Kd_extra_T(i,j,k) = (kt_extra(i,k) - ks_extra(i,k))
420 visc%Kd_extra_S(i,j,k) = 0.0
422 visc%Kd_extra_T(i,j,k) = 0.0
423 visc%Kd_extra_S(i,j,k) = 0.0
426 if (
associated(dd%KT_extra))
then ;
do k=1,nz+1 ;
do i=is,ie
427 dd%KT_extra(i,j,k) = kt_extra(i,k)
428 enddo ;
enddo ;
endif
430 if (
associated(dd%KS_extra))
then ;
do k=1,nz+1 ;
do i=is,ie
431 dd%KS_extra(i,j,k) = ks_extra(i,k)
432 enddo ;
enddo ;
endif
437 if (cs%use_CVMix_ddiff)
then
439 call compute_ddiff_coeffs(h, tv, g, gv, us, j, visc%Kd_extra_T, visc%Kd_extra_S, cs%CVMix_ddiff_csp)
444 if (cs%useKappaShear .or. cs%use_CVMix_shear)
then
445 if (
present(kd_int))
then
446 do k=2,nz ;
do i=is,ie
447 kd_int(i,j,k) = visc%Kd_shear(i,j,k) + 0.5 * (kd_lay(i,j,k-1) + kd_lay(i,j,k))
450 kd_int(i,j,1) = visc%Kd_shear(i,j,1)
451 kd_int(i,j,nz+1) = 0.0
454 do k=1,nz ;
do i=is,ie
455 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * (visc%Kd_shear(i,j,k) + visc%Kd_shear(i,j,k+1))
458 if (
present(kd_int))
then
460 kd_int(i,j,1) = kd_lay(i,j,1) ; kd_int(i,j,nz+1) = 0.0
462 do k=2,nz ;
do i=is,ie
463 kd_int(i,j,k) = 0.5 * (kd_lay(i,j,k-1) + kd_lay(i,j,k))
468 call find_tke_to_kd(h, tv, drho_int, n2_lay, j, dt, g, gv, us, cs, tke_to_kd, &
470 if (
associated(dd%maxTKE))
then ;
do k=1,nz ;
do i=is,ie
471 dd%maxTKE(i,j,k) = maxtke(i,k)
472 enddo ;
enddo ;
endif
473 if (
associated(dd%TKE_to_Kd))
then ;
do k=1,nz ;
do i=is,ie
474 dd%TKE_to_Kd(i,j,k) = tke_to_kd(i,k)
475 enddo ;
enddo ;
endif
478 if (cs%ML_radiation) &
482 call calculate_tidal_mixing(h, n2_bot, j, tke_to_kd, maxtke, g, gv, us, cs%tm_csp, &
483 n2_lay, n2_int, kd_lay, kd_int, cs%Kd_max, visc%Kv_slow)
487 if (cs%bottomdraglaw .and. (cs%BBL_effic>0.0))
then
488 if (cs%use_LOTW_BBL_diffusivity)
then
489 call add_lotw_bbl_diffusivity(h, u, v, tv, fluxes, visc, j, n2_int, g, gv, us, cs, &
490 kd_lay, kd_int, dd%Kd_BBL)
493 maxtke, kb, g, gv, us, cs, kd_lay, kd_int, dd%Kd_BBL)
497 if (cs%limit_dissipation)
then
503 do k=2,nz-1 ;
do i=is,ie
504 dissip = max( cs%dissip_min, &
505 cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_lay(i,k)), &
506 cs%dissip_N2 * n2_lay(i,k))
507 kd_lay(i,j,k) = max(kd_lay(i,j,k) , &
508 dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_lay(i,k) + omega2))))
511 if (
present(kd_int))
then ;
do k=2,nz ;
do i=is,ie
512 dissip = max( cs%dissip_min, &
513 cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_int(i,k)), &
514 cs%dissip_N2 * n2_int(i,k))
515 kd_int(i,j,k) = max(kd_int(i,j,k) , &
516 dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_int(i,k) + omega2))))
517 enddo ;
enddo ;
endif
520 if (
associated(dd%Kd_work))
then
521 do k=1,nz ;
do i=is,ie
522 dd%Kd_Work(i,j,k) = gv%Rho0 * kd_lay(i,j,k) * n2_lay(i,k) * &
529 call hchksum(kd_lay ,
"Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
531 if (cs%useKappaShear)
call hchksum(visc%Kd_shear,
"Turbulent Kd", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
533 if (cs%use_CVMix_ddiff)
then
534 call hchksum(visc%Kd_extra_T,
"MOM_set_diffusivity: Kd_extra_T", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
535 call hchksum(visc%Kd_extra_S,
"MOM_set_diffusivity: Kd_extra_S", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
538 if (
associated(visc%kv_bbl_u) .and.
associated(visc%kv_bbl_v))
then
539 call uvchksum(
"BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, &
540 g%HI, 0, symmetric=.true., scale=us%Z2_T_to_m2_s)
543 if (
associated(visc%bbl_thick_u) .and.
associated(visc%bbl_thick_v))
then
544 call uvchksum(
"BBL bbl_thick_[uv]", visc%bbl_thick_u, &
545 visc%bbl_thick_v, g%HI, 0, symmetric=.true., scale=us%Z_to_m)
548 if (
associated(visc%Ray_u) .and.
associated(visc%Ray_v))
then
549 call uvchksum(
"Ray_[uv]", visc%Ray_u, visc%Ray_v, g%HI, 0, symmetric=.true., scale=us%Z_to_m*us%s_to_T)
554 if (cs%Kd_add > 0.0)
then
555 if (
present(kd_int))
then
557 do k=1,nz ;
do j=js,je ;
do i=is,ie
558 kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add
559 kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add
560 enddo ;
enddo ;
enddo
563 do k=1,nz ;
do j=js,je ;
do i=is,ie
564 kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add
565 enddo ;
enddo ;
enddo
569 if (cs%user_change_diff)
then
570 call user_change_diff(h, tv, g, gv, cs%user_change_diff_CSp, kd_lay, kd_int, &
571 t_f, s_f, dd%Kd_user)
577 if (cs%bkgnd_mixing_csp%id_kd_bkgnd > 0) &
578 call post_data(cs%bkgnd_mixing_csp%id_kd_bkgnd, cs%bkgnd_mixing_csp%kd_bkgnd, cs%bkgnd_mixing_csp%diag)
579 if (cs%bkgnd_mixing_csp%id_kv_bkgnd > 0) &
580 call post_data(cs%bkgnd_mixing_csp%id_kv_bkgnd, cs%bkgnd_mixing_csp%kv_bkgnd, cs%bkgnd_mixing_csp%diag)
583 if (cs%CVMix_ddiff_csp%id_KT_extra > 0) &
584 call post_data(cs%CVMix_ddiff_csp%id_KT_extra, visc%Kd_extra_T, cs%CVMix_ddiff_csp%diag)
585 if (cs%CVMix_ddiff_csp%id_KS_extra > 0) &
586 call post_data(cs%CVMix_ddiff_csp%id_KS_extra, visc%Kd_extra_S, cs%CVMix_ddiff_csp%diag)
587 if (cs%CVMix_ddiff_csp%id_R_rho > 0) &
588 call post_data(cs%CVMix_ddiff_csp%id_R_rho, cs%CVMix_ddiff_csp%R_rho, cs%CVMix_ddiff_csp%diag)
590 if (cs%id_Kd_layer > 0)
call post_data(cs%id_Kd_layer, kd_lay, cs%diag)
595 if (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation .or. &
596 cs%tm_csp%Lowmode_itidal_dissipation)
then
598 if (cs%id_N2 > 0)
call post_data(cs%id_N2, dd%N2_3d, cs%diag)
599 if (cs%id_Kd_user > 0)
call post_data(cs%id_Kd_user, dd%Kd_user, cs%diag)
600 if (cs%id_Kd_Work > 0)
call post_data(cs%id_Kd_Work, dd%Kd_Work, cs%diag)
601 if (cs%id_maxTKE > 0)
call post_data(cs%id_maxTKE, dd%maxTKE, cs%diag)
602 if (cs%id_TKE_to_Kd > 0)
call post_data(cs%id_TKE_to_Kd, dd%TKE_to_Kd, cs%diag)
605 if (cs%id_KT_extra > 0)
call post_data(cs%id_KT_extra, dd%KT_extra, cs%diag)
606 if (cs%id_KS_extra > 0)
call post_data(cs%id_KS_extra, dd%KS_extra, cs%diag)
607 if (cs%id_Kd_BBL > 0)
call post_data(cs%id_Kd_BBL, dd%Kd_BBL, cs%diag)
609 if (
associated(dd%N2_3d))
deallocate(dd%N2_3d)
610 if (
associated(dd%Kd_work))
deallocate(dd%Kd_work)
611 if (
associated(dd%Kd_user))
deallocate(dd%Kd_user)
612 if (
associated(dd%maxTKE))
deallocate(dd%maxTKE)
613 if (
associated(dd%TKE_to_Kd))
deallocate(dd%TKE_to_Kd)
614 if (
associated(dd%KT_extra))
deallocate(dd%KT_extra)
615 if (
associated(dd%KS_extra))
deallocate(dd%KS_extra)
616 if (
associated(dd%Kd_BBL))
deallocate(dd%Kd_BBL)
623 subroutine find_tke_to_kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, &
624 TKE_to_Kd, maxTKE, kb)
628 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
632 real,
dimension(SZI_(G),SZK_(G)+1),
intent(in) :: dRho_int
634 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: N2_lay
636 integer,
intent(in) :: j
637 real,
intent(in) :: dt
639 real,
dimension(SZI_(G),SZK_(G)),
intent(out) :: TKE_to_Kd
644 real,
dimension(SZI_(G),SZK_(G)),
intent(out) :: maxTKE
646 integer,
dimension(SZI_(G)),
intent(out) :: kb
649 real,
dimension(SZI_(G),SZK_(G)) :: &
661 real,
dimension(SZI_(G)) :: &
681 logical :: do_i(SZI_(G))
683 integer :: i, k, is, ie, nz, i_rem, kmb, kb_min
684 is = g%isc ; ie = g%iec ; nz = g%ke
688 h_neglect = gv%H_subroundoff
689 g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
690 if (cs%answers_2018)
then
691 i_rho0 = 1.0 / (gv%Rho0)
692 g_irho0 = (us%L_to_Z**2 * gv%g_Earth) * i_rho0
698 if (cs%simple_TKE_to_Kd)
then
699 do k=1,nz ;
do i=is,ie
700 hn2po2 = (gv%H_to_Z * h(i,j,k)) * (n2_lay(i,k) + omega2)
702 tke_to_kd(i,k) = 1.0 / hn2po2
703 else; tke_to_kd(i,k) = 0.;
endif
706 maxtke(i,k) = hn2po2 * cs%Kd_max
713 if (cs%bulkmixedlayer)
then
714 kmb = gv%nk_rho_varies
715 do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ;
enddo
718 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
721 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
727 do k=kmb+1,nz-1 ;
if (rcv_kmb(i) <= gv%Rlay(k))
exit ;
enddo
732 do k=kb(i)-1,kmb+1,-1
733 if (rho_0(i,kmb) > rho_0(i,k))
exit
734 if (h(i,j,k)>2.0*gv%Angstrom_H) kb(i) = k
741 do i=is,ie ; kb(i) = 1 ;
enddo
747 do k=2,nz-1 ;
do i=is,ie
748 dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
750 do i=is,ie ; dsp1_ds(i,nz) = 0.0 ;
enddo
752 if (cs%bulkmixedlayer)
then
753 kmb = gv%nk_rho_varies
755 htot(i) = gv%H_to_Z*h(i,j,kmb)
758 mfkb(i) = ds_dsp1(i,kb(i)) * (gv%H_to_Z*(h(i,j,kmb) - gv%Angstrom_H))
760 do k=1,kmb-1 ;
do i=is,ie
761 htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
762 mfkb(i) = mfkb(i) + ds_dsp1(i,k+1)*(gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H))
766 maxent(i,1) = 0.0 ; htot(i) = gv%H_to_Z*(h(i,j,1) - gv%Angstrom_H)
769 do k=kb_min,nz-1 ;
do i=is,ie
771 maxent(i,kb(i)) = mfkb(i)
772 elseif (k > kb(i))
then
773 if (cs%answers_2018)
then
774 maxent(i,k) = (1.0/dsp1_ds(i,k))*(maxent(i,k-1) + htot(i))
776 maxent(i,k) = ds_dsp1(i,k)*(maxent(i,k-1) + htot(i))
778 htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
783 htot(i) = gv%H_to_Z*(h(i,j,nz) - gv%Angstrom_H) ; maxent(i,nz) = 0.0
784 do_i(i) = (g%mask2dT(i,j) > 0.5)
788 do i=is,ie ;
if (do_i(i))
then
789 if (k<kb(i))
then ; do_i(i) = .false. ; cycle ;
endif
791 maxent(i,k) = min(maxent(i,k),dsp1_ds(i,k+1)*maxent(i,k+1) + htot(i))
792 htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
799 maxtke(i,1) = 0.0 ; tke_to_kd(i,1) = 0.0
800 maxtke(i,nz) = 0.0 ; tke_to_kd(i,nz) = 0.0
802 do k=2,kmb ;
do i=is,ie
804 tke_to_kd(i,k) = 1.0 / ((n2_lay(i,k) + omega2) * &
805 (gv%H_to_Z*(h(i,j,k) + h_neglect)))
807 do k=kmb+1,kb_min-1 ;
do i=is,ie
810 maxtke(i,k) = 0.0 ; tke_to_kd(i,k) = 0.0
812 do k=kb_min,nz-1 ;
do i=is,ie
822 dh_max = maxent(i,k) * (1.0 + dsp1_ds(i,k))
823 drho_lay = 0.5 * max(drho_int(i,k) + drho_int(i,k+1), 0.0)
824 maxtke(i,k) = i_dt * (g_irho0 * &
825 (0.5*max(drho_int(i,k+1) + dsp1_ds(i,k)*drho_int(i,k), 0.0))) * &
826 ((gv%H_to_Z*h(i,j,k) + dh_max) * maxent(i,k))
830 tke_to_kd(i,k) = 1.0 / (g_rho0 * drho_lay + &
831 cs%omega**2 * gv%H_to_Z*(h(i,j,k) + h_neglect))
838 subroutine find_n2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
839 N2_lay, N2_int, N2_bot)
843 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
847 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
850 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
853 type(
forcing),
intent(in) :: fluxes
854 integer,
intent(in) :: j
856 real,
dimension(SZI_(G),SZK_(G)+1), &
857 intent(out) :: dRho_int
859 real,
dimension(SZI_(G),SZK_(G)+1), &
860 intent(out) :: N2_int
861 real,
dimension(SZI_(G),SZK_(G)), &
862 intent(out) :: N2_lay
863 real,
dimension(SZI_(G)),
intent(out) :: N2_bot
865 real,
dimension(SZI_(G),SZK_(G)+1) :: &
870 real,
dimension(SZI_(G)) :: &
885 logical :: do_i(SZI_(G)), do_any
886 integer :: i, k, is, ie, nz
888 is = g%isc ; ie = g%iec ; nz = g%ke
889 g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
890 h_neglect = gv%H_subroundoff
894 drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
895 drho_int_unfilt(i,1) = 0.0 ; drho_int_unfilt(i,nz+1) = 0.0
897 if (
associated(tv%eqn_of_state))
then
898 if (
associated(fluxes%p_surf))
then
899 do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ;
enddo
901 do i=is,ie ; pres(i) = 0.0 ;
enddo
905 pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
906 temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
907 salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
910 drho_dt(:,k), drho_ds(:,k), is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
912 drho_int(i,k) = max(drho_dt(i,k)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
913 drho_ds(i,k)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
914 drho_int_unfilt(i,k) = max(drho_dt(i,k)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + &
915 drho_ds(i,k)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0)
919 do k=2,nz ;
do i=is,ie
920 drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
925 do k=1,nz ;
do i=is,ie
926 n2_lay(i,k) = g_rho0 * 0.5*(drho_int(i,k) + drho_int(i,k+1)) / &
927 (gv%H_to_Z*(h(i,j,k) + h_neglect))
929 do i=is,ie ; n2_int(i,1) = 0.0 ; n2_int(i,nz+1) = 0.0 ;
enddo
930 do k=2,nz ;
do i=is,ie
931 n2_int(i,k) = g_rho0 * drho_int(i,k) / &
932 (0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k) + h_neglect))
937 hb(i) = 0.0 ; drho_bot(i) = 0.0
938 z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
939 do_i(i) = (g%mask2dT(i,j) > 0.5)
941 if ( (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation) .and. &
942 .not. cs%tm_csp%use_CVMix_tidal )
then
943 h_amp(i) = sqrt(cs%tm_csp%h2(i,j))
951 do i=is,ie ;
if (do_i(i))
then
952 dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
953 z_from_bot(i) = z_from_bot(i) + dz_int
955 hb(i) = hb(i) + dz_int
956 drho_bot(i) = drho_bot(i) + drho_int(i,k)
958 if (z_from_bot(i) > h_amp(i))
then
961 hb(i) = hb(i) + 0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k-2))
962 drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
969 if (.not.do_any)
exit
973 if (hb(i) > 0.0)
then
974 n2_bot(i) = (g_rho0 * drho_bot(i)) / hb(i)
975 else ; n2_bot(i) = 0.0 ;
endif
976 z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
977 do_i(i) = (g%mask2dT(i,j) > 0.5)
982 do i=is,ie ;
if (do_i(i))
then
983 dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
984 z_from_bot(i) = z_from_bot(i) + dz_int
986 n2_int(i,k) = n2_bot(i)
987 if (k>2) n2_lay(i,k-1) = n2_bot(i)
989 if (z_from_bot(i) > h_amp(i))
then
990 if (k>2) n2_int(i,k-1) = n2_bot(i)
996 if (.not.do_any)
exit
999 if (
associated(tv%eqn_of_state))
then
1000 do k=1,nz+1 ;
do i=is,ie
1001 drho_int(i,k) = drho_int_unfilt(i,k)
1014 subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, US, CS, Kd_T_dd, Kd_S_dd)
1020 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1022 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1025 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1028 integer,
intent(in) :: j
1030 real,
dimension(SZI_(G),SZK_(G)+1), &
1031 intent(out) :: Kd_T_dd
1033 real,
dimension(SZI_(G),SZK_(G)+1), &
1034 intent(out) :: Kd_S_dd
1037 real,
dimension(SZI_(G)) :: &
1052 real,
parameter :: Rrho0 = 1.9
1054 integer :: i, k, is, ie, nz
1055 is = g%isc ; ie = g%iec ; nz = g%ke
1057 if (
associated(tv%eqn_of_state))
then
1059 pres(i) = 0.0 ; kd_t_dd(i,1) = 0.0 ; kd_s_dd(i,1) = 0.0
1060 kd_t_dd(i,nz+1) = 0.0 ; kd_s_dd(i,nz+1) = 0.0
1064 pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
1065 temp_int(i) = 0.5 * (t_f(i,j,k-1) + t_f(i,j,k))
1066 salin_int(i) = 0.5 * (s_f(i,j,k-1) + s_f(i,j,k))
1069 drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1072 alpha_dt = -1.0*drho_dt(i) * (t_f(i,j,k-1) - t_f(i,j,k))
1073 beta_ds = drho_ds(i) * (s_f(i,j,k-1) - s_f(i,j,k))
1075 if ((alpha_dt > beta_ds) .and. (beta_ds > 0.0))
then
1076 rrho = min(alpha_dt / beta_ds, rrho0)
1077 diff_dd = 1.0 - ((rrho-1.0)/(rrho0-1.0))
1078 kd_dd = cs%Max_salt_diff_salt_fingers * diff_dd*diff_dd*diff_dd
1079 kd_t_dd(i,k) = 0.7 * kd_dd
1080 kd_s_dd(i,k) = kd_dd
1081 elseif ((alpha_dt < 0.) .and. (beta_ds < 0.) .and. (alpha_dt > beta_ds))
then
1082 rrho = alpha_dt / beta_ds
1083 kd_dd = cs%Kv_molecular * 0.909 * exp(4.6 * exp(-0.54 * (1/rrho - 1)))
1085 if (rrho > 0.5) prandtl = (1.85-0.85/rrho)*rrho
1086 kd_t_dd(i,k) = kd_dd
1087 kd_s_dd(i,k) = prandtl * kd_dd
1089 kd_t_dd(i,k) = 0.0 ; kd_s_dd(i,k) = 0.0
1099 maxTKE, kb, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL)
1103 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1105 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1107 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1111 type(
forcing),
intent(in) :: fluxes
1114 integer,
intent(in) :: j
1115 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: TKE_to_Kd
1120 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: maxTKE
1122 integer,
dimension(SZI_(G)),
intent(in) :: kb
1125 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1126 intent(inout) :: Kd_lay
1128 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1129 intent(inout) :: Kd_int
1131 real,
dimension(:,:,:),
pointer :: Kd_BBL
1135 real,
dimension(SZK_(G)+1) :: &
1137 real,
dimension(SZI_(G)) :: &
1148 real :: TKE_to_layer
1158 logical :: Rayleigh_drag
1162 logical :: domore, do_i(SZI_(G))
1163 logical :: do_diag_Kd_BBL
1165 integer :: i, k, is, ie, nz, i_rem, kb_min
1166 is = g%isc ; ie = g%iec ; nz = g%ke
1168 do_diag_kd_bbl =
associated(kd_bbl)
1170 if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0)))
return
1172 cdrag_sqrt = sqrt(cs%cdrag)
1173 tke_ray = 0.0 ; rayleigh_drag = .false.
1174 if (
associated(visc%Ray_u) .and.
associated(visc%Ray_v)) rayleigh_drag = .true.
1176 i_rho0 = 1.0 / (gv%Rho0)
1177 r0_g = gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)
1179 do k=2,nz ; rint(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ;
enddo
1181 kb_min = max(gv%nk_rho_varies+1,2)
1187 ustar_h = visc%ustar_BBL(i,j)
1188 if (
associated(fluxes%ustar_tidal)) &
1189 ustar_h = ustar_h + fluxes%ustar_tidal(i,j)
1190 absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1191 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))))
1192 if ((ustar_h > 0.0) .and. (absf > 0.5*cs%IMax_decay*ustar_h))
then
1193 i2decay(i) = absf / ustar_h
1197 i2decay(i) = 0.5*cs%IMax_decay
1199 tke(i) = ((cs%BBL_effic * cdrag_sqrt) * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))) ) * &
1202 if (
associated(fluxes%TKE_tidal)) &
1203 tke(i) = tke(i) + fluxes%TKE_tidal(i,j) * i_rho0 * &
1204 (cs%BBL_effic * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))))
1211 gh_sum_top(i) = r0_g * 400.0 * ustar_h**2
1213 do_i(i) = (g%mask2dT(i,j) > 0.5)
1214 htot(i) = gv%H_to_Z*h(i,j,nz)
1215 rho_htot(i) = gv%Rlay(nz)*(gv%H_to_Z*h(i,j,nz))
1216 rho_top(i) = gv%Rlay(1)
1217 if (cs%bulkmixedlayer .and. do_i(i)) rho_top(i) = gv%Rlay(kb(i)-1)
1220 do k=nz-1,2,-1 ; domore = .false.
1221 do i=is,ie ;
if (do_i(i))
then
1222 htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
1223 rho_htot(i) = rho_htot(i) + gv%Rlay(k)*(gv%H_to_Z*h(i,j,k))
1224 if (htot(i)*gv%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i)))
then
1226 rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i)
1228 elseif (k <= kb(i))
then ; do_i(i) = .false.
1229 else ; domore = .true. ;
endif
1231 if (.not.domore)
exit
1234 do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ;
enddo
1237 do i=is,ie ;
if (do_i(i))
then
1238 if (k<kb(i))
then ; do_i(i) = .false. ; cycle ;
endif
1242 tke(i) = tke(i) * exp(-i2decay(i) * (gv%H_to_Z*(h(i,j,k) + h(i,j,k+1))))
1245 if (maxtke(i,k) <= 0.0) cycle
1249 if (tke(i) > 0.0)
then
1250 if (rint(k) <= rho_top(i))
then
1251 tke_to_layer = tke(i)
1253 drl = rint(k+1) - rint(k) ; drbot = rint(k+1) - rho_top(i)
1254 tke_to_layer = tke(i) * drl * &
1255 (3.0*drbot*(rint(k) - rho_top(i)) + drl**2) / (drbot**3)
1257 else ; tke_to_layer = 0.0 ;
endif
1260 if (rayleigh_drag) tke_ray = 0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1261 ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1262 g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1263 (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1264 g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1266 if (tke_to_layer + tke_ray > 0.0)
then
1267 if (cs%BBL_mixing_as_max)
then
1268 if (tke_to_layer + tke_ray > maxtke(i,k)) &
1269 tke_to_layer = maxtke(i,k) - tke_ray
1271 tke(i) = tke(i) - tke_to_layer
1273 if (kd_lay(i,j,k) < (tke_to_layer + tke_ray) * tke_to_kd(i,k))
then
1274 delta_kd = (tke_to_layer + tke_ray) * tke_to_kd(i,k) - kd_lay(i,j,k)
1275 if ((cs%Kd_max >= 0.0) .and. (delta_kd > cs%Kd_max))
then
1276 delta_kd = cs%Kd_max
1277 kd_lay(i,j,k) = kd_lay(i,j,k) + delta_kd
1279 kd_lay(i,j,k) = (tke_to_layer + tke_ray) * tke_to_kd(i,k)
1281 kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * delta_kd
1282 kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * delta_kd
1283 if (do_diag_kd_bbl)
then
1284 kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1285 kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1289 if (kd_lay(i,j,k) >= maxtke(i,k) * tke_to_kd(i,k))
then
1291 tke(i) = tke(i) + tke_ray
1292 elseif (kd_lay(i,j,k) + (tke_to_layer + tke_ray) * tke_to_kd(i,k) > &
1293 maxtke(i,k) * tke_to_kd(i,k))
then
1294 tke_here = ((tke_to_layer + tke_ray) + kd_lay(i,j,k) / tke_to_kd(i,k)) - maxtke(i,k)
1295 tke(i) = (tke(i) - tke_here) + tke_ray
1297 tke_here = tke_to_layer + tke_ray
1298 tke(i) = tke(i) - tke_to_layer
1300 if (tke(i) < 0.0) tke(i) = 0.0
1302 if (tke_here > 0.0)
then
1303 delta_kd = tke_here * tke_to_kd(i,k)
1304 if (cs%Kd_max >= 0.0) delta_kd = min(delta_kd, cs%Kd_max)
1305 kd_lay(i,j,k) = kd_lay(i,j,k) + delta_kd
1306 kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * delta_kd
1307 kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * delta_kd
1308 if (do_diag_kd_bbl)
then
1309 kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1310 kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1320 if ((tke(i)<= 0.0) .and. (tke_ray == 0.0))
then
1321 do_i(i) = .false. ; i_rem = i_rem - 1
1325 if (i_rem == 0)
exit
1334 G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL)
1338 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1340 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1342 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1346 type(
forcing),
intent(in) :: fluxes
1349 integer,
intent(in) :: j
1350 real,
dimension(SZI_(G),SZK_(G)+1), &
1351 intent(in) :: N2_int
1353 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1354 intent(inout) :: Kd_lay
1355 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1356 intent(inout) :: Kd_int
1357 real,
dimension(:,:,:),
pointer :: Kd_BBL
1361 real :: TKE_remaining
1362 real :: TKE_consumed
1371 real :: total_thickness
1378 logical :: Rayleigh_drag
1380 integer :: i, k, km1
1381 real,
parameter :: von_karm = 0.41
1382 logical :: do_diag_Kd_BBL
1384 if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0)))
return
1385 do_diag_kd_bbl =
associated(kd_bbl)
1388 if (cs%LOTW_BBL_use_omega) n2_min = cs%omega**2
1391 rayleigh_drag = .false.
1392 if (
associated(visc%Ray_u) .and.
associated(visc%Ray_v)) rayleigh_drag = .true.
1393 i_rho0 = 1.0 / (gv%Rho0)
1394 cdrag_sqrt = sqrt(cs%cdrag)
1399 absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1400 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))))
1403 ustar = visc%ustar_BBL(i,j)
1408 if (
associated(fluxes%ustar_tidal)) ustar = ustar + fluxes%ustar_tidal(i,j)
1412 idecay = cs%IMax_decay
1413 if ((ustar > 0.0) .and. (absf > cs%IMax_decay * ustar)) idecay = absf / ustar
1418 tke_column = cdrag_sqrt * visc%TKE_BBL(i,j)
1421 if (
associated(fluxes%TKE_tidal)) &
1422 tke_column = tke_column + fluxes%TKE_tidal(i,j) * i_rho0
1423 tke_column = cs%BBL_effic * tke_column
1425 tke_remaining = tke_column
1426 total_thickness = ( sum(h(i,j,:)) + gv%H_subroundoff )* gv%H_to_Z
1427 ustar_d = ustar * total_thickness
1434 dh = gv%H_to_Z * h(i,j,k)
1436 dhm1 = gv%H_to_Z * h(i,j,km1)
1439 if (rayleigh_drag) tke_remaining = tke_remaining + &
1440 0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1441 ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1442 g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1443 (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1444 g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1448 tke_remaining = exp(-idecay*dh) * tke_remaining
1450 z_bot = z_bot + h(i,j,k)*gv%H_to_Z
1451 d_minus_z = max(total_thickness - z_bot, 0.)
1455 if ( ustar_d + absf * ( z_bot * d_minus_z ) == 0.)
then
1458 kd_wall = ((von_karm * ustar2) * (z_bot * d_minus_z)) &
1459 / (ustar_d + absf * (z_bot * d_minus_z))
1464 tke_kd_wall = kd_wall * 0.5 * (dh + dhm1) * max(n2_int(i,k), n2_min)
1467 if (tke_kd_wall > 0.)
then
1468 tke_consumed = min(tke_kd_wall, tke_remaining)
1469 kd_wall = (tke_consumed / tke_kd_wall) * kd_wall
1472 if (tke_remaining > 0.)
then
1481 tke_remaining = tke_remaining - tke_consumed
1484 kd_int(i,j,k) = kd_int(i,j,k) + kd_wall
1485 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * (kd_wall + kd_lower)
1487 if (do_diag_kd_bbl) kd_bbl(i,j,k) = kd_wall
1498 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1500 type(
forcing),
intent(in) :: fluxes
1502 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1503 intent(inout) :: Kd_lay
1504 integer,
intent(in) :: j
1505 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: TKE_to_Kd
1510 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1511 optional,
intent(inout) :: Kd_int
1516 real,
dimension(SZI_(G)) :: h_ml
1517 real,
dimension(SZI_(G)) :: TKE_ml_flux
1518 real,
dimension(SZI_(G)) :: I_decay
1519 real,
dimension(SZI_(G)) :: Kd_mlr_ml
1529 real :: I_decay_len2_TKE
1533 logical :: do_any, do_i(SZI_(G))
1534 integer :: i, k, is, ie, nz, kml
1535 is = g%isc ; ie = g%iec ; nz = g%ke
1537 omega2 = cs%omega**2
1540 h_neglect = gv%H_subroundoff*gv%H_to_Z
1542 if (.not.cs%ML_radiation)
return
1544 do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (g%mask2dT(i,j) > 0.5) ;
enddo
1545 do k=1,kml ;
do i=is,ie ; h_ml(i) = h_ml(i) + gv%H_to_Z*h(i,j,k) ;
enddo ;
enddo
1547 do i=is,ie ;
if (do_i(i))
then
1548 if (cs%ML_omega_frac >= 1.0)
then
1551 f_sq = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1552 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
1553 if (cs%ML_omega_frac > 0.0) &
1554 f_sq = cs%ML_omega_frac * 4.0 * omega2 + (1.0 - cs%ML_omega_frac) * f_sq
1557 ustar_sq = max(fluxes%ustar(i,j), cs%ustar_min)**2
1559 tke_ml_flux(i) = (cs%mstar * cs%ML_rad_coeff) * (ustar_sq * (fluxes%ustar(i,j)))
1560 i_decay_len2_tke = cs%TKE_decay**2 * (f_sq / ustar_sq)
1562 if (cs%ML_rad_TKE_decay) &
1563 tke_ml_flux(i) = tke_ml_flux(i) * exp(-h_ml(i) * sqrt(i_decay_len2_tke))
1566 h_ml_sq = (cs%ML_rad_efold_coeff * (h_ml(i)+h_neglect))**2
1567 i_decay(i) = sqrt((i_decay_len2_tke * h_ml_sq + 1.0) / h_ml_sq)
1571 z1 = (gv%H_to_Z*h(i,j,kml+1)) * i_decay(i)
1573 kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (1.0 - exp(-z1))
1575 kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (z1 * (1.0 - z1 * (0.5 - c1_6 * z1)))
1577 kd_mlr_ml(i) = min(kd_mlr, cs%ML_rad_kd_max)
1578 tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1581 do k=1,kml+1 ;
do i=is,ie ;
if (do_i(i))
then
1582 kd_lay(i,j,k) = kd_lay(i,j,k) + kd_mlr_ml(i)
1583 endif ;
enddo ;
enddo
1584 if (
present(kd_int))
then
1585 do k=2,kml+1 ;
do i=is,ie ;
if (do_i(i))
then
1586 kd_int(i,j,k) = kd_int(i,j,k) + kd_mlr_ml(i)
1587 endif ;
enddo ;
enddo
1588 if (kml<=nz-1)
then ;
do i=is,ie ;
if (do_i(i))
then
1589 kd_int(i,j,kml+2) = kd_int(i,j,kml+2) + 0.5 * kd_mlr_ml(i)
1590 endif ;
enddo ;
endif
1595 do i=is,ie ;
if (do_i(i))
then
1596 dzl = gv%H_to_Z*h(i,j,k) ; z1 = dzl*i_decay(i)
1597 if (cs%ML_Rad_bug)
then
1602 kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * &
1603 us%m_to_Z * ((1.0 - exp(-z1)) / dzl)
1605 kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * &
1606 us%m_to_Z * (i_decay(i) * (1.0 - z1 * (0.5 - c1_6*z1)))
1610 kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (1.0 - exp(-z1))
1612 kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (z1 * (1.0 - z1 * (0.5 - c1_6*z1)))
1615 kd_mlr = min(kd_mlr, cs%ML_rad_kd_max)
1616 kd_lay(i,j,k) = kd_lay(i,j,k) + kd_mlr
1617 if (
present(kd_int))
then
1618 kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * kd_mlr
1619 kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * kd_mlr
1622 tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1623 if (tke_ml_flux(i) * i_decay(i) < 0.1 * cs%Kd_min * omega2)
then
1625 else ; do_any = .true. ;
endif
1627 if (.not.do_any)
exit
1634 subroutine set_bbl_tke(u, v, h, fluxes, visc, G, GV, US, CS)
1638 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1640 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1642 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1644 type(
forcing),
intent(in) :: fluxes
1652 real,
dimension(SZI_(G)) :: &
1656 real,
dimension(SZIB_(G)) :: &
1661 real :: vhtot(szi_(g))
1663 real,
dimension(SZI_(G),SZJB_(G)) :: &
1670 logical :: domore, do_i(szi_(g))
1671 integer :: i, j, k, is, ie, js, je, nz
1672 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1674 if (.not.
associated(cs))
call mom_error(fatal,
"set_BBL_TKE: "//&
1675 "Module must be initialized before it is used.")
1677 if (.not.cs%bottomdraglaw .or. (cs%BBL_effic<=0.0))
then
1678 if (
associated(visc%ustar_BBL))
then
1679 do j=js,je ;
do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ;
enddo ;
enddo
1681 if (
associated(visc%TKE_BBL))
then
1682 do j=js,je ;
do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ;
enddo ;
enddo
1687 cdrag_sqrt = sqrt(cs%cdrag)
1697 do i=is,ie ;
if ((g%mask2dCv(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_v(i,j) > 0.0))
then
1698 do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0
1699 vstar(i,j) = visc%Kv_bbl_v(i,j) / (cdrag_sqrt*visc%bbl_thick_v(i,j))
1701 do_i(i) = .false. ; vstar(i,j) = 0.0 ; htot(i) = 0.0
1705 do i=is,ie ;
if (do_i(i))
then
1706 hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j+1,k))
1707 if ((htot(i) + hvel) >= visc%bbl_thick_v(i,j))
then
1708 vhtot(i) = vhtot(i) + (visc%bbl_thick_v(i,j) - htot(i))*v(i,j,k)
1709 htot(i) = visc%bbl_thick_v(i,j)
1712 vhtot(i) = vhtot(i) + hvel*v(i,j,k)
1713 htot(i) = htot(i) + hvel
1717 if (.not.domore)
exit
1719 do i=is,ie ;
if ((g%mask2dCv(i,j) > 0.5) .and. (htot(i) > 0.0))
then
1720 v2_bbl(i,j) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i))
1727 do i=is-1,ie ;
if ((g%mask2dCu(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_u(i,j) > 0.0))
then
1728 do_i(i) = .true. ; uhtot(i) = 0.0 ; htot(i) = 0.0
1729 ustar(i) = visc%Kv_bbl_u(i,j) / (cdrag_sqrt*visc%bbl_thick_u(i,j))
1731 do_i(i) = .false. ; ustar(i) = 0.0 ; htot(i) = 0.0
1733 do k=nz,1,-1 ; domore = .false.
1734 do i=is-1,ie ;
if (do_i(i))
then
1735 hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i+1,j,k))
1736 if ((htot(i) + hvel) >= visc%bbl_thick_u(i,j))
then
1737 uhtot(i) = uhtot(i) + (visc%bbl_thick_u(i,j) - htot(i))*u(i,j,k)
1738 htot(i) = visc%bbl_thick_u(i,j)
1741 uhtot(i) = uhtot(i) + hvel*u(i,j,k)
1742 htot(i) = htot(i) + hvel
1746 if (.not.domore)
exit
1748 do i=is-1,ie ;
if ((g%mask2dCu(i,j) > 0.5) .and. (htot(i) > 0.0))
then
1749 u2_bbl(i) = (uhtot(i)*uhtot(i))/(htot(i)*htot(i))
1755 visc%ustar_BBL(i,j) = sqrt(0.5*g%IareaT(i,j) * &
1756 ((g%areaCu(i-1,j)*(ustar(i-1)*ustar(i-1)) + &
1757 g%areaCu(i,j)*(ustar(i)*ustar(i))) + &
1758 (g%areaCv(i,j-1)*(vstar(i,j-1)*vstar(i,j-1)) + &
1759 g%areaCv(i,j)*(vstar(i,j)*vstar(i,j))) ) )
1760 visc%TKE_BBL(i,j) = us%L_to_Z**2 * &
1761 (((g%areaCu(i-1,j)*(ustar(i-1)*u2_bbl(i-1)) + &
1762 g%areaCu(i,j) * (ustar(i)*u2_bbl(i))) + &
1763 (g%areaCv(i,j-1)*(vstar(i,j-1)*v2_bbl(i,j-1)) + &
1764 g%areaCv(i,j) * (vstar(i,j)*v2_bbl(i,j))) )*g%IareaT(i,j))
1774 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1779 integer,
dimension(SZI_(G)),
intent(in) :: kb
1784 integer,
intent(in) :: j
1785 real,
dimension(SZI_(G),SZK_(G)),
intent(out) :: ds_dsp1
1789 real,
dimension(SZI_(G),SZK_(G)), &
1790 optional,
intent(in) :: rho_0
1796 real :: a(SZK_(G)), a_0(SZK_(G))
1797 real :: p_ref(SZI_(G))
1798 real :: Rcv(SZI_(G),SZK_(G))
1801 integer :: i, k, k3, is, ie, nz, kmb
1802 is = g%isc ; ie = g%iec ; nz = g%ke
1805 if (gv%g_prime(k+1) /= 0.0)
then
1807 ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
1816 if (cs%bulkmixedlayer)
then
1817 g_r0 = gv%g_Earth / (gv%Rho0)
1818 kmb = gv%nk_rho_varies
1820 do i=is,ie ; p_ref(i) = tv%P_Ref ;
enddo
1823 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1826 if (kb(i) <= nz-1)
then
1831 i_drho = g_r0 / gv%g_prime(k+1)
1834 a(k3+1) = (gv%Rlay(k) - rcv(i,k3)) * i_drho
1836 if ((
present(rho_0)) .and. (a(kmb+1) < 2.0*eps*ds_dsp1(i,k)))
then
1840 if ((rho_0(i,k) > rho_0(i,kmb)) .and. &
1841 (rho_0(i,k+1) > rho_0(i,k)))
then
1842 i_drho = 1.0 / (rho_0(i,k+1)-rho_0(i,k))
1843 a_0(kmb+1) = min((rho_0(i,k)-rho_0(i,kmb)) * i_drho, ds_dsp1(i,k))
1844 if (a_0(kmb+1) > a(kmb+1))
then
1846 a_0(k3) = a_0(kmb+1) + (rho_0(i,kmb)-rho_0(i,k3-1)) * i_drho
1848 if (a(kmb+1) <= eps*ds_dsp1(i,k))
then
1849 do k3=2,kmb+1 ; a(k3) = a_0(k3) ;
enddo
1852 tmp = a(kmb+1)/(eps*ds_dsp1(i,k)) - 1.0
1853 do k3=2,kmb+1 ; a(k3) = tmp*a(k3) + (1.0-tmp)*a_0(k3) ;
enddo
1859 ds_dsp1(i,k) = max(a(kmb+1),1e-5)
1868 ds_dsp1(i,k3) = max(a(k3),ds_dsp1(i,k))
1878 type(time_type),
intent(in) :: time
1884 type(
diag_ctrl),
target,
intent(inout) :: diag
1891 integer,
optional,
intent(out) :: halo_ts
1895 real :: decay_length
1896 logical :: ml_use_omega
1897 logical :: default_2018_answers
1899 # include "version_variable.h"
1900 character(len=40) :: mdl =
"MOM_set_diffusivity"
1901 real :: omega_frac_dflt
1902 integer :: i, j, is, ie, js, je
1903 integer :: isd, ied, jsd, jed
1905 if (
associated(cs))
then
1906 call mom_error(warning,
"diabatic_entrain_init called with an associated "// &
1907 "control structure.")
1912 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1913 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1916 if (
associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp
1917 if (
associated(tm_csp)) cs%tm_csp => tm_csp
1920 cs%BBL_mixing_as_max = .true.
1921 cs%cdrag = 0.003 ; cs%BBL_effic = 0.0
1922 cs%bulkmixedlayer = (gv%nkml > 0)
1927 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".")
1928 cs%inputdir = slasher(cs%inputdir)
1929 call get_param(param_file, mdl,
"FLUX_RI_MAX", cs%FluxRi_max, &
1930 "The flux Richardson number where the stratification is "//&
1931 "large enough that N2 > omega2. The full expression for "//&
1932 "the Flux Richardson number is usually "//&
1933 "FLUX_RI_MAX*N2/(N2+OMEGA2).", default=0.2)
1934 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
1935 "The rotation rate of the earth.", units=
"s-1", &
1936 default=7.2921e-5, scale=us%T_to_s)
1938 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1939 "This sets the default value for the various _2018_ANSWERS parameters.", &
1941 call get_param(param_file, mdl,
"SET_DIFF_2018_ANSWERS", cs%answers_2018, &
1942 "If true, use the order of arithmetic and expressions that recover the "//&
1943 "answers from the end of 2018. Otherwise, use updated and more robust "//&
1944 "forms of the same expressions.", default=default_2018_answers)
1946 call get_param(param_file, mdl,
"ML_RADIATION", cs%ML_radiation, &
1947 "If true, allow a fraction of TKE available from wind "//&
1948 "work to penetrate below the base of the mixed layer "//&
1949 "with a vertical decay scale determined by the minimum "//&
1950 "of: (1) The depth of the mixed layer, (2) an Ekman "//&
1951 "length scale.", default=.false.)
1952 if (cs%ML_radiation)
then
1954 cs%ustar_min = 2e-4 * cs%omega * (gv%Angstrom_Z + gv%H_subroundoff*gv%H_to_Z)
1956 call get_param(param_file, mdl,
"ML_RAD_EFOLD_COEFF", cs%ML_rad_efold_coeff, &
1957 "A coefficient that is used to scale the penetration "//&
1958 "depth for turbulence below the base of the mixed layer. "//&
1959 "This is only used if ML_RADIATION is true.", units=
"nondim", &
1961 call get_param(param_file, mdl,
"ML_RAD_BUG", cs%ML_rad_bug, &
1962 "If true use code with a bug that reduces the energy available "//&
1963 "in the transition layer by a factor of the inverse of the energy "//&
1964 "deposition lenthscale (in m).", default=.true.)
1965 call get_param(param_file, mdl,
"ML_RAD_KD_MAX", cs%ML_rad_kd_max, &
1966 "The maximum diapycnal diffusivity due to turbulence "//&
1967 "radiated from the base of the mixed layer. "//&
1968 "This is only used if ML_RADIATION is true.", &
1969 units=
"m2 s-1", default=1.0e-3, &
1970 scale=us%m2_s_to_Z2_T)
1971 call get_param(param_file, mdl,
"ML_RAD_COEFF", cs%ML_rad_coeff, &
1972 "The coefficient which scales MSTAR*USTAR^3 to obtain "//&
1973 "the energy available for mixing below the base of the "//&
1974 "mixed layer. This is only used if ML_RADIATION is true.", &
1975 units=
"nondim", default=0.2)
1976 call get_param(param_file, mdl,
"ML_RAD_APPLY_TKE_DECAY", cs%ML_rad_TKE_decay, &
1977 "If true, apply the same exponential decay to ML_rad as "//&
1978 "is applied to the other surface sources of TKE in the "//&
1979 "mixed layer code. This is only used if ML_RADIATION is true.", &
1981 call get_param(param_file, mdl,
"MSTAR", cs%mstar, &
1982 "The ratio of the friction velocity cubed to the TKE "//&
1983 "input to the mixed layer.",
"units=nondim", default=1.2)
1984 call get_param(param_file, mdl,
"TKE_DECAY", cs%TKE_decay, &
1985 "The ratio of the natural Ekman depth to the TKE decay scale.", &
1986 units=
"nondim", default=2.5)
1987 call get_param(param_file, mdl,
"ML_USE_OMEGA", ml_use_omega, &
1988 "If true, use the absolute rotation rate instead of the "//&
1989 "vertical component of rotation when setting the decay "//&
1990 "scale for turbulence.", default=.false., do_not_log=.true.)
1991 omega_frac_dflt = 0.0
1992 if (ml_use_omega)
then
1993 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1994 omega_frac_dflt = 1.0
1996 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%ML_omega_frac, &
1997 "When setting the decay scale for turbulence, use this "//&
1998 "fraction of the absolute rotation rate blended with the "//&
1999 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
2000 units=
"nondim", default=omega_frac_dflt)
2003 call get_param(param_file, mdl,
"BOTTOMDRAGLAW", cs%bottomdraglaw, &
2004 "If true, the bottom stress is calculated with a drag "//&
2005 "law of the form c_drag*|u|*u. The velocity magnitude "//&
2006 "may be an assumed value or it may be based on the "//&
2007 "actual velocity in the bottommost HBBL, depending on "//&
2008 "LINEAR_DRAG.", default=.true.)
2009 if (cs%bottomdraglaw)
then
2010 call get_param(param_file, mdl,
"CDRAG", cs%cdrag, &
2011 "The drag coefficient relating the magnitude of the "//&
2012 "velocity field to the bottom stress. CDRAG is only used "//&
2013 "if BOTTOMDRAGLAW is true.", units=
"nondim", default=0.003)
2014 call get_param(param_file, mdl,
"BBL_EFFIC", cs%BBL_effic, &
2015 "The efficiency with which the energy extracted by "//&
2016 "bottom drag drives BBL diffusion. This is only "//&
2017 "used if BOTTOMDRAGLAW is true.", units=
"nondim", default=0.20)
2018 call get_param(param_file, mdl,
"BBL_MIXING_MAX_DECAY", decay_length, &
2019 "The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//&
2020 "to penetrate as far as stratification and rotation permit. The default "//&
2021 "for now is 200 m. This is only used if BOTTOMDRAGLAW is true.", &
2022 units=
"m", default=200.0, scale=us%m_to_Z)
2025 if (decay_length > 0.0) cs%IMax_decay = 1.0/decay_length
2026 call get_param(param_file, mdl,
"BBL_MIXING_AS_MAX", cs%BBL_mixing_as_max, &
2027 "If true, take the maximum of the diffusivity from the "//&
2028 "BBL mixing and the other diffusivities. Otherwise, "//&
2029 "diffusivity from the BBL_mixing is simply added.", &
2031 call get_param(param_file, mdl,
"USE_LOTW_BBL_DIFFUSIVITY", cs%use_LOTW_BBL_diffusivity, &
2032 "If true, uses a simple, imprecise but non-coordinate dependent, model "//&
2033 "of BBL mixing diffusivity based on Law of the Wall. Otherwise, uses "//&
2034 "the original BBL scheme.", default=.false.)
2035 if (cs%use_LOTW_BBL_diffusivity)
then
2036 call get_param(param_file, mdl,
"LOTW_BBL_USE_OMEGA", cs%LOTW_BBL_use_omega, &
2037 "If true, use the maximum of Omega and N for the TKE to diffusion "//&
2038 "calculation. Otherwise, N is N.", default=.true.)
2041 cs%use_LOTW_BBL_diffusivity = .false.
2044 'Bottom Boundary Layer Diffusivity',
'm2 s-1', &
2045 conversion=us%Z2_T_to_m2_s)
2046 call get_param(param_file, mdl,
"SIMPLE_TKE_TO_KD", cs%simple_TKE_to_Kd, &
2047 "If true, uses a simple estimate of Kd/TKE that will "//&
2048 "work for arbitrary vertical coordinates. If false, "//&
2049 "calculates Kd/TKE and bounds based on exact energetics "//&
2050 "for an isopycnal layer-formulation.", &
2054 call bkgnd_mixing_init(time, g, gv, us, param_file, cs%diag, cs%bkgnd_mixing_csp)
2056 call get_param(param_file, mdl,
"KV", cs%Kv, &
2057 "The background kinematic viscosity in the interior. "//&
2058 "The molecular value, ~1e-6 m2 s-1, may be used.", &
2059 units=
"m2 s-1", scale=us%m2_s_to_Z2_T, &
2060 fail_if_missing=.true.)
2062 call get_param(param_file, mdl,
"KD", cs%Kd, &
2063 "The background diapycnal diffusivity of density in the "//&
2064 "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
2065 "may be used.", units=
"m2 s-1", scale=us%m2_s_to_Z2_T, &
2066 fail_if_missing=.true.)
2067 call get_param(param_file, mdl,
"KD_MIN", cs%Kd_min, &
2068 "The minimum diapycnal diffusivity.", &
2069 units=
"m2 s-1", default=0.01*cs%Kd*us%Z2_T_to_m2_s, &
2070 scale=us%m2_s_to_Z2_T)
2071 call get_param(param_file, mdl,
"KD_MAX", cs%Kd_max, &
2072 "The maximum permitted increment for the diapycnal "//&
2073 "diffusivity from TKE-based parameterizations, or a "//&
2074 "negative value for no limit.", units=
"m2 s-1", default=-1.0, &
2075 scale=us%m2_s_to_Z2_T)
2076 if (cs%simple_TKE_to_Kd .and. cs%Kd_max<=0.)
call mom_error(fatal, &
2077 "set_diffusivity_init: To use SIMPLE_TKE_TO_KD, KD_MAX must be set to >0.")
2078 call get_param(param_file, mdl,
"KD_ADD", cs%Kd_add, &
2079 "A uniform diapycnal diffusivity that is added "//&
2080 "everywhere without any filtering or scaling.", &
2081 units=
"m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2082 if (cs%use_LOTW_BBL_diffusivity .and. cs%Kd_max<=0.)
call mom_error(fatal, &
2083 "set_diffusivity_init: KD_MAX must be set (positive) when "// &
2084 "USE_LOTW_BBL_DIFFUSIVITY=True.")
2085 call get_param(param_file, mdl,
"KD_SMOOTH", cs%Kd_smooth, &
2086 "A diapycnal diffusivity that is used to interpolate "//&
2087 "more sensible values of T & S into thin layers.", &
2088 default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
2090 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
2091 "If true, write out verbose debugging data.", &
2092 default=.false., debuggingparam=.true.)
2094 call get_param(param_file, mdl,
"USER_CHANGE_DIFFUSIVITY", cs%user_change_diff, &
2095 "If true, call user-defined code to change the diffusivity.", &
2098 call get_param(param_file, mdl,
"DISSIPATION_MIN", cs%dissip_min, &
2099 "The minimum dissipation by which to determine a lower "//&
2100 "bound of Kd (a floor).", units=
"W m-3", default=0.0, &
2101 scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*(us%T_to_s**2))
2102 call get_param(param_file, mdl,
"DISSIPATION_N0", cs%dissip_N0, &
2103 "The intercept when N=0 of the N-dependent expression "//&
2104 "used to set a minimum dissipation by which to determine "//&
2105 "a lower bound of Kd (a floor): A in eps_min = A + B*N.", &
2106 units=
"W m-3", default=0.0, &
2107 scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*(us%T_to_s**2))
2108 call get_param(param_file, mdl,
"DISSIPATION_N1", cs%dissip_N1, &
2109 "The coefficient multiplying N, following Gargett, used to "//&
2110 "set a minimum dissipation by which to determine a lower "//&
2111 "bound of Kd (a floor): B in eps_min = A + B*N", &
2112 units=
"J m-3", default=0.0, scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*us%T_to_s)
2113 call get_param(param_file, mdl,
"DISSIPATION_KD_MIN", cs%dissip_Kd_min, &
2114 "The minimum vertical diffusivity applied as a floor.", &
2115 units=
"m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2117 cs%limit_dissipation = (cs%dissip_min>0.) .or. (cs%dissip_N1>0.) .or. &
2118 (cs%dissip_N0>0.) .or. (cs%dissip_Kd_min>0.)
2120 if (cs%FluxRi_max > 0.0) &
2121 cs%dissip_N2 = cs%dissip_Kd_min * gv%Rho0 / cs%FluxRi_max
2124 'Diapycnal diffusivity of layers (as set)',
'm2 s-1', &
2125 conversion=us%Z2_T_to_m2_s)
2128 if (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation .or. &
2129 cs%tm_csp%Lowmode_itidal_dissipation)
then
2132 'Work done by Diapycnal Mixing',
'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2134 'Maximum layer TKE',
'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
2136 'Convert TKE to Kd',
's2 m', &
2137 conversion=us%Z2_T_to_m2_s*(us%m_to_Z**3*us%T_to_s**3))
2139 'Buoyancy frequency squared',
's-2', cmor_field_name=
'obvfsq', &
2140 cmor_long_name=
'Square of seawater buoyancy frequency', &
2141 cmor_standard_name=
'square_of_brunt_vaisala_frequency_in_sea_water', &
2142 conversion=us%s_to_T**2)
2144 if (cs%user_change_diff) &
2146 'User-specified Extra Diffusivity',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
2149 call get_param(param_file, mdl,
"DOUBLE_DIFFUSION", cs%double_diffusion, &
2150 "If true, increase diffusivites for temperature or salt "//&
2151 "based on double-diffusive parameterization from MOM4/KPP.", &
2154 if (cs%double_diffusion)
then
2155 call get_param(param_file, mdl,
"MAX_RRHO_SALT_FINGERS", cs%Max_Rrho_salt_fingers, &
2156 "Maximum density ratio for salt fingering regime.", &
2157 default=2.55, units=
"nondim")
2158 call get_param(param_file, mdl,
"MAX_SALT_DIFF_SALT_FINGERS", cs%Max_salt_diff_salt_fingers, &
2159 "Maximum salt diffusivity for salt fingering regime.", &
2160 default=1.e-4, units=
"m2 s-1", scale=us%m2_s_to_Z2_T)
2161 call get_param(param_file, mdl,
"KV_MOLECULAR", cs%Kv_molecular, &
2162 "Molecular viscosity for calculation of fluxes under "//&
2163 "double-diffusive convection.", default=1.5e-6, units=
"m2 s-1", &
2164 scale=us%m2_s_to_Z2_T)
2168 'Double-diffusive diffusivity for temperature',
'm2 s-1', &
2169 conversion=us%Z2_T_to_m2_s)
2172 'Double-diffusive diffusivity for salinity',
'm2 s-1', &
2173 conversion=us%Z2_T_to_m2_s)
2176 if (cs%user_change_diff)
then
2177 call user_change_diff_init(time, g, gv, us, param_file, diag, cs%user_change_diff_CSp)
2180 if (cs%tm_csp%Int_tide_dissipation .and. cs%bkgnd_mixing_csp%Bryan_Lewis_diffusivity) &
2181 call mom_error(fatal,
"MOM_Set_Diffusivity: "// &
2182 "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.")
2184 cs%useKappaShear = kappa_shear_init(time, g, gv, us, param_file, cs%diag, cs%kappaShear_CSp)
2187 if (cs%useKappaShear) &
2191 cs%use_CVMix_shear = cvmix_shear_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_shear_csp)
2194 cs%use_CVMix_ddiff = cvmix_ddiff_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_ddiff_csp)
2195 if (cs%use_CVMix_ddiff) &
2198 if (
present(halo_ts))
then
2200 if (cs%Vertex_Shear) halo_ts = 1
2209 if (.not.
associated(cs))
return
2213 if (cs%user_change_diff) &
2216 if (cs%use_CVMix_shear) &
2219 if (cs%use_CVMix_ddiff) &
2220 call cvmix_ddiff_end(cs%CVMix_ddiff_csp)
2222 if (
associated(cs))
deallocate(cs)