8 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
22 use mom_restart,
only : register_restart_field_as_obsolete
31 implicit none ;
private
33 #include <MOM_memory.h>
35 public set_viscous_bbl, set_viscous_ml, set_visc_init, set_visc_end
36 public set_visc_register_restarts
56 real :: htbl_shelf_min
59 logical :: bottomdraglaw
64 logical :: bbl_use_eos
66 logical :: linear_drag
67 logical :: channel_drag
71 logical :: dynamic_viscous_ml
84 logical :: answers_2018
92 integer :: id_bbl_thick_u = -1, id_kv_bbl_u = -1
93 integer :: id_bbl_thick_v = -1, id_kv_bbl_v = -1
94 integer :: id_ray_u = -1, id_ray_v = -1
95 integer :: id_nkml_visc_u = -1, id_nkml_visc_v = -1
109 subroutine set_viscous_bbl(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
113 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
115 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
117 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
126 logical,
optional,
intent(in) :: symmetrize
131 real,
dimension(SZIB_(G)) :: &
147 real,
dimension(SZIB_(G),SZJ_(G)) :: &
151 real,
dimension(SZI_(G),SZJB_(G)) :: &
155 real,
dimension(SZIB_(G),SZK_(G)) :: &
199 real :: v_at_u, u_at_v
203 real,
dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
205 real :: p_ref(szi_(g))
216 real :: apb_4a, ax2_3apb
217 real :: a2x48_apb3, iapb, ibma_2
257 real :: bbl_visc_frac
259 real,
parameter :: c1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0
262 real :: tmp_val_m1_to_p1
263 logical :: use_bbl_eos, do_i(szib_(g))
264 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, m, n, k2, nkmb, nkml
265 integer :: itt, maxitt=20
268 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
269 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
270 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
271 h_neglect = gv%H_subroundoff
272 rho0x400_g = 400.0*(gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
273 vol_quit = 0.9*gv%Angstrom_H + h_neglect
274 c2pi_3 = 8.0*atan(1.0)/3.0
276 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_vert_friction(BBL): "//&
277 "Module must be initialized before it is used.")
278 if (.not.cs%bottomdraglaw)
return
280 if (
present(symmetrize))
then ;
if (symmetrize)
then
281 jsq = js-1 ; isq = is-1
285 call uvchksum(
"Start set_viscous_BBL [uv]", u, v, g%HI, haloshift=1)
286 call hchksum(h,
"Start set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
287 if (
associated(tv%T))
call hchksum(tv%T,
"Start set_viscous_BBL T", g%HI, haloshift=1)
288 if (
associated(tv%S))
call hchksum(tv%S,
"Start set_viscous_BBL S", g%HI, haloshift=1)
291 use_bbl_eos =
associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
294 u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
295 cdrag_sqrt = sqrt(cs%cdrag)
296 cdrag_sqrt_z = us%m_to_Z * sqrt(cs%cdrag)
302 if ((nkml>0) .and. .not.use_bbl_eos)
then
303 do i=isq,ieq+1 ; p_ref(i) = tv%P_ref ;
enddo
305 do j=jsq,jeq+1 ;
do k=1,nkmb
307 rml(:,j,k), isq, ieq-isq+2, tv%eqn_of_state)
312 do j=js-1,je ;
do i=is-1,ie+1
313 d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
314 mask_v(i,j) = g%mask2dCv(i,j)
317 do j=js-1,je+1 ;
do i=is-1,ie
318 d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
319 mask_u(i,j) = g%mask2dCu(i,j)
322 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
323 if (.not. obc%segment(n)%on_pe) cycle
325 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
326 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je))
then
327 do i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
328 if (obc%segment(n)%direction == obc_direction_n) d_v(i,j) = g%bathyT(i,j)
329 if (obc%segment(n)%direction == obc_direction_s) d_v(i,j) = g%bathyT(i,j+1)
331 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie))
then
332 do j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
333 if (obc%segment(n)%direction == obc_direction_e) d_u(i,j) = g%bathyT(i,j)
334 if (obc%segment(n)%direction == obc_direction_w) d_u(i,j) = g%bathyT(i+1,j)
338 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
341 if (.not. obc%segment(n)%on_pe) cycle
343 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
344 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je))
then
345 do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
346 if (obc%segment(n)%direction == obc_direction_n)
then
347 d_u(i,j+1) = d_u(i,j) ; mask_u(i,j+1) = 0.0
348 elseif (obc%segment(n)%direction == obc_direction_s)
then
349 d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
352 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie))
then
353 do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
354 if (obc%segment(n)%direction == obc_direction_e)
then
355 d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
356 elseif (obc%segment(n)%direction == obc_direction_w)
then
357 d_v(i,j) = d_v(i+1,j) ; mask_v(i,j) = 0.0
363 if (.not.use_bbl_eos) rml_vel(:,:) = 0.0
369 do j=jsq,jeq ;
do m=1,2
376 if (g%mask2dCu(i,j) > 0) do_i(i) = .true.
379 is = g%isc ; ie = g%iec
382 if (g%mask2dCv(i,j) > 0) do_i(i) = .true.
387 do k=1,nz ;
do i=is,ie
389 if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0)
then
390 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
391 (h(i,j,k) + h(i+1,j,k) + h_neglect)
393 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
396 h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
398 if (use_bbl_eos)
then ;
do k=1,nz ;
do i=is,ie
400 t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
401 s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
402 enddo ;
enddo ;
else ;
do k=1,nkmb ;
do i=is,ie
403 rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
404 enddo ;
enddo ;
endif
406 do k=1,nz ;
do i=is,ie
408 if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0)
then
409 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
410 (h(i,j,k) + h(i,j+1,k) + h_neglect)
412 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
415 h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
417 if (use_bbl_eos)
then ;
do k=1,nz ;
do i=is,ie
418 t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
419 s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
420 enddo ;
enddo ;
else ;
do k=1,nkmb ;
do i=is,ie
421 rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
422 enddo ;
enddo ;
endif
425 if (
associated(obc))
then ;
if (obc%number_of_segments > 0)
then
428 do i=is,ie ;
if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none))
then
429 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
431 h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
433 if (use_bbl_eos)
then
435 t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
439 rml_vel(i,k) = rml(i,j,k)
442 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
444 h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
446 if (use_bbl_eos)
then
448 t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
452 rml_vel(i,k) = rml(i+1,j,k)
458 do i=is,ie ;
if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none))
then
459 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
461 h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
463 if (use_bbl_eos)
then
465 t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
469 rml_vel(i,k) = rml(i,j,k)
472 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
474 h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
476 if (use_bbl_eos)
then
478 t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
482 rml_vel(i,k) = rml(i,j+1,k)
490 if (use_bbl_eos .or. .not.cs%linear_drag)
then
491 do i=is,ie ;
if (do_i(i))
then
495 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
496 thtot = 0.0 ; shtot = 0.0
499 if (htot_vel>=cs%Hbbl)
exit
501 hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
502 if (hweight < 1.5*gv%Angstrom_H + h_neglect) cycle
504 htot_vel = htot_vel + h_at_vel(i,k)
505 hwtot = hwtot + hweight
507 if ((.not.cs%linear_drag) .and. (hweight >= 0.0))
then ;
if (m==1)
then
508 v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
509 hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
510 v_at_u*v_at_u + u_bg_sq)
512 u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
513 hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
514 u_at_v*u_at_v + u_bg_sq)
517 if (use_bbl_eos .and. (hweight >= 0.0))
then
518 thtot = thtot + hweight * t_vel(i,k)
519 shtot = shtot + hweight * s_vel(i,k)
523 if (.not.cs%linear_drag .and. (hwtot > 0.0))
then
524 ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
526 ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
529 if (use_bbl_eos)
then ;
if (hwtot > 0.0)
then
530 t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
532 t_eos(i) = 0.0 ; s_eos(i) = 0.0
536 do i=is,ie ; ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel ;
enddo
539 if (use_bbl_eos)
then
542 if (.not.do_i(i))
then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ;
endif
544 do k=1,nz ;
do i=is,ie
545 press(i) = press(i) + gv%H_to_Pa * h_vel(i,k)
548 is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
551 do i=is,ie ;
if (do_i(i))
then
554 ustarsq = rho0x400_g * ustar(i)**2
561 if (use_bbl_eos)
then
562 thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
564 if (h_at_vel(i,k) <= 0.0) cycle
566 oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
567 dr_ds(i)*(shtot - s_vel(i,k)*htot)
568 if (oldfn >= ustarsq)
exit
570 dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
571 dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
572 (h_at_vel(i,k) + htot)
574 if ((oldfn + dfn) <= ustarsq)
then
577 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
581 thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
583 if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0)
then
585 if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
586 dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) &
587 htot = htot + h_at_vel(i,1)
592 oldfn = rhtot - gv%Rlay(k)*htot
593 dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)
595 if (oldfn >= ustarsq)
then
597 elseif ((oldfn + dfn) <= ustarsq)
then
600 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
604 rhtot = rhtot + gv%Rlay(k)*dh
608 oldfn = rhtot - rml_vel(i,k)*htot
609 dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)
611 if (oldfn >= ustarsq)
then
613 elseif ((oldfn + dfn) <= ustarsq)
then
616 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
620 rhtot = rhtot + rml_vel(i,k)*dh
622 if (rhtot - rml_vel(i,1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
624 if (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
632 if (m==1)
then ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
633 else ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ;
endif
635 if (cs%cdrag * u_bg_sq <= 0.0)
then
638 usth = ustar(i)*gv%Z_to_H ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
639 if (htot*usth <= (cs%BBL_thick_min+h_neglect) * (0.5*usth + root))
then
640 bbl_thick = cs%BBL_thick_min
642 bbl_thick = (htot * usth) / (0.5*usth + root)
645 bbl_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
646 ((ustar(i)*ustar(i)) * (gv%Z_to_H**2)) ) )
648 if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
655 if ((bbl_thick > 0.5*cs%Hbbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%Hbbl
657 if (cs%Channel_drag)
then
663 tmp = g%mask2dCu(i,j+1) * d_u(i,j+1)
664 dp = 2.0 * d_vel * tmp / (d_vel + tmp)
665 tmp = g%mask2dCu(i,j-1) * d_u(i,j-1)
666 dm = 2.0 * d_vel * tmp / (d_vel + tmp)
669 tmp = g%mask2dCv(i+1,j) * d_v(i+1,j)
670 dp = 2.0 * d_vel * tmp / (d_vel + tmp)
671 tmp = g%mask2dCv(i-1,j) * d_v(i-1,j)
672 dm = 2.0 * d_vel * tmp / (d_vel + tmp)
674 if (dm > dp)
then ; tmp = dp ; dp = dm ; dm = tmp ;
endif
677 dp = gv%Z_to_H*dp ; dm = gv%Z_to_H*dm ; d_vel = gv%Z_to_H*d_vel
679 a_3 = (dp + dm - 2.0*d_vel) ; a = 3.0*a_3 ; a_12 = 0.25*a_3
683 if (abs(a) < 1e-2*(slope + cs%BBL_thick_min)) a = 0.0
691 vol_open = d_vel - dm ; vol_2_reg = vol_open
694 vol_open = 0.25*slope*tmp + c1_12*a
695 vol_2_reg = 0.5*tmp**2 * (a - c1_3*slope)
698 c24_a = 24.0/a ; iapb = 1.0/(a+slope)
699 apb_4a = (slope+a)/(4.0*a) ; a2x48_apb3 = (48.0*(a*a))*(iapb**3)
700 ax2_3apb = 2.0*c1_3*a*iapb
701 elseif (a == 0.0)
then
703 if (slope > 0) iapb = 1.0/slope
705 vol_open = d_vel - dm
706 if (slope >= -a)
then
707 iapb = 1.0e30 ;
if (slope+a /= 0.0) iapb = 1.0/(a+slope)
708 vol_direct = 0.0 ; l_direct = 0.0 ; c24_a = 0.0
710 c24_a = 24.0/a ; iapb = 1.0/(a+slope)
711 l_direct = 1.0 + slope/a
712 vol_direct = -c1_6*a*l_direct**3
714 ibma_2 = 2.0 / (slope - a)
717 l(nz+1) = 0.0 ; vol = 0.0 ; vol_err = 0.0 ; bbl_visc_frac = 0.0
722 vol = vol + h_vel(i,k)
723 h_vel_pos = h_vel(i,k) + h_neglect
725 if (vol >= vol_open)
then ; l(k) = 1.0
727 l(k) = sqrt(2.0*vol*iapb)
731 if (vol < vol_2_reg)
then
734 if (a2x48_apb3*vol < 1e-8)
then
736 l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2_3apb*l0)
738 l(k) = apb_4a * (1.0 - &
739 2.0 * cos(c1_3*acos(a2x48_apb3*vol - 1.0) - c2pi_3))
747 tmp_val_m1_to_p1 = 1.0 - c24_a*(vol_open - vol)
748 tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
749 l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
754 if (vol <= vol_direct)
then
756 l(k) = (-0.25*c24_a*vol)**c1_3
764 if (vol_below + vol_err <= vol_direct)
then
765 l0 = l_direct ; vol_0 = vol_direct
767 l0 = l(k+1) ; vol_0 = vol_below + vol_err
773 dv_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol_0)
777 if (.not.cs%answers_2018)
then
778 vol_tol = max(0.5*gv%Angstrom_H + gv%H_subroundoff, 1e-14*vol)
779 vol_quit = max(0.9*gv%Angstrom_H + gv%H_subroundoff, 1e-14*vol)
782 if ((.not.cs%answers_2018) .and. (dvol <= 0.0))
then
784 vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
785 elseif ( ((.not.cs%answers_2018) .and. &
786 (a*a*dvol**3 < vol_tol*dv_dl2**2 *(dv_dl2*vol_tol - 2.0*a*l0*dvol))) .or. &
787 (cs%answers_2018 .and. (a*a*dvol**3 < gv%Angstrom_H*dv_dl2**2 * &
788 (0.25*dv_dl2*gv%Angstrom_H - a*l0*dvol) )) )
then
791 l(k) = sqrt(l0*l0 + dvol / dv_dl2)
792 vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
794 if (dv_dl2*(1.0-l0*l0) < dvol + &
795 dv_dl2 * (vol_open - vol)*ibma_2)
then
796 l_max = sqrt(1.0 - (vol_open - vol)*ibma_2)
798 l_max = sqrt(l0*l0 + dvol / dv_dl2)
800 l_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l_max))
802 vol_err_min = 0.5*(l_min**2)*(slope + a_3*(3.0-4.0*l_min)) - vol
803 vol_err_max = 0.5*(l_max**2)*(slope + a_3*(3.0-4.0*l_max)) - vol
805 if (abs(vol_err_min) <= vol_quit)
then
806 l(k) = l_min ; vol_err = vol_err_min
808 l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
809 (vol_err_max - vol_err_min))
811 vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
812 if (abs(vol_err) <= vol_quit)
exit
815 l(k) = l(k) - vol_err / (l(k)* (slope + a - 2.0*a*l(k)))
826 if (l(k) > l(k+1))
then
827 if (vol_below < bbl_thick)
then
828 bbl_frac = (1.0-vol_below/bbl_thick)**2
829 bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
834 if (m==1)
then ; cell_width = g%dy_Cu(i,j)
835 else ; cell_width = g%dx_Cv(i,j) ;
endif
836 gam = 1.0 - l(k+1)/l(k)
837 rayleigh = us%m_to_Z * cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
838 (12.0*cs%c_Smag*h_vel_pos) / (12.0*cs%c_Smag*h_vel_pos + &
839 gv%m_to_H * cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
845 if (rayleigh > 0.0)
then
846 v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
847 visc%Ray_u(i,j,k) = rayleigh*us%T_to_s*sqrt(u(i,j,k)*u(i,j,k) + &
848 v_at_u*v_at_u + u_bg_sq)
849 else ; visc%Ray_u(i,j,k) = 0.0 ;
endif
851 if (rayleigh > 0.0)
then
852 u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
853 visc%Ray_v(i,j,k) = rayleigh*us%T_to_s*sqrt(v(i,j,k)*v(i,j,k) + &
854 u_at_v*u_at_v + u_bg_sq)
855 else ; visc%Ray_v(i,j,k) = 0.0 ;
endif
860 bbl_thick_z = bbl_thick * gv%H_to_Z
862 visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, &
863 cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
864 visc%bbl_thick_u(i,j) = bbl_thick_z
866 visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, &
867 cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
868 visc%bbl_thick_v(i,j) = bbl_thick_z
874 bbl_thick_z = bbl_thick * gv%H_to_Z
876 visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
877 visc%bbl_thick_u(i,j) = bbl_thick_z
879 visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
880 visc%bbl_thick_v(i,j) = bbl_thick_z
887 if (cs%id_bbl_thick_u > 0) &
888 call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
889 if (cs%id_kv_bbl_u > 0) &
890 call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
891 if (cs%id_bbl_thick_v > 0) &
892 call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
893 if (cs%id_kv_bbl_v > 0) &
894 call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
895 if (cs%id_Ray_u > 0) &
896 call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
897 if (cs%id_Ray_v > 0) &
898 call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)
901 if (
associated(visc%Ray_u) .and.
associated(visc%Ray_v)) &
902 call uvchksum(
"Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI, haloshift=0, scale=us%Z_to_m)
903 if (
associated(visc%kv_bbl_u) .and.
associated(visc%kv_bbl_v)) &
904 call uvchksum(
"kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
905 if (
associated(visc%bbl_thick_u) .and.
associated(visc%bbl_thick_v)) &
906 call uvchksum(
"bbl_thick_[uv]", visc%bbl_thick_u, &
907 visc%bbl_thick_v, g%HI, haloshift=0, scale=us%Z_to_m)
910 end subroutine set_viscous_bbl
913 function set_v_at_u(v, h, G, i, j, k, mask2dCv, OBC)
915 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
917 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
919 integer,
intent(in) :: i
920 integer,
intent(in) :: j
921 integer,
intent(in) :: k
922 real,
dimension(SZI_(G),SZJB_(G)),&
923 intent(in) :: mask2dcv
928 real :: hwt(0:1,-1:0)
930 integer :: i0, j0, i1, j1
932 do j0 = -1,0 ;
do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
933 hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
936 if (
associated(obc))
then ;
if (obc%number_of_segments > 0)
then
937 do j0 = -1,0 ;
do i0 = 0,1 ;
if ((obc%segnum_v(i+i0,j+j0) /= obc_none))
then
938 i1 = i+i0 ; j1 = j+j0
939 if (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_n)
then
940 hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
941 elseif (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_s)
then
942 hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
944 endif ;
enddo ;
enddo
947 hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
949 if (hwt_tot > 0.0) set_v_at_u = &
950 ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
951 (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt_tot
953 end function set_v_at_u
956 function set_u_at_v(u, h, G, i, j, k, mask2dCu, OBC)
958 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
960 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
962 integer,
intent(in) :: i
963 integer,
intent(in) :: j
964 integer,
intent(in) :: k
965 real,
dimension(SZIB_(G),SZJ_(G)), &
966 intent(in) :: mask2dcu
971 real :: hwt(-1:0,0:1)
973 integer :: i0, j0, i1, j1
975 do j0 = 0,1 ;
do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
976 hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
979 if (
associated(obc))
then ;
if (obc%number_of_segments > 0)
then
980 do j0 = 0,1 ;
do i0 = -1,0 ;
if ((obc%segnum_u(i+i0,j+j0) /= obc_none))
then
981 i1 = i+i0 ; j1 = j+j0
982 if (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_e)
then
983 hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
984 elseif (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_w)
then
985 hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
987 endif ;
enddo ;
enddo
990 hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
992 if (hwt_tot > 0.0) set_u_at_v = &
993 ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
994 (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt_tot
996 end function set_u_at_v
1003 subroutine set_viscous_ml(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetrize)
1007 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1009 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1011 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1019 real,
intent(in) :: dt
1022 logical,
optional,
intent(in) :: symmetrize
1026 real,
dimension(SZIB_(G)) :: &
1047 real,
dimension(SZIB_(G),SZJ_(G)) :: &
1050 real,
dimension(SZI_(G),SZJB_(G)) :: &
1053 real :: h_at_vel(szib_(g),szk_(g))
1056 integer :: k_massive(szib_(g))
1093 real :: cdrag_sqrt_z
1117 logical :: use_eos, do_any, do_any_shelf, do_i(szib_(g))
1118 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, k2, nkmb, nkml, n
1121 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1122 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1123 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
1125 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_vert_friction(visc_ML): "//&
1126 "Module must be initialized before it is used.")
1127 if (.not.(cs%dynamic_viscous_ML .or.
associated(forces%frac_shelf_u) .or. &
1128 associated(forces%frac_shelf_v)) )
return
1130 if (
present(symmetrize))
then ;
if (symmetrize)
then
1131 jsq = js-1 ; isq = is-1
1134 rho0x400_g = 400.0*(gv%Rho0/(us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
1135 u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
1136 cdrag_sqrt = sqrt(cs%cdrag)
1137 cdrag_sqrt_z = us%m_to_Z * sqrt(cs%cdrag)
1140 use_eos =
associated(tv%eqn_of_state)
1141 dt_rho0 = dt/gv%H_to_kg_m2
1142 h_neglect = gv%H_subroundoff
1143 h_tiny = 2.0*gv%Angstrom_H + h_neglect
1144 g_h_rho0 = (gv%g_Earth*gv%H_to_Z) / gv%Rho0
1146 if (
associated(forces%frac_shelf_u) .neqv.
associated(forces%frac_shelf_v)) &
1147 call mom_error(fatal,
"set_viscous_ML: one of forces%frac_shelf_u and "//&
1148 "forces%frac_shelf_v is associated, but the other is not.")
1150 if (
associated(forces%frac_shelf_u))
then
1153 call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1154 call safe_alloc_ptr(visc%tbl_thick_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1155 call safe_alloc_ptr(visc%tbl_thick_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1156 call safe_alloc_ptr(visc%kv_tbl_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1157 call safe_alloc_ptr(visc%kv_tbl_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1164 do j=js-1,je ;
do i=is-1,ie+1
1165 mask_v(i,j) = g%mask2dCv(i,j)
1168 do j=js-1,je+1 ;
do i=is-1,ie
1169 mask_u(i,j) = g%mask2dCu(i,j)
1172 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
1175 if (.not. obc%segment(n)%on_pe) cycle
1177 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
1178 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je))
then
1179 do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
1180 if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
1181 if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
1183 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= je))
then
1184 do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
1185 if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
1186 if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
1195 if (cs%dynamic_viscous_ML)
then
1199 if (g%mask2dCu(i,j) < 0.5)
then
1200 do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
1202 do_i(i) = .true. ; do_any = .true.
1204 thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1205 uhtot(i) = dt_rho0 * forces%taux(i,j)
1206 vhtot(i) = 0.25 * dt_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
1207 (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
1209 if (cs%omega_frac >= 1.0)
then ; absf = 2.0*cs%omega ;
else
1210 absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
1211 if (cs%omega_frac > 0.0) &
1212 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1214 u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i+1,j)))
1215 idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1219 if (do_any)
then ;
do k=1,nz
1222 if (use_eos .and. (k==nkml+1))
then
1225 press(i) = gv%H_to_Pa * htot(i)
1227 i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
1228 t_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i+1,j,k2)*tv%T(i+1,j,k2)) * i_2hlay
1229 s_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i+1,j,k2)*tv%S(i+1,j,k2)) * i_2hlay
1232 isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
1235 do i=isq,ieq ;
if (do_i(i))
then
1237 hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
1238 if (hlay > h_tiny)
then
1239 i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
1240 v_at_u = 0.5 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1241 h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i_2hlay
1242 uh2 = us%m_s_to_L_T**2*((uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2)
1245 t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i_2hlay
1246 s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i_2hlay
1247 ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1248 dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1250 ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1253 if (ghprime > 0.0)
then
1254 ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1255 if (ribulk * uh2 <= (htot(i)**2) * ghprime)
then
1256 visc%nkml_visc_u(i,j) = real(k_massive(i))
1258 elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime)
then
1259 visc%nkml_visc_u(i,j) = real(k-1) + &
1260 ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1267 if (do_i(i)) do_any = .true.
1270 if (.not.do_any)
exit
1273 do i=isq,ieq ;
if (do_i(i))
then
1274 htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
1275 uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
1276 vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1277 h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
1279 thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k))
1280 shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k))
1282 rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
1287 if (do_any)
then ;
do i=isq,ieq ;
if (do_i(i))
then
1288 visc%nkml_visc_u(i,j) = k_massive(i)
1289 endif ;
enddo ;
endif
1292 do_any_shelf = .false.
1293 if (
associated(forces%frac_shelf_u))
then
1295 if (forces%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0)
then
1297 visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
1299 do_i(i) = .true. ; do_any_shelf = .true.
1304 if (do_any_shelf)
then
1305 do k=1,nz ;
do i=isq,ieq ;
if (do_i(i))
then
1306 if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0)
then
1307 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
1308 (h(i,j,k) + h(i+1,j,k) + h_neglect)
1310 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
1313 h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1314 endif ;
enddo ;
enddo
1316 do i=isq,ieq ;
if (do_i(i))
then
1317 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1318 thtot(i) = 0.0 ; shtot(i) = 0.0
1319 if (use_eos .or. .not.cs%linear_drag)
then ;
do k=1,nz
1320 if (htot_vel>=cs%Htbl_shelf)
exit
1321 hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1322 if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1324 htot_vel = htot_vel + h_at_vel(i,k)
1325 hwtot = hwtot + hweight
1327 if (.not.cs%linear_drag)
then
1328 v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
1329 hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
1330 v_at_u**2 + u_bg_sq)
1333 thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1334 shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1338 if ((.not.cs%linear_drag) .and. (hwtot > 0.0))
then
1339 ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
1341 ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
1344 if (use_eos)
then ;
if (hwtot > 0.0)
then
1345 t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1347 t_eos(i) = 0.0 ; s_eos(i) = 0.0
1353 dr_dt, dr_ds, isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
1356 do i=isq,ieq ;
if (do_i(i))
then
1359 ustarsq = rho0x400_g * ustar(i)**2
1362 thtot(i) = 0.0 ; shtot(i) = 0.0
1364 if (h_at_vel(i,k) <= 0.0) cycle
1365 t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1366 s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1367 oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1368 if (oldfn >= ustarsq)
exit
1370 dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
1371 dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
1372 (h_at_vel(i,k)+htot(i))
1373 if ((oldfn + dfn) <= ustarsq)
then
1376 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1379 htot(i) = htot(i) + dh
1380 thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1382 if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0))
then
1383 t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
1384 s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
1385 if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1386 dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1387 htot(i) = htot(i) + h_at_vel(i,nz)
1392 rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1394 oldfn = rlay*htot(i) - rhtot(i)
1395 if (oldfn >= ustarsq)
exit
1397 dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1398 if ((oldfn + dfn) <= ustarsq)
then
1401 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1404 htot(i) = htot(i) + dh
1405 rhtot(i) = rhtot(i) + rlay*dh
1407 if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1408 htot(i) = htot(i) + h_at_vel(i,nz)
1415 ustar1 = ustar(i)*gv%Z_to_H
1416 h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1417 tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1418 ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1419 visc%tbl_thick_shelf_u(i,j) = tbl_thick_z
1420 visc%Kv_tbl_shelf_u(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1430 if (cs%dynamic_viscous_ML)
then
1434 if (g%mask2dCv(i,j) < 0.5)
then
1435 do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
1437 do_i(i) = .true. ; do_any = .true.
1439 thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1440 vhtot(i) = dt_rho0 * forces%tauy(i,j)
1441 uhtot(i) = 0.25 * dt_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
1442 (forces%taux(i-1,j) + forces%taux(i,j+1)))
1444 if (cs%omega_frac >= 1.0)
then ; absf = 2.0*cs%omega ;
else
1445 absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1446 if (cs%omega_frac > 0.0) &
1447 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1450 u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i,j+1)))
1451 idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1456 if (do_any)
then ;
do k=1,nz
1459 if (use_eos .and. (k==nkml+1))
then
1462 press(i) = gv%H_to_Pa * htot(i)
1464 i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
1465 t_eos(i) = (h(i,j,k2)*tv%T(i,j,k2) + h(i,j+1,k2)*tv%T(i,j+1,k2)) * i_2hlay
1466 s_eos(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i,j+1,k2)*tv%S(i,j+1,k2)) * i_2hlay
1469 is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
1472 do i=is,ie ;
if (do_i(i))
then
1474 hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
1475 if (hlay > h_tiny)
then
1476 i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
1477 u_at_v = 0.5 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1478 h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i_2hlay
1479 uh2 = us%m_s_to_L_T**2*((uhtot(i) - htot(i)*u_at_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2)
1482 t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i_2hlay
1483 s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i_2hlay
1484 ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1485 dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1487 ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1490 if (ghprime > 0.0)
then
1491 ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1492 if (ribulk * uh2 <= htot(i)**2 * ghprime)
then
1493 visc%nkml_visc_v(i,j) = real(k_massive(i))
1495 elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime)
then
1496 visc%nkml_visc_v(i,j) = real(k-1) + &
1497 ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1504 if (do_i(i)) do_any = .true.
1507 if (.not.do_any)
exit
1510 do i=is,ie ;
if (do_i(i))
then
1511 htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
1512 vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
1513 uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1514 h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
1516 thtot(i) = thtot(i) + 0.5 * (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k))
1517 shtot(i) = shtot(i) + 0.5 * (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k))
1519 rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
1524 if (do_any)
then ;
do i=is,ie ;
if (do_i(i))
then
1525 visc%nkml_visc_v(i,j) = k_massive(i)
1526 endif ;
enddo ;
endif
1529 do_any_shelf = .false.
1530 if (
associated(forces%frac_shelf_v))
then
1532 if (forces%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0)
then
1534 visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
1536 do_i(i) = .true. ; do_any_shelf = .true.
1541 if (do_any_shelf)
then
1542 do k=1,nz ;
do i=is,ie ;
if (do_i(i))
then
1543 if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0)
then
1544 h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
1545 (h(i,j,k) + h(i,j+1,k) + h_neglect)
1547 h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
1550 h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1551 endif ;
enddo ;
enddo
1553 do i=is,ie ;
if (do_i(i))
then
1554 htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1555 thtot(i) = 0.0 ; shtot(i) = 0.0
1556 if (use_eos .or. .not.cs%linear_drag)
then ;
do k=1,nz
1557 if (htot_vel>=cs%Htbl_shelf)
exit
1558 hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1559 if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1561 htot_vel = htot_vel + h_at_vel(i,k)
1562 hwtot = hwtot + hweight
1564 if (.not.cs%linear_drag)
then
1565 u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
1566 hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
1567 u_at_v**2 + u_bg_sq)
1570 thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1571 shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1575 if (.not.cs%linear_drag)
then ;
if (hwtot > 0.0)
then
1576 ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
1578 ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
1581 if (use_eos)
then ;
if (hwtot > 0.0)
then
1582 t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1584 t_eos(i) = 0.0 ; s_eos(i) = 0.0
1590 dr_dt, dr_ds, is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
1593 do i=is,ie ;
if (do_i(i))
then
1596 ustarsq = rho0x400_g * ustar(i)**2
1599 thtot(i) = 0.0 ; shtot(i) = 0.0
1601 if (h_at_vel(i,k) <= 0.0) cycle
1602 t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1603 s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1604 oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1605 if (oldfn >= ustarsq)
exit
1607 dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
1608 dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
1609 (h_at_vel(i,k)+htot(i))
1610 if ((oldfn + dfn) <= ustarsq)
then
1613 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1616 htot(i) = htot(i) + dh
1617 thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1619 if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0))
then
1620 t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
1621 s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
1622 if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1623 dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1624 htot(i) = htot(i) + h_at_vel(i,nz)
1629 rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1631 oldfn = rlay*htot(i) - rhtot(i)
1632 if (oldfn >= ustarsq)
exit
1634 dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1635 if ((oldfn + dfn) <= ustarsq)
then
1638 dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1641 htot(i) = htot(i) + dh
1642 rhtot = rhtot + rlay*dh
1644 if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1645 htot(i) = htot(i) + h_at_vel(i,nz)
1652 ustar1 = ustar(i)*gv%Z_to_H
1653 h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1654 tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1655 ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1656 visc%tbl_thick_shelf_v(i,j) = tbl_thick_z
1657 visc%Kv_tbl_shelf_v(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1665 if (
associated(visc%nkml_visc_u) .and.
associated(visc%nkml_visc_v)) &
1666 call uvchksum(
"nkml_visc_[uv]", visc%nkml_visc_u, &
1667 visc%nkml_visc_v, g%HI,haloshift=0)
1669 if (cs%id_nkml_visc_u > 0) &
1670 call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
1671 if (cs%id_nkml_visc_v > 0) &
1672 call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)
1674 end subroutine set_viscous_ml
1677 subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
1687 logical :: use_kappa_shear, ks_at_vertex
1688 logical :: adiabatic, usekpp, useepbl
1689 logical :: use_cvmix_shear, mle_use_pbl_mld, use_cvmix_conv
1690 integer :: isd, ied, jsd, jed, nz
1692 character(len=40) :: mdl =
"MOM_set_visc"
1693 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1695 call get_param(param_file, mdl,
"ADIABATIC", adiabatic, default=.false., &
1698 use_kappa_shear = .false. ; ks_at_vertex = .false. ; use_cvmix_shear = .false.
1699 usekpp = .false. ; useepbl = .false. ; use_cvmix_conv = .false.
1701 if (.not.adiabatic)
then
1702 use_kappa_shear = kappa_shear_is_used(param_file)
1703 ks_at_vertex = kappa_shear_at_vertex(param_file)
1704 use_cvmix_shear = cvmix_shear_is_used(param_file)
1705 use_cvmix_conv = cvmix_conv_is_used(param_file)
1706 call get_param(param_file, mdl,
"USE_KPP", usekpp, &
1707 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1708 "to calculate diffusivities and non-local transport in the OBL.", &
1709 default=.false., do_not_log=.true.)
1710 call get_param(param_file, mdl,
"ENERGETICS_SFC_PBL", useepbl, &
1711 "If true, use an implied energetics planetary boundary "//&
1712 "layer scheme to determine the diffusivity and viscosity "//&
1713 "in the surface boundary layer.", default=.false., do_not_log=.true.)
1716 if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv)
then
1719 "Shear-driven turbulent diffusivity at interfaces",
"m2 s-1", z_grid=
'i')
1721 if (usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv .or. &
1722 (use_kappa_shear .and. .not.ks_at_vertex ))
then
1725 "Shear-driven turbulent viscosity at interfaces",
"m2 s-1", z_grid=
'i')
1727 if (use_kappa_shear .and. ks_at_vertex)
then
1728 call safe_alloc_ptr(visc%TKE_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1729 call safe_alloc_ptr(visc%Kv_shear_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1731 "Shear-driven turbulent viscosity at vertex interfaces",
"m2 s-1", &
1732 hor_grid=
"Bu", z_grid=
'i')
1733 elseif (use_kappa_shear)
then
1740 "Vertical turbulent viscosity at interfaces due to slow processes", &
1741 "m2 s-1", z_grid=
'i')
1744 call get_param(param_file, mdl,
"MLE_USE_PBL_MLD", mle_use_pbl_mld, &
1745 default=.false., do_not_log=.true.)
1747 call get_param(param_file, mdl,
"HFREEZE", hfreeze, &
1748 default=-1.0, do_not_log=.true.)
1750 if (mle_use_pbl_mld)
then
1753 "Instantaneous active mixing layer depth",
"m")
1756 if (hfreeze >= 0.0 .and. .not.mle_use_pbl_mld)
then
1761 end subroutine set_visc_register_restarts
1764 subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
1765 type(time_type),
target,
intent(in) :: time
1771 type(
diag_ctrl),
target,
intent(inout) :: diag
1781 real :: csmag_chan_dflt, smag_const1, tke_decay_dflt, bulk_ri_ml_dflt
1782 real :: kv_background
1783 real :: omega_frac_dflt
1788 real :: z2_t_rescale
1790 integer :: i, j, k, is, ie, js, je, n
1791 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
1792 logical :: default_2018_answers
1793 logical :: use_kappa_shear, adiabatic, use_omega
1794 logical :: use_cvmix_ddiff, differential_diffusion, use_kpp
1797 # include "version_variable.h"
1798 character(len=40) :: mdl =
"MOM_set_visc"
1800 if (
associated(cs))
then
1801 call mom_error(warning,
"set_visc_init called with an associated "// &
1802 "control structure.")
1809 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1810 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1811 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1817 cs%RiNo_mix = .false. ; use_cvmix_ddiff = .false.
1818 differential_diffusion = .false.
1819 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1820 "This sets the default value for the various _2018_ANSWERS parameters.", &
1822 call get_param(param_file, mdl,
"SET_VISC_2018_ANSWERS", cs%answers_2018, &
1823 "If true, use the order of arithmetic and expressions that recover the "//&
1824 "answers from the end of 2018. Otherwise, use updated and more robust "//&
1825 "forms of the same expressions.", default=default_2018_answers)
1826 call get_param(param_file, mdl,
"BOTTOMDRAGLAW", cs%bottomdraglaw, &
1827 "If true, the bottom stress is calculated with a drag "//&
1828 "law of the form c_drag*|u|*u. The velocity magnitude "//&
1829 "may be an assumed value or it may be based on the "//&
1830 "actual velocity in the bottommost HBBL, depending on "//&
1831 "LINEAR_DRAG.", default=.true.)
1832 call get_param(param_file, mdl,
"CHANNEL_DRAG", cs%Channel_drag, &
1833 "If true, the bottom drag is exerted directly on each "//&
1834 "layer proportional to the fraction of the bottom it "//&
1835 "overlies.", default=.false.)
1836 call get_param(param_file, mdl,
"LINEAR_DRAG", cs%linear_drag, &
1837 "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
1838 "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
1839 call get_param(param_file, mdl,
"ADIABATIC", adiabatic, default=.false., &
1842 call log_param(param_file, mdl,
"ADIABATIC",adiabatic, &
1843 "There are no diapycnal mass fluxes if ADIABATIC is "//&
1844 "true. This assumes that KD = KDML = 0.0 and that "//&
1845 "there is no buoyancy forcing, but makes the model "//&
1846 "faster by eliminating subroutine calls.", default=.false.)
1849 if (.not.adiabatic)
then
1850 cs%RiNo_mix = kappa_shear_is_used(param_file)
1851 call get_param(param_file, mdl,
"DOUBLE_DIFFUSION", differential_diffusion, &
1852 "If true, increase diffusivites for temperature or salt "//&
1853 "based on double-diffusive parameterization from MOM4/KPP.", &
1855 use_cvmix_ddiff = cvmix_ddiff_is_used(param_file)
1858 call get_param(param_file, mdl,
"PRANDTL_TURB", visc%Prandtl_turb, &
1859 "The turbulent Prandtl number applied to shear "//&
1860 "instability.", units=
"nondim", default=1.0)
1861 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false.)
1863 call get_param(param_file, mdl,
"DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
1864 "If true, use a bulk Richardson number criterion to "//&
1865 "determine the mixed layer thickness for viscosity.", &
1867 if (cs%dynamic_viscous_ML)
then
1868 call get_param(param_file, mdl,
"BULK_RI_ML", bulk_ri_ml_dflt, default=0.0)
1869 call get_param(param_file, mdl,
"BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
1870 "The efficiency with which mean kinetic energy released "//&
1871 "by mechanically forced entrainment of the mixed layer "//&
1872 "is converted to turbulent kinetic energy. By default, "//&
1873 "BULK_RI_ML_VISC = BULK_RI_ML or 0.", units=
"nondim", &
1874 default=bulk_ri_ml_dflt)
1875 call get_param(param_file, mdl,
"TKE_DECAY", tke_decay_dflt, default=0.0)
1876 call get_param(param_file, mdl,
"TKE_DECAY_VISC", cs%TKE_decay, &
1877 "TKE_DECAY_VISC relates the vertical rate of decay of "//&
1878 "the TKE available for mechanical entrainment to the "//&
1879 "natural Ekman depth for use in calculating the dynamic "//&
1880 "mixed layer viscosity. By default, "//&
1881 "TKE_DECAY_VISC = TKE_DECAY or 0.", units=
"nondim", &
1882 default=tke_decay_dflt)
1883 call get_param(param_file, mdl,
"ML_USE_OMEGA", use_omega, &
1884 "If true, use the absolute rotation rate instead of the "//&
1885 "vertical component of rotation when setting the decay "//&
1886 "scale for turbulence.", default=.false., do_not_log=.true.)
1887 omega_frac_dflt = 0.0
1889 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1890 omega_frac_dflt = 1.0
1892 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%omega_frac, &
1893 "When setting the decay scale for turbulence, use this "//&
1894 "fraction of the absolute rotation rate blended with the "//&
1895 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1896 units=
"nondim", default=omega_frac_dflt)
1897 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
1898 "The rotation rate of the earth.", units=
"s-1", &
1899 default=7.2921e-5, scale=us%T_to_s)
1901 cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
1903 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
1904 "The rotation rate of the earth.", units=
"s-1", &
1905 default=7.2921e-5, scale=us%T_to_s)
1908 call get_param(param_file, mdl,
"HBBL", cs%Hbbl, &
1909 "The thickness of a bottom boundary layer with a "//&
1910 "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
1911 "the thickness over which near-bottom velocities are "//&
1912 "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
1913 "but LINEAR_DRAG is not.", units=
"m", fail_if_missing=.true.)
1914 if (cs%bottomdraglaw)
then
1915 call get_param(param_file, mdl,
"CDRAG", cs%cdrag, &
1916 "CDRAG is the drag coefficient relating the magnitude of "//&
1917 "the velocity field to the bottom stress. CDRAG is only "//&
1918 "used if BOTTOMDRAGLAW is defined.", units=
"nondim", &
1920 call get_param(param_file, mdl,
"DRAG_BG_VEL", cs%drag_bg_vel, &
1921 "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1922 "LINEAR_DRAG) or an unresolved velocity that is "//&
1923 "combined with the resolved velocity to estimate the "//&
1924 "velocity magnitude. DRAG_BG_VEL is only used when "//&
1925 "BOTTOMDRAGLAW is defined.", units=
"m s-1", default=0.0)
1926 call get_param(param_file, mdl,
"BBL_USE_EOS", cs%BBL_use_EOS, &
1927 "If true, use the equation of state in determining the "//&
1928 "properties of the bottom boundary layer. Otherwise use "//&
1929 "the layer target potential densities.", default=.false.)
1931 call get_param(param_file, mdl,
"BBL_THICK_MIN", cs%BBL_thick_min, &
1932 "The minimum bottom boundary layer thickness that can be "//&
1933 "used with BOTTOMDRAGLAW. This might be "//&
1934 "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1935 "near-bottom viscosity.", units=
"m", default=0.0)
1936 call get_param(param_file, mdl,
"HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
1937 "The minimum top boundary layer thickness that can be "//&
1938 "used with BOTTOMDRAGLAW. This might be "//&
1939 "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1940 "near-top viscosity.", units=
"m", default=cs%BBL_thick_min, scale=gv%m_to_H)
1941 call get_param(param_file, mdl,
"HTBL_SHELF", cs%Htbl_shelf, &
1942 "The thickness over which near-surface velocities are "//&
1943 "averaged for the drag law under an ice shelf. By "//&
1944 "default this is the same as HBBL", units=
"m", default=cs%Hbbl, scale=gv%m_to_H)
1946 cs%Hbbl = cs%Hbbl * gv%m_to_H
1947 cs%BBL_thick_min = cs%BBL_thick_min * gv%m_to_H
1949 call get_param(param_file, mdl,
"KV", kv_background, &
1950 "The background kinematic viscosity in the interior. "//&
1951 "The molecular value, ~1e-6 m2 s-1, may be used.", &
1952 units=
"m2 s-1", fail_if_missing=.true.)
1954 call get_param(param_file, mdl,
"ADD_KV_SLOW", visc%add_Kv_slow, &
1955 "If true, the background vertical viscosity in the interior "//&
1956 "(i.e., tidal + background + shear + convection) is added "//&
1957 "when computing the coupling coefficient. The purpose of this "//&
1958 "flag is to be able to recover previous answers and it will likely "//&
1959 "be removed in the future since this option should always be true.", &
1962 call get_param(param_file, mdl,
"USE_KPP", use_kpp, &
1963 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1964 "to calculate diffusivities and non-local transport in the OBL.", &
1965 do_not_log=.true., default=.false.)
1967 if (use_kpp .and. visc%add_Kv_slow)
call mom_error(fatal,
"set_visc_init: "//&
1968 "When USE_KPP=True, ADD_KV_SLOW must be false. Otherwise vertical "//&
1969 "viscosity due to slow processes will be double counted. Please set "//&
1970 "ADD_KV_SLOW=False.")
1972 call get_param(param_file, mdl,
"KV_BBL_MIN", cs%KV_BBL_min, &
1973 "The minimum viscosities in the bottom boundary layer.", &
1974 units=
"m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1975 call get_param(param_file, mdl,
"KV_TBL_MIN", cs%KV_TBL_min, &
1976 "The minimum viscosities in the top boundary layer.", &
1977 units=
"m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1979 if (cs%Channel_drag)
then
1980 call get_param(param_file, mdl,
"SMAG_LAP_CONST", smag_const1, default=-1.0)
1982 csmag_chan_dflt = 0.15
1983 if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1
1985 call get_param(param_file, mdl,
"SMAG_CONST_CHANNEL", cs%c_Smag, &
1986 "The nondimensional Laplacian Smagorinsky constant used "//&
1987 "in calculating the channel drag if it is enabled. The "//&
1988 "default is to use the same value as SMAG_LAP_CONST if "//&
1989 "it is defined, or 0.15 if it is not. The value used is "//&
1990 "also 0.15 if the specified value is negative.", &
1991 units=
"nondim", default=csmag_chan_dflt)
1992 if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
1995 if (cs%RiNo_mix .and. kappa_shear_at_vertex(param_file))
then
1997 call pass_var(visc%Kv_shear_Bu, g%Domain, position=corner, complete=.true.)
2000 if (cs%bottomdraglaw)
then
2001 allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed)) ; visc%bbl_thick_u = 0.0
2002 allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed)) ; visc%kv_bbl_u = 0.0
2003 allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb)) ; visc%bbl_thick_v = 0.0
2004 allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb)) ; visc%kv_bbl_v = 0.0
2005 allocate(visc%ustar_bbl(isd:ied,jsd:jed)) ; visc%ustar_bbl = 0.0
2006 allocate(visc%TKE_bbl(isd:ied,jsd:jed)) ; visc%TKE_bbl = 0.0
2008 cs%id_bbl_thick_u = register_diag_field(
'ocean_model',
'bbl_thick_u', &
2009 diag%axesCu1, time,
'BBL thickness at u points',
'm', conversion=us%Z_to_m)
2010 cs%id_kv_bbl_u = register_diag_field(
'ocean_model',
'kv_bbl_u', diag%axesCu1, &
2011 time,
'BBL viscosity at u points',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
2012 cs%id_bbl_thick_v = register_diag_field(
'ocean_model',
'bbl_thick_v', &
2013 diag%axesCv1, time,
'BBL thickness at v points',
'm', conversion=us%Z_to_m)
2014 cs%id_kv_bbl_v = register_diag_field(
'ocean_model',
'kv_bbl_v', diag%axesCv1, &
2015 time,
'BBL viscosity at v points',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
2017 if (cs%Channel_drag)
then
2018 allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray_u = 0.0
2019 allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray_v = 0.0
2020 cs%id_Ray_u = register_diag_field(
'ocean_model',
'Rayleigh_u', diag%axesCuL, &
2021 time,
'Rayleigh drag velocity at u points',
'm s-1', conversion=us%Z_to_m*us%s_to_T)
2022 cs%id_Ray_v = register_diag_field(
'ocean_model',
'Rayleigh_v', diag%axesCvL, &
2023 time,
'Rayleigh drag velocity at v points',
'm s-1', conversion=us%Z_to_m*us%s_to_T)
2026 if (use_cvmix_ddiff .or. differential_diffusion)
then
2027 allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0
2028 allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0
2031 if (cs%dynamic_viscous_ML)
then
2032 allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed)) ; visc%nkml_visc_u = 0.0
2033 allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb)) ; visc%nkml_visc_v = 0.0
2034 cs%id_nkml_visc_u = register_diag_field(
'ocean_model',
'nkml_visc_u', &
2035 diag%axesCu1, time,
'Number of layers in viscous mixed layer at u points',
'm')
2036 cs%id_nkml_visc_v = register_diag_field(
'ocean_model',
'nkml_visc_v', &
2037 diag%axesCv1, time,
'Number of layers in viscous mixed layer at v points',
'm')
2040 call register_restart_field_as_obsolete(
'Kd_turb',
'Kd_shear', restart_cs)
2041 call register_restart_field_as_obsolete(
'Kv_turb',
'Kv_shear', restart_cs)
2046 if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) &
2047 z_rescale = us%m_to_Z / us%m_to_Z_restart
2049 if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
2050 i_t_rescale = us%s_to_T_restart / us%s_to_T
2051 z2_t_rescale = z_rescale**2*i_t_rescale
2053 if (z2_t_rescale /= 1.0)
then
2054 if (
associated(visc%Kd_shear))
then ;
if (
query_initialized(visc%Kd_shear,
"Kd_shear", restart_cs))
then
2055 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2056 visc%Kd_shear(i,j,k) = z2_t_rescale * visc%Kd_shear(i,j,k)
2057 enddo ;
enddo ;
enddo
2060 if (
associated(visc%Kv_shear))
then ;
if (
query_initialized(visc%Kv_shear,
"Kv_shear", restart_cs))
then
2061 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2062 visc%Kv_shear(i,j,k) = z2_t_rescale * visc%Kv_shear(i,j,k)
2063 enddo ;
enddo ;
enddo
2066 if (
associated(visc%Kv_shear_Bu))
then ;
if (
query_initialized(visc%Kv_shear_Bu,
"Kv_shear_Bu", restart_cs))
then
2067 do k=1,nz+1 ;
do j=js-1,je ;
do i=is-1,ie
2068 visc%Kv_shear_Bu(i,j,k) = z2_t_rescale * visc%Kv_shear_Bu(i,j,k)
2069 enddo ;
enddo ;
enddo
2072 if (
associated(visc%Kv_slow))
then ;
if (
query_initialized(visc%Kv_slow,
"Kv_slow", restart_cs))
then
2073 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2074 visc%Kv_slow(i,j,k) = z2_t_rescale * visc%Kv_slow(i,j,k)
2075 enddo ;
enddo ;
enddo
2079 end subroutine set_visc_init
2082 subroutine set_visc_end(visc, CS)
2087 if (cs%bottomdraglaw)
then
2088 deallocate(visc%bbl_thick_u) ;
deallocate(visc%bbl_thick_v)
2089 deallocate(visc%kv_bbl_u) ;
deallocate(visc%kv_bbl_v)
2091 if (cs%Channel_drag)
then
2092 deallocate(visc%Ray_u) ;
deallocate(visc%Ray_v)
2094 if (cs%dynamic_viscous_ML)
then
2095 deallocate(visc%nkml_visc_u) ;
deallocate(visc%nkml_visc_v)
2097 if (
associated(visc%Kd_shear))
deallocate(visc%Kd_shear)
2098 if (
associated(visc%Kv_slow))
deallocate(visc%Kv_slow)
2099 if (
associated(visc%TKE_turb))
deallocate(visc%TKE_turb)
2100 if (
associated(visc%Kv_shear))
deallocate(visc%Kv_shear)
2101 if (
associated(visc%Kv_shear_Bu))
deallocate(visc%Kv_shear_Bu)
2102 if (
associated(visc%ustar_bbl))
deallocate(visc%ustar_bbl)
2103 if (
associated(visc%TKE_bbl))
deallocate(visc%TKE_bbl)
2104 if (
associated(visc%taux_shelf))
deallocate(visc%taux_shelf)
2105 if (
associated(visc%tauy_shelf))
deallocate(visc%tauy_shelf)
2106 if (
associated(visc%tbl_thick_shelf_u))
deallocate(visc%tbl_thick_shelf_u)
2107 if (
associated(visc%tbl_thick_shelf_v))
deallocate(visc%tbl_thick_shelf_v)
2108 if (
associated(visc%kv_tbl_shelf_u))
deallocate(visc%kv_tbl_shelf_u)
2109 if (
associated(visc%kv_tbl_shelf_v))
deallocate(visc%kv_tbl_shelf_v)
2112 end subroutine set_visc_end