6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
7 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
19 implicit none ;
private
21 #include <MOM_memory.h>
67 integer :: max_rino_it
71 logical :: dkdq_iteration_bug
74 logical :: ks_at_vertex
76 logical :: eliminate_massless
83 logical :: debug = .false.
87 integer :: id_kd_shear = -1, id_tke = -1, id_ild2 = -1, id_dz_int = -1
94 #undef ADD_DIAGNOSTICS
100 kv_io, dt, G, GV, US, CS, initialize_all)
104 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
106 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
108 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
113 real,
dimension(:,:),
pointer :: p_surf
114 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
115 intent(inout) :: kappa_io
119 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
120 intent(out) :: tke_io
122 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
123 intent(inout) :: kv_io
127 real,
intent(in) :: dt
130 logical,
optional,
intent(in) :: initialize_all
134 real,
dimension(SZI_(G),SZK_(GV)) :: &
138 real,
dimension(SZI_(G),SZK_(GV)+1) :: &
141 real,
dimension(SZK_(GV)) :: &
148 real,
dimension(SZK_(GV)+1) :: &
159 logical :: use_temperature
161 logical :: new_kappa = .true.
164 integer,
dimension(SZK_(GV)+1) :: kc
167 real,
dimension(SZK_(GV)+1) :: kf
169 integer :: is, ie, js, je, i, j, k, nz, nzc
172 #ifdef ADD_DIAGNOSTICS
173 real,
dimension(SZK_(GV)+1) :: &
175 real,
dimension(SZI_(G),SZK_(GV)+1) :: &
177 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
180 is = g%isc ; ie = g%iec; js = g%jsc ; je = g%jec ; nz = gv%ke
182 use_temperature = .false. ;
if (
associated(tv%T)) use_temperature = .true.
183 new_kappa = .true. ;
if (
present(initialize_all)) new_kappa = initialize_all
186 dz_massless = 0.1*sqrt(k0dt)
189 #ifdef ADD_DIAGNOSTICS
194 do k=1,nz ;
do i=is,ie
195 h_2d(i,k) = h(i,j,k)*gv%H_to_Z
196 u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
198 if (use_temperature)
then ;
do k=1,nz ;
do i=is,ie
199 t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
200 enddo ;
enddo ;
else ;
do k=1,nz ;
do i=is,ie
201 rho_2d(i,k) = gv%Rlay(k)
202 enddo ;
enddo ;
endif
203 if (.not.new_kappa)
then ;
do k=1,nz+1 ;
do i=is,ie
204 kappa_2d(i,k) = kappa_io(i,j,k)
205 enddo ;
enddo ;
endif
210 do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
214 if (cs%eliminate_massless)
then
218 dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
219 t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
223 if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
224 (h_2d(i,k) > dz_massless)) nzc = nzc+1
231 dz(nzc) = dz(nzc) + h_2d(i,k)
232 u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
233 v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
234 if (use_temperature)
then
235 t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
236 s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
238 t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
239 s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
245 do k=1,nzc ; idz(k) = 1.0 / dz(k) ;
enddo
249 kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
251 if (kc(k) > kc(k-1))
then
252 kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
254 kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
261 u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
263 if (use_temperature)
then
265 t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
269 t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
273 do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ;
enddo
275 f2 = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
276 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
277 surface_pres = 0.0 ;
if (
associated(p_surf)) surface_pres = p_surf(i,j)
284 do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ;
enddo
286 do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k) ;
enddo
289 #ifdef ADD_DIAGNOSTICS
291 dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
292 tke_avg, tv, cs, gv, us, i_ld2_1d, dz_int_1d)
295 dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
296 tke_avg, tv, cs, gv, us)
303 kappa_2d(i,k) = kappa_avg(k)
309 if (kf(k) == 0.0)
then
310 kappa_2d(i,k) = kappa_avg(kc(k))
311 tke_2d(i,k) = tke_avg(kc(k))
313 kappa_2d(i,k) = (1.0-kf(k)) * kappa_avg(kc(k)) + &
314 kf(k) * kappa_avg(kc(k)+1)
315 tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + &
316 kf(k) * tke_avg(kc(k)+1)
320 #ifdef ADD_DIAGNOSTICS
322 i_ld2_2d(i,k) = i_ld2_1d(k) ; dz_int_2d(i,k) = dz_int_1d(k)
328 kappa_2d(i,k) = 0.0 ; tke_2d(i,k) = 0.0
329 #ifdef ADD_DIAGNOSTICS
330 i_ld2_2d(i,k) = 0.0 ; dz_int_2d(i,k) = 0.0
335 do k=1,nz+1 ;
do i=is,ie
336 kappa_io(i,j,k) = g%mask2dT(i,j) * kappa_2d(i,k)
337 tke_io(i,j,k) = g%mask2dT(i,j) * tke_2d(i,k)
338 kv_io(i,j,k) = ( g%mask2dT(i,j) * kappa_2d(i,k) ) * cs%Prandtl_turb
339 #ifdef ADD_DIAGNOSTICS
340 i_ld2_3d(i,j,k) = i_ld2_2d(i,k) ; dz_int_3d(i,j,k) = dz_int_2d(i,k)
347 call hchksum(kappa_io,
"kappa", g%HI, scale=us%Z2_T_to_m2_s)
348 call hchksum(tke_io,
"tke", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
351 if (cs%id_Kd_shear > 0)
call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
352 if (cs%id_TKE > 0)
call post_data(cs%id_TKE, tke_io, cs%diag)
353 #ifdef ADD_DIAGNOSTICS
354 if (cs%id_ILd2 > 0)
call post_data(cs%id_ILd2, i_ld2_3d, cs%diag)
355 if (cs%id_dz_Int > 0)
call post_data(cs%id_dz_Int, dz_int_3d, cs%diag)
363 kv_io, dt, G, GV, US, CS, initialize_all)
367 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
369 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
371 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
373 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
375 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
380 real,
dimension(:,:),
pointer :: p_surf
382 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
383 intent(out) :: kappa_io
385 real,
dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
386 intent(out) :: tke_io
388 real,
dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
389 intent(inout) :: kv_io
393 real,
intent(in) :: dt
396 logical,
optional,
intent(in) :: initialize_all
400 real,
dimension(SZIB_(G),SZK_(GV)) :: &
404 real,
dimension(SZIB_(G),SZK_(GV)+1,2) :: &
406 real,
dimension(SZIB_(G),SZK_(GV)+1) :: &
408 real,
dimension(SZK_(GV)) :: &
415 real,
dimension(SZK_(GV)+1) :: &
428 logical :: use_temperature
430 logical :: new_kappa = .true.
434 integer,
dimension(SZK_(GV)+1) :: kc
437 real,
dimension(SZK_(GV)+1) :: kf
439 integer :: isb, ieb, jsb, jeb, i, j, k, nz, nzc, j2, j2m1
442 #ifdef ADD_DIAGNOSTICS
443 real,
dimension(SZK_(GV)+1) :: &
445 real,
dimension(SZI_(G),SZK_(GV)+1) :: &
447 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
450 isb = g%isc-1 ; ieb = g%iecB ; jsb = g%jsc-1 ; jeb = g%jecB ; nz = gv%ke
452 use_temperature = .false. ;
if (
associated(tv%T)) use_temperature = .true.
453 new_kappa = .true. ;
if (
present(initialize_all)) new_kappa = initialize_all
456 dz_massless = 0.1*sqrt(k0dt)
457 i_prandtl = 0.0 ;
if (cs%Prandtl_turb > 0.0) i_prandtl = 1.0 / cs%Prandtl_turb
460 #ifdef ADD_DIAGNOSTICS
465 j2 = mod(j,2)+1 ; j2m1 = 3-j2
468 do k=1,nz ;
do i=isb,ieb
469 u_2d(i,k) = (u_in(i,j,k) * (g%mask2dCu(i,j) * (h(i,j,k) + h(i+1,j,k))) + &
470 u_in(i,j+1,k) * (g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) ) / &
471 ((g%mask2dCu(i,j) * (h(i,j,k) + h(i+1,j,k)) + &
472 g%mask2dCu(i,j+1) * (h(i,j+1,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
473 v_2d(i,k) = (v_in(i,j,k) * (g%mask2dCv(i,j) * (h(i,j,k) + h(i,j+1,k))) + &
474 v_in(i+1,j,k) * (g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) ) / &
475 ((g%mask2dCv(i,j) * (h(i,j,k) + h(i,j+1,k)) + &
476 g%mask2dCv(i+1,j) * (h(i+1,j,k) + h(i+1,j+1,k))) + gv%H_subroundoff)
477 i_hwt = 1.0 / (((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
478 (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k))) + &
480 if (use_temperature)
then
481 t_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * t_in(i,j,k) + &
482 (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * t_in(i+1,j+1,k)) + &
483 ((g%mask2dT(i+1,j) * h(i+1,j,k)) * t_in(i+1,j,k) + &
484 (g%mask2dT(i,j+1) * h(i,j+1,k)) * t_in(i,j+1,k)) ) * i_hwt
485 s_2d(i,k) = ( ((g%mask2dT(i,j) * h(i,j,k)) * s_in(i,j,k) + &
486 (g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) * s_in(i+1,j+1,k)) + &
487 ((g%mask2dT(i+1,j) * h(i+1,j,k)) * s_in(i+1,j,k) + &
488 (g%mask2dT(i,j+1) * h(i,j+1,k)) * s_in(i,j+1,k)) ) * i_hwt
490 h_2d(i,k) = gv%H_to_Z * ((g%mask2dT(i,j) * h(i,j,k) + g%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
491 (g%mask2dT(i+1,j) * h(i+1,j,k) + g%mask2dT(i,j+1) * h(i,j+1,k)) ) / &
492 ((g%mask2dT(i,j) + g%mask2dT(i+1,j+1)) + &
493 (g%mask2dT(i+1,j) + g%mask2dT(i,j+1)) + 1.0e-36 )
498 if (.not.use_temperature)
then ;
do k=1,nz ;
do i=isb,ieb
499 rho_2d(i,k) = gv%Rlay(k)
500 enddo ;
enddo ;
endif
501 if (.not.new_kappa)
then ;
do k=1,nz+1 ;
do i=isb,ieb
502 kappa_2d(i,k,j2) = kv_io(i,j,k) * i_prandtl
503 enddo ;
enddo ;
endif
508 do i=isb,ieb ;
if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
509 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0)
then
513 if (cs%eliminate_massless)
then
517 dz(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
518 t0xdz(k) = 0.0 ; s0xdz(k) = 0.0
522 if ((k>cs%nkml) .and. (dz(nzc) > 0.0) .and. &
523 (h_2d(i,k) > dz_massless)) nzc = nzc+1
530 dz(nzc) = dz(nzc) + h_2d(i,k)
531 u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
532 v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
533 if (use_temperature)
then
534 t0xdz(nzc) = t0xdz(nzc) + t_2d(i,k)*h_2d(i,k)
535 s0xdz(nzc) = s0xdz(nzc) + s_2d(i,k)*h_2d(i,k)
537 t0xdz(nzc) = t0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
538 s0xdz(nzc) = s0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
544 do k=1,nzc ; idz(k) = 1.0 / dz(k) ;
enddo
548 kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
550 if (kc(k) > kc(k-1))
then
551 kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
553 kf(k) = dz_in_lay*idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
560 u0xdz(k) = u_2d(i,k)*dz(k) ; v0xdz(k) = v_2d(i,k)*dz(k)
562 if (use_temperature)
then
564 t0xdz(k) = t_2d(i,k)*dz(k) ; s0xdz(k) = s_2d(i,k)*dz(k)
568 t0xdz(k) = rho_2d(i,k)*dz(k) ; s0xdz(k) = rho_2d(i,k)*dz(k)
572 do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ;
enddo
574 f2 = g%CoriolisBu(i,j)**2
575 surface_pres = 0.0 ;
if (
associated(p_surf)) &
576 surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + &
577 (p_surf(i+1,j) + p_surf(i,j+1)))
583 do k=1,nzc+1 ; kappa(k) = us%m2_s_to_Z2_T*1.0 ;
enddo
585 do k=1,nzc+1 ; kappa(k) = kappa_2d(i,k,j2) ;
enddo
588 #ifdef ADD_DIAGNOSTICS
590 dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
591 tke_avg, tv, cs, gv, us, i_ld2_1d, dz_int_1d)
594 dz, u0xdz, v0xdz, t0xdz, s0xdz, kappa_avg, &
595 tke_avg, tv, cs, gv, us)
601 kappa_2d(i,k,j2) = kappa_avg(k)
607 if (kf(k) == 0.0)
then
608 kappa_2d(i,k,j2) = kappa_avg(kc(k))
609 tke_2d(i,k) = tke_avg(kc(k))
611 kappa_2d(i,k,j2) = (1.0-kf(k)) * kappa_avg(kc(k)) + kf(k) * kappa_avg(kc(k)+1)
612 tke_2d(i,k) = (1.0-kf(k)) * tke_avg(kc(k)) + kf(k) * tke_avg(kc(k)+1)
616 #ifdef ADD_DIAGNOSTICS
618 i_ld2_2d(i,k) = i_ld2_1d(k) ; dz_int_2d(i,k) = dz_int_1d(k)
624 kappa_2d(i,k,j2) = 0.0 ; tke_2d(i,k) = 0.0
625 #ifdef ADD_DIAGNOSTICS
626 i_ld2_2d(i,k) = 0.0 ; dz_int_2d(i,k) = 0.0
631 do k=1,nz+1 ;
do i=isb,ieb
632 tke_io(i,j,k) = g%mask2dBu(i,j) * tke_2d(i,k)
633 kv_io(i,j,k) = ( g%mask2dBu(i,j) * kappa_2d(i,k,j2) ) * cs%Prandtl_turb
634 #ifdef ADD_DIAGNOSTICS
635 i_ld2_3d(i,j,k) = i_ld2_2d(i,k) ; dz_int_3d(i,j,k) = dz_int_2d(i,k)
638 if (j>=g%jsc)
then ;
do k=1,nz+1 ;
do i=g%isc,g%iec
640 kappa_io(i,j,k) = g%mask2dT(i,j) * 0.25 * &
641 ((kappa_2d(i-1,k,j2m1) + kappa_2d(i,k,j2)) + &
642 (kappa_2d(i-1,k,j2) + kappa_2d(i,k,j2m1)))
643 enddo ;
enddo ;
endif
648 call hchksum(kappa_io,
"kappa", g%HI, scale=us%Z2_T_to_m2_s)
649 call bchksum(tke_io,
"tke", g%HI)
652 if (cs%id_Kd_shear > 0)
call post_data(cs%id_Kd_shear, kappa_io, cs%diag)
653 if (cs%id_TKE > 0)
call post_data(cs%id_TKE, tke_io, cs%diag)
654 #ifdef ADD_DIAGNOSTICS
655 if (cs%id_ILd2 > 0)
call post_data(cs%id_ILd2, i_ld2_3d, cs%diag)
656 if (cs%id_dz_Int > 0)
call post_data(cs%id_dz_Int, dz_int_3d, cs%diag)
664 dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
665 tke_avg, tv, CS, GV, US, I_Ld2_1d, dz_Int_1d)
668 real,
dimension(SZK_(GV)+1), &
669 intent(inout) :: kappa
670 real,
dimension(SZK_(GV)+1), &
673 integer,
intent(in) :: nzc
674 real,
intent(in) :: f2
675 real,
intent(in) :: surface_pres
676 real,
dimension(SZK_(GV)), &
678 real,
dimension(SZK_(GV)), &
680 real,
dimension(SZK_(GV)), &
682 real,
dimension(SZK_(GV)), &
684 real,
dimension(SZK_(GV)), &
686 real,
dimension(SZK_(GV)+1), &
687 intent(out) :: kappa_avg
688 real,
dimension(SZK_(GV)+1), &
689 intent(out) :: tke_avg
690 real,
intent(in) :: dt
696 real,
dimension(SZK_(GV)+1), &
697 optional,
intent(out) :: I_Ld2_1d
698 real,
dimension(SZK_(GV)+1), &
699 optional,
intent(out) :: dz_Int_1d
702 real,
dimension(nzc) :: &
711 real,
dimension(nzc+1) :: &
746 real :: dist_from_bot
758 real :: tol_dksrc_low
772 logical :: use_temperature
774 integer :: ks_kappa, ke_kappa
775 integer :: dt_halvings
778 integer :: dt_refinements
781 integer :: k, itt, itt_dt
783 integer :: max_debug_itt ;
parameter(max_debug_itt=20)
784 real :: wt(SZK_(GV)+1), wt_tot, I_wt_tot, wt_itt
785 real,
dimension(SZK_(GV)+1) :: &
786 Ri_k, tke_prev, dtke, dkappa, dtke_norm, &
789 real,
dimension(SZK_(GV)+1,0:max_debug_itt) :: &
790 tke_it1, N2_it1, Sh2_it1, ksrc_it1, kappa_it1, kprev_it1
791 real,
dimension(SZK_(GV)+1,1:max_debug_itt) :: &
792 dkappa_it1, wt_it1, K_Q_it1, d_dkappa_it1, dkappa_norm
793 real,
dimension(SZK_(GV),0:max_debug_itt) :: &
794 u_it1, v_it1, rho_it1, T_it1, S_it1
795 real,
dimension(0:max_debug_itt) :: &
796 dk_wt_it1, dkpos_wt_it1, dkneg_wt_it1, k_mag
797 real,
dimension(max_debug_itt) :: dt_it1
800 ri_crit = cs%Rino_crit
801 gr0 = gv%z_to_H*gv%H_to_Pa
802 g_r0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
806 tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*cs%kappa_tol_err
808 use_temperature = .false. ;
if (
associated(tv%T)) use_temperature = .true.
812 do k=1,nzc ; idz(k) = 1.0 / dz(k) ;
enddo
815 i_dz_int(1) = 2.0 / dz(1)
816 dist_from_top(1) = 0.0
818 i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
819 dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
821 i_dz_int(nzc+1) = 2.0 / dz(nzc)
826 a1(2) = k0dt*i_dz_int(2)
827 b1 = 1.0 / (dz(1) + a1(2))
828 u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
829 t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
830 c1(2) = a1(2) * b1 ; d1 = dz(1) * b1
832 bd1 = dz(k) + d1*a1(k)
833 a1(k+1) = k0dt*i_dz_int(k+1)
834 b1 = 1.0 / (bd1 + a1(k+1))
835 u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
836 v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
837 t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
838 sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
839 c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1
844 b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
845 u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
846 v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
848 b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
849 t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
850 sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
852 u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
853 t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
857 b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
858 u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
860 t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
869 dz_int(1) = 0.0 ; dz_int(2) = dz(1)
871 norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
872 dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
873 dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
875 dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
879 dist_from_bot = dist_from_bot + dz(k)
880 i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
881 (dist_from_top(k) * dist_from_bot)**2
885 if (use_temperature)
then
886 pressure(1) = surface_pres
888 pressure(k) = pressure(k-1) + gr0*dz(k-1)
889 t_int(k) = 0.5*(t(k-1) + t(k))
890 sal_int(k) = 0.5*(sal(k-1) + sal(k))
893 dbuoy_ds, 2, nzc-1, tv%eqn_of_state, scale=-g_r0*us%kg_m3_to_R)
895 do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ;
enddo
899 n2_debug(1) = 0.0 ; n2_debug(nzc+1) = 0.0
901 n2_debug(k) = max((dbuoy_dt(k) * (t0xdz(k-1)*idz(k-1) - t0xdz(k)*idz(k)) + &
902 dbuoy_ds(k) * (s0xdz(k-1)*idz(k-1) - s0xdz(k)*idz(k))) * &
906 u_it1(k,0) = u0xdz(k)*idz(k) ; v_it1(k,0) = v0xdz(k)*idz(k)
907 t_it1(k,0) = t0xdz(k)*idz(k) ; s_it1(k,0) = s0xdz(k)*idz(k)
910 kprev_it1(k,0) = kappa(k) ; kappa_it1(k,0) = kappa(k)
912 n2_it1(k,0) = n2_debug(k) ; sh2_it1(k,0) = s2(k) ; ksrc_it1(k,0) = k_src(k)
915 u_it1(k,0) = 0.0 ; v_it1(k,0) = 0.0
916 t_it1(k,0) = 0.0 ; s_it1(k,0) = 0.0
917 kprev_it1(k+1,0) = 0.0 ; kappa_it1(k+1,0) = 0.0 ; tke_it1(k+1,0) = 0.0
918 n2_it1(k+1,0) = 0.0 ; sh2_it1(k+1,0) = 0.0 ; ksrc_it1(k+1,0) = 0.0
920 do itt=1,max_debug_itt
923 u_it1(k,itt) = 0.0 ; v_it1(k,itt) = 0.0
924 t_it1(k,itt) = 0.0 ; s_it1(k,itt) = 0.0
928 kprev_it1(k,itt) = 0.0 ; kappa_it1(k,itt) = 0.0 ; tke_it1(k,itt) = 0.0
929 n2_it1(k,itt) = 0.0 ; sh2_it1(k,itt) = 0.0
930 ksrc_it1(k,itt) = 0.0
931 dkappa_it1(k,itt) = 0.0 ; wt_it1(k,itt) = 0.0
932 k_q_it1(k,itt) = 0.0 ; d_dkappa_it1(k,itt) = 0.0
935 do k=1,gv%ke+1 ; ksrc_av(k) = 0.0 ;
enddo
940 dbuoy_dt, dbuoy_ds, u, v, t, sal, gv, us, &
941 n2=n2, s2=s2, vel_underflow=cs%vel_underflow)
948 kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
949 local_src_avg(k) = 0.0
951 if ( dz_int(k) > 0.0 ) &
952 local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
958 do itt=1,cs%max_KS_it
965 ri_k(k) = 1e3 ;
if (s2(k) > 1e-3*n2(k)) ri_k(k) = n2(k) / s2(k)
966 if (itt > 1)
then ; tke_prev(k) = tke(k) ;
else ; tke_prev(k) = 0.0 ;
endif
972 nzc, cs, gv, us, k_q, tke, kappa_out, kappa_src, local_src)
977 ks_kappa = gv%ke+1 ; ke_kappa = 0
978 do k=2,nzc ;
if (kappa_out(k) > 0.0)
then
981 do k=nzc,ks_kappa,-1 ;
if (kappa_out(k) > 0.0)
then
984 if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
990 if ((ke_kappa < ks_kappa) .or. (itt==cs%max_RiNo_it))
then
996 tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
997 tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
998 tol_chg(k) = tol2 * local_src_avg(k)
1001 do itt_dt=1,(cs%max_KS_it+1-itt)/2
1010 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1011 gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1012 vel_underflow=cs%vel_underflow)
1014 idtt = 1.0 / dt_test
1015 do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1016 if (n2(k) < ri_crit * s2(k))
then
1017 k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1018 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1019 if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1020 (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))))
then
1021 valid_dt = .false. ;
exit
1024 if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))
then
1025 valid_dt = .false. ; k_src(k) = 0.0 ;
exit
1031 dt_test = 0.5*dt_test
1033 if ((dt_test < dt_rem) .and. valid_dt)
then
1034 dt_inc = 0.5*dt_test
1035 do itt_dt=1,dt_refinements
1037 nzc, dz, i_dz_int, dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1038 gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=cs%vel_underflow)
1040 idtt = 1.0 / (dt_test+dt_inc)
1041 do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1042 if (n2(k) < ri_crit * s2(k))
then
1043 k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1044 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1045 if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1046 (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))))
then
1047 valid_dt = .false. ;
exit
1050 if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))
then
1051 valid_dt = .false. ; k_src(k) = 0.0 ;
exit
1056 if (valid_dt) dt_test = dt_test + dt_inc
1063 dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc, dt_rem)
1065 local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
1072 if (ke_kappa < ks_kappa)
then
1075 dt_wt = dt_rem / dt ; dt_rem = 0.0
1080 tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
1082 tke_pred(k) = tke(k) ; kappa_pred(k) = 0.0 ; kappa(k) = 0.0
1089 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1090 gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1091 vel_underflow=cs%vel_underflow)
1095 do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ;
enddo
1096 call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1097 nzc, cs, gv, us, k_q_tmp, tke_pred, kappa_pred)
1100 ks_kappa = gv%ke+1 ; ke_kappa = 0
1102 kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1103 if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
1104 if (kappa_mid(k) > 0.0) ke_kappa = k
1109 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1110 gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1111 vel_underflow=cs%vel_underflow)
1115 call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1116 nzc, cs, gv, us, k_q, tke_pred, kappa_pred)
1120 dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1122 kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1123 kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
1124 tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
1125 kappa(k) = kappa_pred(k)
1130 if (dt_rem > 0.0)
then
1134 dz, i_dz_int, dbuoy_dt, dbuoy_ds, u, v, t, sal, &
1135 gv, us, n2, s2, vel_underflow=cs%vel_underflow)
1140 if (itt <= max_debug_itt)
then
1141 dt_it1(itt) = dt_now
1142 dk_wt_it1(itt) = 0.0 ; dkpos_wt_it1(itt) = 0.0 ; dkneg_wt_it1(itt) = 0.0
1144 wt_itt = 1.0/real(itt) ; wt_tot = 0.0
1146 ksrc_av(k) = (1.0-wt_itt)*ksrc_av(k) + wt_itt*k_src(k)
1147 wt_tot = wt_tot + dz_int(k) * ksrc_av(k)
1150 i_wt_tot = 0.0 ;
if (wt_tot > 0.0) i_wt_tot = 1.0/wt_tot
1153 wt(k) = (dz_int(k)*ksrc_av(k)) * i_wt_tot
1154 k_mag(itt) = k_mag(itt) + wt(k)*kappa_mid(k)
1155 dkappa_it1(k,itt) = kappa_pred(k) - kappa_out(k)
1156 dk_wt_it1(itt) = dk_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1157 if (dkappa_it1(k,itt) > 0.0)
then
1158 dkpos_wt_it1(itt) = dkpos_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1160 dkneg_wt_it1(itt) = dkneg_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1162 wt_it1(k,itt) = wt(k)
1166 ri_k(k) = 1e3 ;
if (n2(k) < 1e3 * s2(k)) ri_k(k) = n2(k) / s2(k)
1167 dtke(k) = tke_pred(k) - tke(k)
1168 dtke_norm(k) = dtke(k) / (0.5*(tke(k) + tke_pred(k)))
1169 dkappa(k) = kappa_pred(k) - kappa_out(k)
1171 if (itt <= max_debug_itt)
then
1173 u_it1(k,itt) = u(k) ; v_it1(k,itt) = v(k)
1174 t_it1(k,itt) = t(k) ; s_it1(k,itt) = sal(k)
1177 kprev_it1(k,itt) = kappa_out(k)
1178 kappa_it1(k,itt) = kappa_mid(k) ; tke_it1(k,itt) = 0.5*(tke(k)+tke_pred(k))
1179 n2_it1(k,itt)=n2(k) ; sh2_it1(k,itt)=s2(k)
1180 ksrc_it1(k,itt) = kappa_src(k)
1181 k_q_it1(k,itt) = kappa_out(k) / (tke(k))
1183 if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
1184 d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1186 dkappa_norm(k,itt) = dkappa(k) / max(0.5*(kappa_pred(k) + kappa_out(k)), us%m2_s_to_Z2_T*1e-100)
1191 if (dt_rem <= 0.0)
exit
1195 #ifdef ADD_DIAGNOSTICS
1196 if (
present(i_ld2_1d))
then
1197 do k=1,gv%ke+1 ; i_ld2_1d(k) = 0.0 ;
enddo
1198 do k=2,nzc ;
if (tke(k) > 0.0) &
1199 i_ld2_1d(k) = i_l2_bdry(k) + (n2(k) / cs%lambda**2 + f2) / tke(k)
1202 if (
present(dz_int_1d))
then
1203 do k=1,nzc+1 ; dz_int_1d(k) = dz_int(k) ;
enddo
1204 do k=nzc+2,gv%ke ; dz_int_1d(k) = 0.0 ;
enddo
1214 dz, I_dz_int, dbuoy_dT, dbuoy_dS, &
1215 u, v, T, Sal, GV, US, N2, S2, ks_int, ke_int, vel_underflow)
1216 integer,
intent(in) :: nz
1218 real,
dimension(nz+1),
intent(in) :: kappa
1220 real,
dimension(nz),
intent(in) :: u0
1221 real,
dimension(nz),
intent(in) :: v0
1222 real,
dimension(nz),
intent(in) :: T0
1223 real,
dimension(nz),
intent(in) :: S0
1224 real,
dimension(nz),
intent(in) :: dz
1225 real,
dimension(nz+1),
intent(in) :: I_dz_int
1227 real,
dimension(nz+1),
intent(in) :: dbuoy_dT
1229 real,
dimension(nz+1),
intent(in) :: dbuoy_dS
1231 real,
intent(in) :: dt
1232 real,
dimension(nz),
intent(inout) :: u
1233 real,
dimension(nz),
intent(inout) :: v
1234 real,
dimension(nz),
intent(inout) :: T
1235 real,
dimension(nz),
intent(inout) :: Sal
1238 real,
dimension(nz+1),
optional, &
1240 real,
dimension(nz+1),
optional, &
1242 integer,
optional,
intent(in) :: ks_int
1243 integer,
optional,
intent(in) :: ke_int
1245 real,
optional,
intent(in) :: vel_underflow
1250 real,
dimension(nz+1) :: c1
1253 real :: underflow_vel
1254 real :: a_a, a_b, b1, d1, bd1, b1nz_0
1255 integer :: k, ks, ke
1258 if (
present(ks_int)) ks = max(ks_int-1,1)
1259 if (
present(ke_int)) ke = min(ke_int,nz)
1260 underflow_vel = 0.0 ;
if (
present(vel_underflow)) underflow_vel = vel_underflow
1265 a_b = dt*(kappa(ks+1)*i_dz_int(ks+1))
1266 b1 = 1.0 / (dz(ks) + a_b)
1267 c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1
1269 u(ks) = (b1 * dz(ks))*u0(ks) ; v(ks) = (b1 * dz(ks))*v0(ks)
1270 t(ks) = (b1 * dz(ks))*t0(ks) ; sal(ks) = (b1 * dz(ks))*s0(ks)
1273 a_b = dt*(kappa(k+1)*i_dz_int(k+1))
1274 bd1 = dz(k) + d1*a_a
1275 b1 = 1.0 / (bd1 + a_b)
1276 c1(k+1) = a_b * b1 ; d1 = bd1 * b1
1278 u(k) = b1 * (dz(k)*u0(k) + a_a*u(k-1))
1279 v(k) = b1 * (dz(k)*v0(k) + a_a*v(k-1))
1280 t(k) = b1 * (dz(k)*t0(k) + a_a*t(k-1))
1281 sal(k) = b1 * (dz(k)*s0(k) + a_a*sal(k-1))
1288 b1 = 1.0 / (dz(ke) + d1*a_a)
1289 t(ke) = b1 * (dz(ke)*t0(ke) + a_a*t(ke-1))
1290 sal(ke) = b1 * (dz(ke)*s0(ke) + a_a*sal(ke-1))
1296 a_b = dt*(kappa(nz+1)*i_dz_int(nz+1))
1297 b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1301 u(ke) = b1nz_0 * (dz(ke)*u0(ke) + a_a*u(ke-1))
1302 v(ke) = b1nz_0 * (dz(ke)*v0(ke) + a_a*v(ke-1))
1303 if (abs(u(ke)) < underflow_vel) u(ke) = 0.0
1304 if (abs(v(ke)) < underflow_vel) v(ke) = 0.0
1307 u(k) = u(k) + c1(k+1)*u(k+1)
1308 v(k) = v(k) + c1(k+1)*v(k+1)
1309 if (abs(u(k)) < underflow_vel) u(k) = 0.0
1310 if (abs(v(k)) < underflow_vel) v(k) = 0.0
1311 t(k) = t(k) + c1(k+1)*t(k+1)
1312 sal(k) = sal(k) + c1(k+1)*sal(k+1)
1316 u(k) = u0(k) ; v(k) = v0(k) ; t(k) = t0(k) ; sal(k) = s0(k)
1317 if (abs(u(k)) < underflow_vel) u(k) = 0.0
1318 if (abs(v(k)) < underflow_vel) v(k) = 0.0
1322 if (
present(s2))
then
1324 l2_to_z2 = us%L_to_Z**2
1325 s2(1) = 0.0 ; s2(nz+1) = 0.0
1327 s2(ks) = ((u(ks)-u0(ks-1))**2 + (v(ks)-v0(ks-1))**2) * (l2_to_z2*i_dz_int(ks)**2)
1329 s2(k) = ((u(k)-u(k-1))**2 + (v(k)-v(k-1))**2) * (l2_to_z2*i_dz_int(k)**2)
1332 s2(ke+1) = ((u0(ke+1)-u(ke))**2 + (v0(ke+1)-v(ke))**2) * (l2_to_z2*i_dz_int(ke+1)**2)
1335 if (
present(n2))
then
1336 n2(1) = 0.0 ; n2(nz+1) = 0.0
1338 n2(ks) = max(0.0, i_dz_int(ks) * &
1339 (dbuoy_dt(ks) * (t0(ks-1)-t(ks)) + dbuoy_ds(ks) * (s0(ks-1)-sal(ks))))
1341 n2(k) = max(0.0, i_dz_int(k) * &
1342 (dbuoy_dt(k) * (t(k-1)-t(k)) + dbuoy_ds(k) * (sal(k-1)-sal(k))))
1345 n2(ke+1) = max(0.0, i_dz_int(ke+1) * &
1346 (dbuoy_dt(ke+1) * (t(ke)-t0(ke+1)) + dbuoy_ds(ke+1) * (sal(ke)-s0(ke+1))))
1352 subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, &
1353 nz, CS, GV, US, K_Q, tke, kappa, kappa_src, local_src)
1354 integer,
intent(in) :: nz
1355 real,
dimension(nz+1),
intent(in) :: N2
1356 real,
dimension(nz+1),
intent(in) :: S2
1357 real,
dimension(nz+1),
intent(in) :: kappa_in
1359 real,
dimension(nz+1),
intent(in) :: dz_Int
1361 real,
dimension(nz+1),
intent(in) :: I_L2_bdry
1363 real,
dimension(nz),
intent(in) :: Idz
1364 real,
intent(in) :: f2
1368 real,
dimension(nz+1),
intent(inout) :: K_Q
1371 real,
dimension(nz+1),
intent(out) :: tke
1373 real,
dimension(nz+1),
intent(out) :: kappa
1375 real,
dimension(nz+1),
optional, &
1376 intent(out) :: kappa_src
1377 real,
dimension(nz+1),
optional, &
1378 intent(out) :: local_src
1383 real,
dimension(nz) :: &
1386 real,
dimension(nz+1) :: &
1407 real :: cQcomp, cKcomp
1421 real :: eden1, eden2, I_eden, ome
1422 real :: diffusive_src
1430 real :: decay_term_k
1431 real :: decay_term_Q
1441 real,
parameter :: roundoff = 1.0e-16
1444 logical :: tke_noflux_bottom_BC = .false.
1445 logical :: tke_noflux_top_BC = .false.
1446 logical :: do_Newton
1447 logical :: abort_Newton
1449 logical :: was_Newton
1450 logical :: within_tolerance
1452 integer :: ks_src, ke_src
1453 integer :: ks_kappa, ke_kappa, ke_tke
1454 integer :: ks_kappa_prev, ke_kappa_prev
1455 integer :: itt, k, k2
1457 integer :: max_debug_itt ;
parameter(max_debug_itt=20)
1458 real :: K_err_lin, Q_err_lin, TKE_src_norm
1459 real,
dimension(nz+1) :: &
1463 real,
dimension(nz+1,1:max_debug_itt) :: &
1464 tke_it1, kappa_it1, kprev_it1, &
1465 dkappa_it1, K_Q_it1, d_dkappa_it1, dkappa_norm_it1
1469 c_n2 = cs%C_N**2 ; c_s2 = cs%C_S**2
1470 q0 = cs%TKE_bg ; kappa0 = cs%kappa_0
1471 tke_min = max(cs%TKE_bg, 1.0e-20*us%m_to_Z**2*us%T_to_s**2)
1472 ri_crit = cs%Rino_crit
1473 ilambda2 = 1.0 / cs%lambda**2
1474 kappa_trunc = cs%kappa_trunc
1475 do_newton = .false. ; abort_newton = .false.
1476 tol_err = cs%kappa_tol_err
1479 ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1481 ke_src = 0 ; ks_src = nz+1
1483 if (n2(k) < ri_crit * s2(k))
then
1487 k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1488 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1490 if (ks_src > k) ks_src = k
1497 if (ks_src > ke_src)
then
1499 kappa(k) = 0.0 ; k_q(k) = 0.0 ; tke(k) = tke_min
1501 if (
present(kappa_src))
then ;
do k=1,nz+1 ; kappa_src(k) = 0.0 ;
enddo ;
endif
1502 if (
present(local_src))
then ;
do k=1,nz+1 ; local_src(k) = 0.0 ;
enddo ;
endif
1507 kappa(k) = kappa_in(k)
1509 tke_decay(k) = sqrt(c_n2*n2(k) + c_s2*s2(k))
1510 if ((kappa(k) > 0.0) .and. (k_q(k) > 0.0))
then
1511 tke(k) = kappa(k) / k_q(k)
1517 kappa(1) = 0.0 ; kappa(nz+1) = 0.0
1522 eden2 = kappa0 * idz(nz)
1523 if (tke_noflux_bottom_bc)
then
1524 eden1 = dz_int(nz+1)*tke_decay(nz+1)
1525 i_eden = 1.0 / (eden2 + eden1)
1526 e1(nz+1) = eden2 * i_eden ; ome = eden1 * i_eden
1528 e1(nz+1) = 0.0 ; ome = 1.0
1531 eden1 = dz_int(k)*tke_decay(k) + ome * eden2
1532 eden2 = kappa0 * idz(k-1)
1533 i_eden = 1.0 / (eden2 + eden1)
1534 e1(k) = eden2 * i_eden ; ome = eden1 * i_eden
1540 do itt=1,cs%max_RiNo_it
1547 do k=1,nz+1 ; kappa_prev(k) = kappa(k) ; tke_prev(k) = tke(k) ;
enddo
1550 if (.not.do_newton)
then
1556 ke_tke = max(ke_kappa,ke_kappa_prev)+1
1558 do k=1,min(ke_tke,nz)
1559 aq(k) = (0.5*(kappa(k)+kappa(k+1)) + kappa0) * idz(k)
1562 if (tke_noflux_top_bc)
then
1563 tke_src = kappa0*s2(1) + q0 * tke_decay(1)
1564 bqd1 = dz_int(1) * tke_decay(1)
1565 bq = 1.0 / (bqd1 + aq(1))
1566 tke(1) = bq * (dz_int(1)*tke_src)
1567 cq(2) = aq(1) * bq ; cqcomp = bqd1 * bq
1569 tke(1) = q0 ; cq(2) = 0.0 ; cqcomp = 1.0
1573 tke_src = (kappa(k) + kappa0)*s2(k) + q0*tke_decay(k)
1574 bqd1 = dz_int(k)*(tke_decay(k) + n2(k)*k_q(k)) + cqcomp*aq(k-1)
1575 bq = 1.0 / (bqd1 + aq(k))
1576 tke(k) = bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1))
1577 cq(k+1) = aq(k) * bq ; cqcomp = bqd1 * bq
1579 if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_bc))
then
1584 tke_src = kappa0*s2(k) + q0*tke_decay(k)
1587 bq = 1.0 / (dz_int(k)*tke_decay(k) + cqcomp*aq(k-1))
1588 tke(k) = max(tke_min, bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)))
1589 dq(k) = tke(k) + dq(k)
1591 bq = 1.0 / ((dz_int(k)*tke_decay(k) + cqcomp*aq(k-1)) + aq(k))
1592 cq(k+1) = aq(k) * bq
1595 tke(k) = max((bq * (dz_int(k)*tke_src + aq(k-1)*tke(k-1)) + &
1596 cq(k+1)*(tke(k+1) - e1(k+1)*tke(k))) / (1.0 - cq(k+1)*e1(k+1)), tke_min)
1597 dq(k) = tke(k) + dq(k)
1602 dq(k) = e1(k)*dq(k-1)
1603 tke(k) = max(tke(k) + dq(k), tke_min)
1604 if (abs(dq(k)) < roundoff*tke(k))
exit
1607 if (dq(k2) == 0.0)
exit
1613 tke(k) = max(tke(k) + cq(k+1)*tke(k+1), tke_min)
1614 dq(k) = tke(k) + dq(k)
1621 ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1624 ck(2) = 0.0 ; ckcomp = 1.0
1625 if (itt == 1)
then ;
dO k=2,nz
1626 i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1631 i_ld2(k) = (n2(k)*ilambda2 + f2) / tke(k) + i_l2_bdry(k)
1632 bkd1 = dz_int(k)*i_ld2(k) + ckcomp*idz(k-1)
1633 bk = 1.0 / (bkd1 + idz(k))
1635 kappa(k) = bk * (idz(k-1)*kappa(k-1) + dz_int(k) * k_src(k))
1636 ck(k+1) = idz(k) * bk ; ckcomp = bkd1 * bk
1639 if (kappa(k) < ckcomp*kappa_trunc)
then
1641 if (k > ke_src)
then ; ke_kappa = k-1 ; k_q(k) = 0.0 ;
exit ;
endif
1642 elseif (kappa(k) < 2.0*ckcomp*kappa_trunc)
then
1643 kappa(k) = 2.0 * (kappa(k) - ckcomp*kappa_trunc)
1646 k_q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1647 dk(ke_kappa) = dk(ke_kappa) + kappa(ke_kappa)
1648 do k=ke_kappa+2,ke_kappa_prev
1649 dk(k) = -kappa(k) ; kappa(k) = 0.0 ; k_q(k) = 0.0
1651 do k=ke_kappa-1,2,-1
1652 kappa(k) = kappa(k) + ck(k+1)*kappa(k+1)
1654 if (kappa(k) <= kappa_trunc)
then
1656 if (k < ks_src)
then ; ks_kappa = k+1 ; k_q(k) = 0.0 ;
exit ;
endif
1657 elseif (kappa(k) < 2.0*kappa_trunc)
then
1658 kappa(k) = 2.0 * (kappa(k) - kappa_trunc)
1661 dk(k) = dk(k) + kappa(k)
1662 k_q(k) = kappa(k) / tke(k)
1664 do k=ks_kappa_prev,ks_kappa-2 ; kappa(k) = 0.0 ; k_q(k) = 0.0 ;
enddo
1669 ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1671 dk(1) = 0.0 ; ck(2) = 0.0 ; ckcomp = 1.0 ; dkdq(1) = 0.0
1672 aq(1) = (0.5*(kappa(1)+kappa(2))+kappa0) * idz(1)
1673 dqdz(1) = 0.5*(tke(1) - tke(2))*idz(1)
1674 if (tke_noflux_top_bc)
then
1675 tke_src = dz_int(1) * (kappa0*s2(1) - (tke(1) - q0)*tke_decay(1)) - &
1676 aq(1) * (tke(1) - tke(2))
1678 bq = 1.0 / (aq(1) + dz_int(1)*tke_decay(1))
1680 cqcomp = (dz_int(1)*tke_decay(1)) * bq
1681 dqmdk(2) = -dqdz(1) * bq
1682 dq(1) = bq * tke_src
1684 dq(1) = 0.0 ; cq(2) = 0.0 ; cqcomp = 1.0 ; dqmdk(2) = 0.0
1688 i_ld2(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1690 kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa(k)) + &
1691 idz(k-1)*(kappa(k-1)-kappa(k)) - idz(k)*(kappa(k)-kappa(k+1))
1695 decay_term_k = -idz(k-1)*dqmdk(k)*dkdq(k-1) + dz_int(k)*i_ld2(k)
1696 if (decay_term_k < 0.0)
then ; abort_newton = .true. ;
exit ;
endif
1697 bk = 1.0 / (idz(k) + idz(k-1)*ckcomp + decay_term_k)
1699 ck(k+1) = bk * idz(k)
1700 ckcomp = bk * (idz(k-1)*ckcomp + decay_term_k)
1701 if (cs%dKdQ_iteration_bug)
then
1702 dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1703 us%m_to_Z*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1705 dkdq(k) = bk * (idz(k-1)*dkdq(k-1)*cq(k) + &
1706 dz_int(k)*(n2(k)*ilambda2 + f2) * i_q**2 * kappa(k) )
1708 dk(k) = bk * (kap_src + idz(k-1)*dk(k-1) + idz(k-1)*dkdq(k-1)*dq(k-1))
1711 if (dk(k) <= ckcomp*(kappa_trunc - kappa(k)))
then
1712 dk(k) = -ckcomp*kappa(k)
1714 elseif (dk(k) < ckcomp*(2.0*kappa_trunc - kappa(k)))
then
1715 dk(k) = 2.0 * dk(k) - ckcomp*(2.0*kappa_trunc - kappa(k))
1719 aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1720 dqdz(k) = 0.5*(tke(k) - tke(k+1))*idz(k)
1721 tke_src = dz_int(k) * (((kappa(k) + kappa0)*s2(k) - kappa(k)*n2(k)) - &
1722 (tke(k) - q0)*tke_decay(k)) - &
1723 (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k)))
1724 v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1725 v2 = (v1*dqmdk(k) + dqdz(k-1)*ck(k)) + &
1726 ((dqdz(k-1) - dqdz(k)) + dz_int(k)*(s2(k) - n2(k)))
1730 decay_term_q = dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k) - v2*dkdq(k)
1731 if (decay_term_q < 0.0)
then ; abort_newton = .true. ;
exit ;
endif
1732 bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1734 cq(k+1) = aq(k) * bq
1735 cqcomp = (cqcomp*aq(k-1) + decay_term_q) * bq
1736 dqmdk(k+1) = (v2 * ck(k+1) - dqdz(k)) * bq
1739 dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + &
1740 (v2 * dk(k) + tke_src)), cqcomp*(-0.5*tke(k)))
1743 if ((itt > 1) .and. (k > ke_src) .and. (dk(k) == 0.0) .and. &
1744 ((kappa(k) + kappa(k+1)) == 0.0))
then
1746 ke_kappa = k-1 ;
exit
1749 if ((ke_kappa == nz) .and. (.not. abort_newton))
then
1750 dk(nz+1) = 0.0 ; dkdq(nz+1) = 0.0
1751 if (tke_noflux_bottom_bc)
then
1753 tke_src = dz_int(k) * (kappa0*s2(k) - (tke(k) - q0)*tke_decay(k)) + &
1754 aq(k-1) * (tke(k-1) - tke(k))
1756 v1 = aq(k-1) + dqdz(k-1)*dkdq(k-1)
1757 decay_term_q = max(0.0, dz_int(k)*tke_decay(k) - dqdz(k-1)*dkdq(k-1)*cq(k))
1758 if (decay_term_q < 0.0)
then
1759 abort_newton = .true.
1761 bq = 1.0 / (aq(k) + (cqcomp*aq(k-1) + decay_term_q))
1763 dq(k) = max(bq * ((v1 * dq(k-1) + dqdz(k-1)*dk(k-1)) + tke_src), -0.5*tke(k))
1764 tke(k) = max(tke(k) + dq(k), tke_min)
1769 elseif (.not. abort_newton)
then
1771 dq(ke_kappa+1) = dq(ke_kappa+1) / (1.0 - cq(ke_kappa+2)*e1(ke_kappa+2))
1772 tke(ke_kappa+1) = max(tke(ke_kappa+1) + dq(ke_kappa+1), tke_min)
1773 do k=ke_kappa+2,nz+1
1777 aq(k) = (0.5*(kappa(k)+kappa(k+1))+kappa0) * idz(k)
1778 tke_src_norm = (dz_int(k) * (kappa0*s2(k) - (tke(k)-q0)*tke_decay(k)) - &
1779 (aq(k) * (tke(k) - tke(k+1)) - aq(k-1) * (tke(k-1) - tke(k))) ) / &
1780 (aq(k) + (aq(k-1) + dz_int(k)*tke_decay(k)))
1785 dq(k) = max(e1(k)*dq(k-1),-0.5*tke(k))
1786 tke(k) = max(tke(k) + dq(k), tke_min)
1787 if (abs(dq(k)) < roundoff*tke(k))
exit
1790 do k2=k+1,ke_kappa_prev+1 ; dq(k2) = 0.0 ; dk(k2) = 0.0 ;
enddo
1791 do k=k2,nz+1 ;
if (dq(k) == 0.0)
exit ; dq(k) = 0.0 ; dk(k) = 0.0 ;
enddo
1794 if (.not. abort_newton)
then
1797 dq(k) = max(dq(k) + (cq(k+1)*dq(k+1) + dqmdk(k+1) * dk(k+1)), -0.5*tke(k))
1798 tke(k) = max(tke(k) + dq(k), tke_min)
1799 dk(k) = dk(k) + (ck(k+1)*dk(k+1) + dkdq(k) * dq(k))
1801 if (dk(k) <= kappa_trunc - kappa(k))
then
1804 if ((k < ks_src) .and. (k+1 > ks_kappa)) ks_kappa = k+1
1805 elseif (dk(k) < 2.0*kappa_trunc - kappa(k))
then
1806 dk(k) = 2.0*dk(k) - (2.0*kappa_trunc - kappa(k))
1807 kappa(k) = max(kappa(k) + dk(k), 0.0)
1808 if (k<=ks_kappa) ks_kappa = 2
1810 kappa(k) = kappa(k) + dk(k)
1811 if (k<=ks_kappa) ks_kappa = 2
1814 dq(1) = max(dq(1) + cq(2)*dq(2) + dqmdk(2) * dk(2), tke_min - tke(1))
1815 tke(1) = max(tke(1) + dq(1), tke_min)
1829 i_ld2_debug(k) = (n2(k)*ilambda2 + f2) * i_q + i_l2_bdry(k)
1831 kap_src = dz_int(k) * (k_src(k) - i_ld2(k)*kappa_prev(k)) + &
1832 (idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
1833 idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
1834 k_err_lin = -idz(k-1)*(dk(k-1)-dk(k)) + idz(k)*(dk(k)-dk(k+1)) + &
1835 dz_int(k)*i_ld2_debug(k)*dk(k) - kap_src - &
1836 dz_int(k)*(n2(k)*ilambda2 + f2)*i_q**2*kappa_prev(k) * dq(k)
1838 tke_src = dz_int(k) * ((kappa_prev(k) + kappa0)*s2(k) - &
1839 kappa_prev(k)*n2(k) - (tke_prev(k) - q0)*tke_decay(k)) - &
1840 (aq(k) * (tke_prev(k) - tke_prev(k+1)) - aq(k-1) * (tke_prev(k-1) - tke_prev(k)))
1841 q_err_lin = tke_src + (aq(k-1) * (dq(k-1)-dq(k)) - aq(k) * (dq(k)-dq(k+1))) - &
1842 0.5*(tke_prev(k)-tke_prev(k+1))*idz(k) * (dk(k) + dk(k+1)) - &
1843 0.5*(tke_prev(k)-tke_prev(k-1))*idz(k-1)* (dk(k-1) + dk(k)) + &
1844 dz_int(k) * (dk(k) * (s2(k) - n2(k)) - dq(k)*tke_decay(k))
1850 if ((tol_err < newton_err) .and. (.not.abort_newton))
then
1852 newton_test = newton_err ;
if (do_newton) newton_test = 2.0*newton_err
1853 was_newton = do_newton
1854 within_tolerance = .true. ; do_newton = .true.
1855 do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1856 kappa_mean = kappa0 + (kappa(k) - 0.5*dk(k))
1857 if (abs(dk(k)) > newton_test * kappa_mean)
then
1858 if (do_newton) abort_newton = .true.
1859 within_tolerance = .false. ; do_newton = .false. ;
exit
1860 elseif (abs(dk(k)) > tol_err * kappa_mean)
then
1861 within_tolerance = .false. ;
if (.not.do_newton)
exit
1863 if (abs(dq(k)) > newton_test*(tke(k) - 0.5*dq(k)))
then
1864 if (do_newton) abort_newton = .true.
1865 do_newton = .false. ;
if (.not.within_tolerance)
exit
1870 within_tolerance = .true.
1871 do k=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
1872 if (abs(dk(k)) > tol_err * (kappa0 + (kappa(k) - 0.5*dk(k))))
then
1873 within_tolerance = .false. ;
exit
1878 if (abort_newton)
then
1879 do_newton = .false. ; abort_newton = .false.
1881 newton_err = 0.5*newton_err
1883 do k=2,nz ; k_q(k) = kappa(k) / max(tke(k), tke_min) ;
enddo
1887 if (itt <= max_debug_itt)
then
1889 kprev_it1(k,itt) = kappa_prev(k)
1890 kappa_it1(k,itt) = kappa(k) ; tke_it1(k,itt) = tke(k)
1891 dkappa_it1(k,itt) = kappa(k) - kappa_prev(k)
1892 dkappa_norm_it1(k,itt) = (kappa(k) - kappa_prev(k)) / &
1893 (kappa0 + 0.5*(kappa(k) + kappa_prev(k)))
1894 k_q_it1(k,itt) = kappa(k) / max(tke(k),tke_min)
1895 d_dkappa_it1(k,itt) = 0.0
1896 if (itt > 1)
then ;
if (abs(dkappa_it1(k,itt-1)) > 1e-20*us%m2_s_to_Z2_T) &
1897 d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1903 if (within_tolerance)
exit
1908 do it2=itt+1,max_debug_itt ;
do k=1,nz+1
1909 kprev_it1(k,it2) = 0.0 ; kappa_it1(k,it2) = 0.0 ; tke_it1(k,it2) = 0.0
1910 dkappa_it1(k,it2) = 0.0 ; k_q_it1(k,it2) = 0.0 ; d_dkappa_it1(k,it2) = 0.0
1915 do k=1,ks_kappa-1 ; k_q(k) = 0.0 ;
enddo
1916 do k=ks_kappa,ke_kappa ; k_q(k) = kappa(k) / tke(k) ;
enddo
1917 do k=ke_kappa+1,nz+1 ; k_q(k) = 0.0 ;
enddo
1920 if (
present(local_src))
then
1921 local_src(1) = 0.0 ; local_src(nz+1) = 0.0
1923 diffusive_src = idz(k-1)*(kappa(k-1)-kappa(k)) + idz(k)*(kappa(k+1)-kappa(k))
1924 chg_by_k0 = kappa0 * ((idz(k-1)+idz(k)) / dz_int(k) + i_ld2(k))
1925 if (diffusive_src <= 0.0)
then
1926 local_src(k) = k_src(k) + chg_by_k0
1928 local_src(k) = (k_src(k) + chg_by_k0) + diffusive_src / dz_int(k)
1932 if (
present(kappa_src))
then
1933 kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
1935 kappa_src(k) = k_src(k)
1943 type(time_type),
intent(in) :: time
1949 type(
diag_ctrl),
target,
intent(inout) :: diag
1956 logical :: merge_mixedlayer
1958 #include "version_variable.h"
1959 character(len=40) :: mdl =
"MOM_kappa_shear"
1960 real :: kappa_0_unscaled
1964 if (
associated(cs))
then
1965 call mom_error(warning,
"kappa_shear_init called with an associated "// &
1966 "control structure.")
1982 "Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008")
1984 "If true, use the Jackson-Hallberg-Legg (JPO 2008) "//&
1985 "shear mixing parameterization.", default=.false.)
1986 call get_param(param_file, mdl,
"VERTEX_SHEAR", cs%KS_at_vertex, &
1987 "If true, do the calculations of the shear-driven mixing "//&
1988 "at the cell vertices (i.e., the vorticity points).", &
1990 call get_param(param_file, mdl,
"RINO_CRIT", cs%RiNo_crit, &
1991 "The critical Richardson number for shear mixing.", &
1992 units=
"nondim", default=0.25)
1993 call get_param(param_file, mdl,
"SHEARMIX_RATE", cs%Shearmix_rate, &
1994 "A nondimensional rate scale for shear-driven entrainment. "//&
1995 "Jackson et al find values in the range of 0.085-0.089.", &
1996 units=
"nondim", default=0.089)
1997 call get_param(param_file, mdl,
"MAX_RINO_IT", cs%max_RiNo_it, &
1998 "The maximum number of iterations that may be used to "//&
1999 "estimate the Richardson number driven mixing.", &
2000 units=
"nondim", default=50)
2001 call get_param(param_file, mdl,
"KD", kd_normal, default=1.0e-7, do_not_log=.true.)
2002 call get_param(param_file, mdl,
"KD_KAPPA_SHEAR_0", cs%kappa_0, &
2003 "The background diffusivity that is used to smooth the "//&
2004 "density and shear profiles before solving for the "//&
2005 "diffusivities. Defaults to value of KD.", &
2006 units=
"m2 s-1", default=kd_normal, scale=us%m2_s_to_Z2_T, unscaled=kappa_0_unscaled)
2007 call get_param(param_file, mdl,
"KD_TRUNC_KAPPA_SHEAR", cs%kappa_trunc, &
2008 "The value of shear-driven diffusivity that is considered negligible "//&
2009 "and is rounded down to 0. The default is 1% of KD_KAPPA_SHEAR_0.", &
2010 units=
"m2 s-1", default=0.01*kappa_0_unscaled, scale=us%m2_s_to_Z2_T)
2011 call get_param(param_file, mdl,
"FRI_CURVATURE", cs%FRi_curvature, &
2012 "The nondimensional curvature of the function of the "//&
2013 "Richardson number in the kappa source term in the "//&
2014 "Jackson et al. scheme.", units=
"nondim", default=-0.97)
2015 call get_param(param_file, mdl,
"TKE_N_DECAY_CONST", cs%C_N, &
2016 "The coefficient for the decay of TKE due to "//&
2017 "stratification (i.e. proportional to N*tke). "//&
2018 "The values found by Jackson et al. are 0.24-0.28.", &
2019 units=
"nondim", default=0.24)
2022 call get_param(param_file, mdl,
"TKE_SHEAR_DECAY_CONST", cs%C_S, &
2023 "The coefficient for the decay of TKE due to shear (i.e. "//&
2024 "proportional to |S|*tke). The values found by Jackson "//&
2025 "et al. are 0.14-0.12.", units=
"nondim", default=0.14)
2026 call get_param(param_file, mdl,
"KAPPA_BUOY_SCALE_COEF", cs%lambda, &
2027 "The coefficient for the buoyancy length scale in the "//&
2028 "kappa equation. The values found by Jackson et al. are "//&
2029 "in the range of 0.81-0.86.", units=
"nondim", default=0.82)
2030 call get_param(param_file, mdl,
"KAPPA_N_OVER_S_SCALE_COEF2", cs%lambda2_N_S, &
2031 "The square of the ratio of the coefficients of the "//&
2032 "buoyancy and shear scales in the diffusivity equation, "//&
2033 "Set this to 0 (the default) to eliminate the shear scale. "//&
2034 "This is only used if USE_JACKSON_PARAM is true.", &
2035 units=
"nondim", default=0.0)
2036 call get_param(param_file, mdl,
"KAPPA_SHEAR_TOL_ERR", cs%kappa_tol_err, &
2037 "The fractional error in kappa that is tolerated. "//&
2038 "Iteration stops when changes between subsequent "//&
2039 "iterations are smaller than this everywhere in a "//&
2040 "column. The peak diffusivities usually converge most "//&
2041 "rapidly, and have much smaller errors than this.", &
2042 units=
"nondim", default=0.1)
2043 call get_param(param_file, mdl,
"TKE_BACKGROUND", cs%TKE_bg, &
2044 "A background level of TKE used in the first iteration "//&
2045 "of the kappa equation. TKE_BACKGROUND could be 0.", &
2046 units=
"m2 s-2", default=0.0, scale=us%m_to_Z**2*us%T_to_s**2)
2047 call get_param(param_file, mdl,
"KAPPA_SHEAR_ELIM_MASSLESS", cs%eliminate_massless, &
2048 "If true, massless layers are merged with neighboring "//&
2049 "massive layers in this calculation. The default is "//&
2050 "true and I can think of no good reason why it should "//&
2051 "be false. This is only used if USE_JACKSON_PARAM is true.", &
2053 call get_param(param_file, mdl,
"MAX_KAPPA_SHEAR_IT", cs%max_KS_it, &
2054 "The maximum number of iterations that may be used to "//&
2055 "estimate the time-averaged diffusivity.", units=
"nondim", &
2057 call get_param(param_file, mdl,
"PRANDTL_TURB", cs%Prandtl_turb, &
2058 "The turbulent Prandtl number applied to shear "//&
2059 "instability.", units=
"nondim", default=1.0, do_not_log=.true.)
2060 call get_param(param_file, mdl,
"VEL_UNDERFLOW", cs%vel_underflow, &
2061 "A negligibly small velocity magnitude below which velocity "//&
2062 "components are set to 0. A reasonable value might be "//&
2063 "1e-30 m/s, which is less than an Angstrom divided by "//&
2064 "the age of the universe.", units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
2066 call get_param(param_file, mdl,
"DEBUG_KAPPA_SHEAR", cs%debug, &
2067 "If true, write debugging data for the kappa-shear code. \n"//&
2068 "Caution: this option is _very_ verbose and should only "//&
2069 "be used in single-column mode!", &
2070 default=.false., debuggingparam=.true.)
2071 call get_param(param_file, mdl,
"KAPPA_SHEAR_ITER_BUG", cs%dKdQ_iteration_bug, &
2072 "If true. use an older, dimensionally inconsistent estimate of the "//&
2073 "derivative of diffusivity with energy in the Newton's method iteration. "//&
2074 "The bug causes undercorrections when dz > 1m.", default=.true.)
2082 call get_param(param_file, mdl,
"KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
2083 "If true, combine the mixed layers together before "//&
2084 "solving the kappa-shear equations.", default=.true.)
2085 if (merge_mixedlayer) cs%nkml = gv%nkml
2093 cs%id_Kd_shear = register_diag_field(
'ocean_model',
'Kd_shear',diag%axesTi,time, &
2094 'Shear-driven Diapycnal Diffusivity',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
2095 cs%id_TKE = register_diag_field(
'ocean_model',
'TKE_shear',diag%axesTi,time, &
2096 'Shear-driven Turbulent Kinetic Energy',
'm2 s-2', conversion=us%Z_to_m**2*us%s_to_T**2)
2097 #ifdef ADD_DIAGNOSTICS
2098 cs%id_ILd2 = register_diag_field(
'ocean_model',
'ILd2_shear',diag%axesTi,time, &
2099 'Inverse kappa decay scale at interfaces',
'm-2', conversion=us%m_to_Z**2)
2100 cs%id_dz_Int = register_diag_field(
'ocean_model',
'dz_Int_shear',diag%axesTi,time, &
2101 'Finite volume thickness of interfaces',
'm', conversion=us%Z_to_m)
2111 character(len=40) :: mdl =
"MOM_kappa_shear"
2114 default=.false., do_not_log=.true.)
2122 character(len=40) :: mdl =
"MOM_kappa_shear"
2124 logical :: do_kappa_shear
2126 call get_param(param_file, mdl,
"USE_JACKSON_PARAM", do_kappa_shear, &
2127 default=.false., do_not_log=.true.)
2129 if (do_kappa_shear) &
2131 "If true, do the calculations of the shear-driven mixing "//&
2132 "at the cell vertices (i.e., the vorticity points).", &
2133 default=.false., do_not_log=.true.)