Calculates the acceleration due to the horizontal viscosity.
A combination of biharmonic and Laplacian forms can be used. The coefficient may either be a constant or a shear-dependent form. The biharmonic is determined by twice taking the divergence of an appropriately defined stress tensor. The Laplacian is determined by doing so once.
To work, the following fields must be set outside of the usual is:ie range before this subroutine is called: u[is-2:ie+2,js-2:je+2] v[is-2:ie+2,js-2:je+2] h[is-1:ie+1,js-1:je+1]
200 type(ocean_grid_type),
intent(in) :: G
201 type(verticalGrid_type),
intent(in) :: GV
202 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
204 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
206 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
208 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
211 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
214 type(MEKE_type),
pointer :: MEKE
216 type(VarMix_CS),
pointer :: VarMix
218 type(unit_scale_type),
intent(in) :: US
219 type(hor_visc_CS),
pointer :: CS
221 type(ocean_OBC_type),
optional,
pointer :: OBC
222 type(barotropic_CS),
optional,
pointer :: BT
226 real,
dimension(SZIB_(G),SZJ_(G)) :: &
232 real,
dimension(SZI_(G),SZJB_(G)) :: &
238 real,
dimension(SZI_(G),SZJ_(G)) :: &
252 grad_vort_mag_h_2d, &
261 real,
dimension(SZIB_(G),SZJB_(G)) :: &
275 grad_vort_mag_q_2d, &
282 real,
dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
293 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
298 target_diss_rate_GME, &
306 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
319 real :: vert_vort_mag
331 real :: visc_bound_rem
340 real :: GME_coeff_limiter
345 real :: backscat_subround
348 logical :: rescale_Kh, legacy_bound
349 logical :: find_FrictWork
350 logical :: apply_OBC = .false.
351 logical :: use_MEKE_Ku
352 logical :: use_MEKE_Au
353 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
354 integer :: i, j, k, n
355 real :: inv_PI3, inv_PI2, inv_PI5
356 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
357 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
359 h_neglect = gv%H_subroundoff
360 h_neglect3 = h_neglect**3
361 inv_pi3 = 1.0/((4.0*atan(1.0))**3)
362 inv_pi2 = 1.0/((4.0*atan(1.0))**2)
363 inv_pi5 = inv_pi3 * inv_pi2
368 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then
369 apply_obc = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally
371 endif ;
endif ;
endif
373 if (.not.
associated(cs))
call mom_error(fatal, &
374 "MOM_hor_visc: Module must be initialized before it is used.")
375 if (.not.(cs%Laplacian .or. cs%biharmonic))
return
377 find_frictwork = (cs%id_FrictWork > 0)
378 if (cs%id_FrictWorkIntz > 0) find_frictwork = .true.
379 if (
associated(meke))
then
380 if (
associated(meke%mom_src)) find_frictwork = .true.
381 backscat_subround = 0.0
382 if (find_frictwork .and.
associated(meke%mom_src) .and. (meke%backscatter_Ro_c > 0.0) .and. &
383 (meke%backscatter_Ro_Pow /= 0.0)) &
384 backscat_subround = (1.0e-16/meke%backscatter_Ro_c)**(1.0/meke%backscatter_Ro_Pow)
388 if (
associated(varmix))
then
389 rescale_kh = varmix%Resoln_scaled_Kh
390 if (rescale_kh .and. &
391 (.not.
associated(varmix%Res_fn_h) .or. .not.
associated(varmix%Res_fn_q))) &
392 call mom_error(fatal,
"MOM_hor_visc: VarMix%Res_fn_h and " //&
393 "VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
395 legacy_bound = (cs%Smagorinsky_Kh .or. cs%Leith_Kh) .and. &
396 (cs%bound_Kh .and. .not.cs%better_bound_Kh)
399 use_meke_ku =
associated(meke%Ku)
400 use_meke_au =
associated(meke%Au)
402 do j=js,je ;
do i=is,ie
403 boundary_mask(i,j) = (g%mask2dCu(i,j) * g%mask2dCv(i,j) * g%mask2dCu(i-1,j) * g%mask2dCv(i,j-1))
408 h0_gme = 1000.0*us%m_to_Z
410 gme_coeff_limiter = 1e7*us%T_to_s
413 gme_coeff_h(:,:,:) = 0.0
414 gme_coeff_q(:,:,:) = 0.0
415 str_xx_gme(:,:) = 0.0
416 str_xy_gme(:,:) = 0.0
421 call barotropic_get_tav(bt, ubtav, vbtav, g, us)
422 call pass_vector(ubtav, vbtav, g%Domain)
425 do j=js,je ;
do i=is,ie
426 dudx_bt(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * ubtav(i,j) - &
427 g%IdyCu(i-1,j) * ubtav(i-1,j))
428 dvdy_bt(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * vbtav(i,j) - &
429 g%IdxCv(i,j-1) * vbtav(i,j-1))
433 call pass_var(dudx_bt, g%Domain, complete=.true.)
434 call pass_var(dvdy_bt, g%Domain, complete=.true.)
438 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
439 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
443 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
444 dvdx_bt(i,j) = cs%DY_dxBu(i,j)*(vbtav(i+1,j)*g%IdyCv(i+1,j) &
445 - vbtav(i,j)*g%IdyCv(i,j))
446 dudy_bt(i,j) = cs%DX_dyBu(i,j)*(ubtav(i,j+1)*g%IdxCu(i,j+1) &
447 - ubtav(i,j)*g%IdxCu(i,j))
451 call pass_var(dvdx_bt, g%Domain, position=corner, complete=.true.)
452 call pass_var(dudy_bt, g%Domain, position=corner, complete=.true.)
457 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
458 sh_xy_bt(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
463 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
464 sh_xy_bt(i,j) = g%mask2dBu(i,j) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
473 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
474 grad_vel_mag_bt_h(i,j) = boundary_mask(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
475 (0.25*(dvdx_bt(i,j)+dvdx_bt(i-1,j)+dvdx_bt(i,j-1)+dvdx_bt(i-1,j-1)) )**2 + &
476 (0.25*(dudy_bt(i,j)+dudy_bt(i-1,j)+dudy_bt(i,j-1)+dudy_bt(i-1,j-1)) )**2)
480 if (
associated(meke))
then ;
if (
associated(meke%mom_src))
then
482 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
483 max_diss_rate_bt(i,j) = 2.0*meke%MEKE(i,j) * grad_vel_mag_bt_h(i,j)
491 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
492 grad_vel_mag_bt_q(i,j) = boundary_mask(i,j) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
493 (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
494 (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2)
520 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
521 dudx(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * u(i,j,k) - &
522 g%IdyCu(i-1,j) * u(i-1,j,k))
523 dvdy(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * v(i,j,k) - &
524 g%IdxCv(i,j-1) * v(i,j-1,k))
525 sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
529 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
530 dvdx(i,j) = cs%DY_dxBu(i,j)*(v(i+1,j,k)*g%IdyCv(i+1,j) - v(i,j,k)*g%IdyCv(i,j))
531 dudy(i,j) = cs%DX_dyBu(i,j)*(u(i,j+1,k)*g%IdxCu(i,j+1) - u(i,j,k)*g%IdxCu(i,j))
539 if (cs%use_land_mask)
then
540 do j=js-2,je+2 ;
do i=isq-1,ieq+1
541 h_u(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k))
543 do j=jsq-1,jeq+1 ;
do i=is-2,ie+2
544 h_v(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k))
547 do j=js-2,je+2 ;
do i=isq-1,ieq+1
548 h_u(i,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
550 do j=jsq-1,jeq+1 ;
do i=is-2,ie+2
551 h_v(i,j) = 0.5 * (h(i,j,k) + h(i,j+1,k))
557 if (apply_obc)
then ;
do n=1,obc%number_of_segments
558 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
559 if (obc%zero_strain .or. obc%freeslip_strain .or. obc%computed_strain)
then
560 if (obc%segment(n)%is_N_or_S .and. (j >= js-2) .and. (j <= jeq+1))
then
561 do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
562 if (obc%zero_strain)
then
563 dvdx(i,j) = 0. ; dudy(i,j) = 0.
564 elseif (obc%freeslip_strain)
then
566 elseif (obc%computed_strain)
then
567 if (obc%segment(n)%direction == obc_direction_n)
then
568 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)*(obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%IdxCu(i,j)
570 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)*(u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdxCu(i,j+1)
572 elseif (obc%specified_strain)
then
573 if (obc%segment(n)%direction == obc_direction_n)
then
574 dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j)*g%dxBu(i,j)
576 dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j+1)*g%dxBu(i,j)
580 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-2) .and. (i <= ieq+1))
then
581 do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
582 if (obc%zero_strain)
then
583 dvdx(i,j) = 0. ; dudy(i,j) = 0.
584 elseif (obc%freeslip_strain)
then
586 elseif (obc%computed_strain)
then
587 if (obc%segment(n)%direction == obc_direction_e)
then
588 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)*(obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%IdyCv(i,j)
590 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)*(v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdyCv(i+1,j)
592 elseif (obc%specified_strain)
then
593 if (obc%segment(n)%direction == obc_direction_e)
then
594 dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i,j)*g%dxBu(i,j)
596 dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i+1,j)*g%dxBu(i,j)
602 if (obc%segment(n)%direction == obc_direction_n)
then
607 if ((j >= jsq-1) .and. (j <= jeq+1))
then
608 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
612 elseif (obc%segment(n)%direction == obc_direction_s)
then
613 if ((j >= jsq-1) .and. (j <= jeq+1))
then
614 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
615 h_v(i,j) = h(i,j+1,k)
618 elseif (obc%segment(n)%direction == obc_direction_e)
then
619 if ((i >= isq-1) .and. (i <= ieq+1))
then
620 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
624 elseif (obc%segment(n)%direction == obc_direction_w)
then
625 if ((i >= isq-1) .and. (i <= ieq+1))
then
626 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
627 h_u(i,j) = h(i+1,j,k)
633 if (apply_obc)
then ;
do n=1,obc%number_of_segments
634 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
635 if (obc%segment(n)%direction == obc_direction_n)
then
636 if ((j >= js-2) .and. (j <= je))
then
637 do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
638 h_u(i,j+1) = h_u(i,j)
641 elseif (obc%segment(n)%direction == obc_direction_s)
then
642 if ((j >= js-1) .and. (j <= je+1))
then
643 do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+1,obc%segment(n)%HI%ied)
644 h_u(i,j) = h_u(i,j+1)
647 elseif (obc%segment(n)%direction == obc_direction_e)
then
648 if ((i >= is-2) .and. (i <= ie))
then
649 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
650 h_v(i+1,j) = h_v(i,j)
653 elseif (obc%segment(n)%direction == obc_direction_w)
then
654 if ((i >= is-1) .and. (i <= ie+1))
then
655 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
656 h_v(i,j) = h_v(i+1,j)
665 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
666 sh_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) + dudy(i,j) )
669 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
670 sh_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) + dudy(i,j) )
675 if (cs%biharmonic)
then
676 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
677 u0(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*sh_xx(i+1,j) - cs%DY2h(i,j)*sh_xx(i,j)) + &
678 cs%IDX2dyCu(i,j)*(cs%DX2q(i,j)*sh_xy(i,j) - cs%DX2q(i,j-1)*sh_xy(i,j-1))
680 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
681 v0(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j)*sh_xy(i,j) - cs%DY2q(i-1,j)*sh_xy(i-1,j)) - &
682 cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*sh_xx(i,j+1) - cs%DX2h(i,j)*sh_xx(i,j))
684 if (apply_obc) then;
if (obc%zero_biharmonic)
then
685 do n=1,obc%number_of_segments
686 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
687 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1))
then
688 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
691 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1))
then
692 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
700 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
706 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
707 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)
708 dvdx(i,j) = dy_dxbu * (v(i+1,j,k) * g%IdyCv(i+1,j) - v(i,j,k) * g%IdyCv(i,j))
709 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)
710 dudy(i,j) = dx_dybu * (u(i,j+1,k) * g%IdxCu(i,j+1) - u(i,j,k) * g%IdxCu(i,j))
715 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
716 vort_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) - dudy(i,j) )
717 dudy(i,j) = (2.0-g%mask2dBu(i,j)) * dudy(i,j)
718 dvdx(i,j) = (2.0-g%mask2dBu(i,j)) * dvdx(i,j)
721 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
722 vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )
723 dudy(i,j) = g%mask2dBu(i,j) * dudy(i,j)
724 dvdx(i,j) = g%mask2dBu(i,j) * dvdx(i,j)
728 call pass_var(vort_xy, g%Domain, position=corner, complete=.true.)
731 do j=js-2,jeq+1 ;
do i=is-1,ieq+1
732 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)
733 vort_xy_dx(i,j) = dy_dxbu * (vort_xy(i,j) * g%IdyCu(i,j) - vort_xy(i-1,j) * g%IdyCu(i-1,j))
736 do j=js-1,jeq+1 ;
do i=is-2,ieq+1
737 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)
738 vort_xy_dy(i,j) = dx_dybu * (vort_xy(i,j) * g%IdxCv(i,j) - vort_xy(i,j-1) * g%IdxCv(i,j-1))
741 call pass_vector(vort_xy_dy, vort_xy_dx, g%Domain)
743 if (cs%modified_Leith)
then
745 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
746 div_xx(i,j) = 0.5*((g%dyCu(i,j) * u(i,j,k) * (h(i+1,j,k)+h(i,j,k)) - &
747 g%dyCu(i-1,j) * u(i-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + &
748 (g%dxCv(i,j) * v(i,j,k) * (h(i,j,k)+h(i,j+1,k)) - &
749 g%dxCv(i,j-1)*v(i,j-1,k)*(h(i,j,k)+h(i,j-1,k))))*g%IareaT(i,j)/ &
750 (h(i,j,k) + gv%H_subroundoff)
755 call pass_var(div_xx, g%Domain, complete=.true.)
759 do j=jsq-1,jeq+2 ;
do i=is-2,ieq+1
760 div_xx_dx(i,j) = g%IdxCu(i,j)*(div_xx(i+1,j) - div_xx(i,j))
763 do j=js-2,jeq+1 ;
do i=isq-1,ieq+2
764 div_xx_dy(i,j) = g%IdyCv(i,j)*(div_xx(i,j+1) - div_xx(i,j))
768 call pass_vector(div_xx_dx, div_xx_dy, g%Domain)
773 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
774 grad_div_mag_h(i,j) = sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i-1,j)))**2 + &
775 (0.5*(div_xx_dy(i,j) + div_xx_dy(i,j-1)))**2)
778 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
779 grad_div_mag_q(i,j) = sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i,j+1)))**2 + &
780 (0.5*(div_xx_dy(i,j) + div_xx_dy(i+1,j)))**2)
785 do j=jsq-1,jeq+2 ;
do i=is-2,ieq+1
788 do j=js-2,jeq+1 ;
do i=isq-1,ieq+2
792 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
793 grad_div_mag_h(i,j) = 0.0
796 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
797 grad_div_mag_q(i,j) = 0.0
803 if (cs%use_beta_in_Leith)
then
805 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
806 beta_h(i,j) = sqrt( g%dF_dx(i,j)**2 + g%dF_dy(i,j)**2 )
808 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
809 beta_q(i,j) = sqrt( (0.25*(g%dF_dx(i,j)+g%dF_dx(i+1,j)+g%dF_dx(i,j+1)+g%dF_dx(i+1,j+1))**2) + &
810 (0.25*(g%dF_dy(i,j)+g%dF_dy(i+1,j)+g%dF_dy(i,j+1)+g%dF_dy(i+1,j+1))**2) )
813 do j=js-2,jeq+1 ;
do i=is-1,ieq+1
814 vort_xy_dx(i,j) = vort_xy_dx(i,j) + 0.5 * ( g%dF_dx(i,j) + g%dF_dx(i,j+1))
816 do j=js-1,jeq+1 ;
do i=is-2,ieq+1
817 vort_xy_dy(i,j) = vort_xy_dy(i,j) + 0.5 * ( g%dF_dy(i,j) + g%dF_dy(i+1,j))
821 if (cs%use_QG_Leith_visc)
then
824 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
825 grad_vort_mag_h_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
826 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
829 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
830 grad_vort_mag_q_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
831 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
834 call calc_qg_leith_viscosity(varmix, g, gv, h, k, div_xx_dx, div_xx_dy, &
835 vort_xy_dx, vort_xy_dy)
840 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
841 grad_vort_mag_h(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
842 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
845 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
846 grad_vort_mag_q(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
847 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
854 do j=jsq,jeq+1 ;
do i=isq,ieq+1
855 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah))
then
856 shear_mag = us%T_to_s * sqrt(sh_xx(i,j)*sh_xx(i,j) + &
857 0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
858 (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
860 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
861 if (cs%use_QG_Leith_visc)
then
862 vert_vort_mag = us%T_to_s*min(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3*grad_vort_mag_h_2d(i,j))
864 vert_vort_mag = us%T_to_s*(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j))
867 if (cs%better_bound_Ah .or. cs%better_bound_Kh)
then
868 hrat_min = min(1.0, min(h_u(i,j), h_u(i-1,j), h_v(i,j), h_v(i,j-1)) / &
869 (h(i,j,k) + h_neglect) )
873 if (cs%Laplacian)
then
876 kh = cs%Kh_bg_xx(i,j)
877 if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xx(i,j) * shear_mag )
878 if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3)
880 if (rescale_kh) kh = varmix%Res_fn_h(i,j) * kh
881 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_h(i,j)
883 if (legacy_bound) kh = min(kh, cs%Kh_Max_xx(i,j))
884 kh = max( kh, cs%Kh_bg_min )
886 kh = kh + meke%Ku(i,j) * meke_res_fn
887 if (cs%anisotropic) kh = kh + cs%Kh_aniso * ( 1. - cs%n1n2_h(i,j)**2 )
891 if (cs%better_bound_Kh)
then
892 if (kh >= hrat_min*cs%Kh_Max_xx(i,j))
then
894 kh = hrat_min*cs%Kh_Max_xx(i,j)
897 visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xx(i,j))
901 if ((cs%id_Kh_h>0) .or. find_frictwork) kh_h(i,j,k) = kh
902 if (cs%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)
904 str_xx(i,j) = -kh * sh_xx(i,j)
909 if (cs%anisotropic)
then
911 local_strain = 0.25 * ( (sh_xy(i,j) + sh_xy(i-1,j-1)) + (sh_xy(i-1,j) + sh_xy(i,j-1)) )
913 str_xx(i,j) = str_xx(i,j) - cs%Kh_aniso * cs%n1n2_h(i,j) * cs%n1n1_m_n2n2_h(i,j) * local_strain
916 if (cs%biharmonic)
then
919 ahsm = 0.0; ahlth = 0.0
920 if ((cs%Smagorinsky_Ah) .or. (cs%Leith_Ah))
then
921 if (cs%Smagorinsky_Ah)
then
922 if (cs%bound_Coriolis)
then
923 ahsm = shear_mag * (cs%Biharm_const_xx(i,j) + &
924 cs%Biharm_const2_xx(i,j)*shear_mag)
926 ahsm = cs%Biharm_const_xx(i,j) * shear_mag
929 if (cs%Leith_Ah) ahlth = cs%biharm5_const_xx(i,j) * vert_vort_mag * inv_pi5
930 ah = max(max(cs%Ah_bg_xx(i,j), ahsm), ahlth)
931 if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
932 ah = min(ah, cs%Ah_Max_xx(i,j))
934 ah = cs%Ah_bg_xx(i,j)
937 if (use_meke_au) ah = ah + meke%Au(i,j)
939 if (cs%better_bound_Ah)
then
940 ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xx(i,j))
943 if ((cs%id_Ah_h>0) .or. find_frictwork) ah_h(i,j,k) = ah
945 str_xx(i,j) = str_xx(i,j) + ah * &
946 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
947 cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(i,j-1)))
950 bhstr_xx(i,j) = ah * &
951 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
952 cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(i,j-1)))
953 bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
959 if (cs%biharmonic)
then
961 do j=js-1,jeq ;
do i=is-1,ieq
962 dvdx(i,j) = cs%DY_dxBu(i,j)*(v0(i+1,j)*g%IdyCv(i+1,j) - v0(i,j)*g%IdyCv(i,j))
963 dudy(i,j) = cs%DX_dyBu(i,j)*(u0(i,j+1)*g%IdxCu(i,j+1) - u0(i,j)*g%IdxCu(i,j))
966 if (apply_obc)
then ;
if (obc%zero_strain .or. obc%freeslip_strain)
then
967 do n=1,obc%number_of_segments
968 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
969 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= jeq))
then
970 do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
971 if (obc%zero_strain)
then
972 dvdx(i,j) = 0. ; dudy(i,j) = 0.
973 elseif (obc%freeslip_strain)
then
977 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ieq))
then
978 do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
979 if (obc%zero_strain)
then
980 dvdx(i,j) = 0. ; dudy(i,j) = 0.
981 elseif (obc%freeslip_strain)
then
992 do j=js-1,jeq ;
do i=is-1,ieq
993 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah))
then
994 shear_mag = us%T_to_s * sqrt(sh_xy(i,j)*sh_xy(i,j) + &
995 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + &
996 (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j))))
998 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
999 if (cs%use_QG_Leith_visc)
then
1000 vert_vort_mag = us%T_to_s*min(grad_vort_mag_q(i,j) + grad_div_mag_q(i,j), 3*grad_vort_mag_q_2d(i,j))
1002 vert_vort_mag = us%T_to_s*(grad_vort_mag_q(i,j) + grad_div_mag_q(i,j))
1005 h2uq = 4.0 * h_u(i,j) * h_u(i,j+1)
1006 h2vq = 4.0 * h_v(i,j) * h_v(i+1,j)
1007 hq(i,j) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
1008 ((h_u(i,j) + h_u(i,j+1)) + (h_v(i,j) + h_v(i+1,j))))
1010 if (cs%better_bound_Ah .or. cs%better_bound_Kh)
then
1011 hrat_min = min(1.0, min(h_u(i,j), h_u(i,j+1), h_v(i,j), h_v(i+1,j)) / &
1012 (hq(i,j) + h_neglect) )
1013 visc_bound_rem = 1.0
1016 if (cs%no_slip .and. (g%mask2dBu(i,j) < 0.5))
then
1017 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
1018 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0)
then
1021 hu = g%mask2dCu(i,j) * h_u(i,j) + g%mask2dCu(i,j+1) * h_u(i,j+1)
1022 hv = g%mask2dCv(i,j) * h_v(i,j) + g%mask2dCv(i+1,j) * h_v(i+1,j)
1023 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) * &
1024 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) == 0.0)
then
1030 hq(i,j) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
1031 hrat_min = min(1.0, min(hu, hv) / (hq(i,j) + h_neglect) )
1036 if (cs%Laplacian)
then
1039 kh = cs%Kh_bg_xy(i,j)
1040 if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xy(i,j) * shear_mag )
1041 if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xy(i,j) * vert_vort_mag*inv_pi3)
1043 if (rescale_kh) kh = varmix%Res_fn_q(i,j) * kh
1044 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_q(i,j)
1046 if (legacy_bound) kh = min(kh, cs%Kh_Max_xy(i,j))
1047 kh = max( kh, cs%Kh_bg_min )
1048 if (use_meke_ku)
then
1049 kh = kh + 0.25*( (meke%Ku(i,j) + meke%Ku(i+1,j+1)) + &
1050 (meke%Ku(i+1,j) + meke%Ku(i,j+1)) ) * meke_res_fn
1053 if (cs%anisotropic) kh = kh + cs%Kh_aniso * cs%n1n2_q(i,j)**2
1057 if (cs%better_bound_Kh)
then
1058 if (kh >= hrat_min*cs%Kh_Max_xy(i,j))
then
1059 visc_bound_rem = 0.0
1060 kh = hrat_min*cs%Kh_Max_xy(i,j)
1061 elseif (cs%Kh_Max_xy(i,j)>0.)
then
1063 visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xy(i,j))
1067 if (cs%id_Kh_q>0) kh_q(i,j,k) = kh
1068 if (cs%id_vort_xy_q>0) vort_xy_q(i,j,k) = vort_xy(i,j)
1070 str_xy(i,j) = -kh * sh_xy(i,j)
1075 if (cs%anisotropic)
then
1077 local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )
1079 str_xy(i,j) = str_xy(i,j) - cs%Kh_aniso * cs%n1n2_q(i,j) * cs%n1n1_m_n2n2_q(i,j) * local_strain
1082 if (cs%biharmonic)
then
1085 ahsm = 0.0 ; ahlth = 0.0
1086 if (cs%Smagorinsky_Ah .or. cs%Leith_Ah)
then
1087 if (cs%Smagorinsky_Ah)
then
1088 if (cs%bound_Coriolis)
then
1089 ahsm = shear_mag * (cs%Biharm_const_xy(i,j) + &
1090 cs%Biharm_const2_xy(i,j)*shear_mag)
1092 ahsm = cs%Biharm_const_xy(i,j) * shear_mag
1095 if (cs%Leith_Ah) ahlth = cs%Biharm5_const_xy(i,j) * vert_vort_mag * inv_pi5
1096 ah = max(max(cs%Ah_bg_xy(i,j), ahsm), ahlth)
1097 if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
1098 ah = min(ah, cs%Ah_Max_xy(i,j))
1100 ah = cs%Ah_bg_xy(i,j)
1103 if (use_meke_au)
then
1104 ah = ah + 0.25*( (meke%Au(i,j) + meke%Au(i+1,j+1)) + &
1105 (meke%Au(i+1,j) + meke%Au(i,j+1)) )
1108 if (cs%better_bound_Ah)
then
1109 ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xy(i,j))
1112 if (cs%id_Ah_q>0) ah_q(i,j,k) = ah
1114 str_xy(i,j) = str_xy(i,j) + ah * ( dvdx(i,j) + dudy(i,j) )
1117 bhstr_xy(i,j) = ah * ( dvdx(i,j) + dudy(i,j) ) * &
1118 (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
1125 if (find_frictwork)
then
1126 if (cs%Laplacian)
then
1127 if (cs%answers_2018)
then
1128 do j=js,je ;
do i=is,ie
1129 grad_vel_mag_h(i,j) = boundary_mask(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
1130 (0.25*(dvdx(i,j)+dvdx(i-1,j)+dvdx(i,j-1)+dvdx(i-1,j-1)) )**2 + &
1131 (0.25*(dudy(i,j)+dudy(i-1,j)+dudy(i,j-1)+dudy(i-1,j-1)) )**2)
1134 do j=js,je ;
do i=is,ie
1135 grad_vel_mag_h(i,j) = boundary_mask(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + &
1136 ((0.25*((dvdx(i,j) + dvdx(i-1,j-1)) + (dvdx(i-1,j) + dvdx(i,j-1))) )**2 + &
1137 (0.25*((dudy(i,j) + dudy(i-1,j-1)) + (dudy(i-1,j) + dudy(i,j-1))) )**2))
1141 do j=js,je ;
do i=is,ie
1142 grad_vel_mag_h(i,j) = 0.0
1146 if (cs%biharmonic)
then
1147 do j=js,je ;
do i=is,ie
1148 grad_d2vel_mag_h(i,j) = boundary_mask(i,j) * ((0.5*(u0(i,j) + u0(i-1,j)))**2 + &
1149 (0.5*(v0(i,j) + v0(i,j-1)))**2)
1152 do j=js,je ;
do i=is,ie
1153 grad_d2vel_mag_h(i,j) = 0.0
1157 do j=js,je ;
do i=is,ie
1159 diss_rate(i,j,k) = -us%s_to_T*kh_h(i,j,k) * grad_vel_mag_h(i,j) - &
1160 us%s_to_T*ah_h(i,j,k) * grad_d2vel_mag_h(i,j)
1162 if (
associated(meke))
then ;
if (
associated(meke%mom_src))
then
1165 max_diss_rate(i,j,k) = 2.0*meke%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
1166 frictwork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * gv%H_to_kg_m2
1167 frictworkmax(i,j,k) = -max_diss_rate(i,j,k) * h(i,j,k) * gv%H_to_kg_m2
1173 if (cs%use_GME)
then
1174 target_diss_rate_gme(i,j,k) = fwfrac * max_diss_rate(i,j,k) - diss_rate(i,j,k)
1178 frictwork_diss(i,j,k) = diss_rate(i,j,k) * h(i,j,k) * gv%H_to_kg_m2
1184 if (cs%use_GME)
then
1185 if (.not. (
associated(meke)))
call mom_error(fatal,
"MEKE must be enabled for GME to be used.")
1187 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1189 if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_h(i,j)>0) )
then
1190 gme_coeff = fwfrac*us%T_to_s*max_diss_rate(i,j,k) / grad_vel_mag_bt_h(i,j)
1193 if ((g%bathyT(i,j) < h0_gme) .and. (h0_gme > 0.0)) &
1194 gme_coeff = (g%bathyT(i,j) / h0_gme)**2 * gme_coeff
1197 gme_coeff = min(gme_coeff * boundary_mask(i,j), gme_coeff_limiter)
1200 if ((cs%id_GME_coeff_h>0) .or. find_frictwork) gme_coeff_h(i,j,k) = gme_coeff
1202 str_xx_gme(i,j) = gme_coeff * sh_xx_bt(i,j)
1206 do j=js-1,jeq ;
do i=is-1,ieq
1208 if ((max_diss_rate(i,j,k) > 0) .and. (grad_vel_mag_bt_q(i,j)>0) )
then
1210 gme_coeff = fwfrac*us%T_to_s*max_diss_rate(i,j,k) / grad_vel_mag_bt_q(i,j)
1212 if ((g%bathyT(i,j) < h0_gme) .and. (h0_gme > 0.0)) &
1213 gme_coeff = (g%bathyT(i,j) / h0_gme)**2 * gme_coeff
1217 gme_coeff = min(gme_coeff * boundary_mask(i,j), gme_coeff_limiter)
1220 if (cs%id_GME_coeff_q>0) gme_coeff_q(i,j,k) = gme_coeff
1221 str_xy_gme(i,j) = gme_coeff * sh_xy_bt(i,j)
1226 call smooth_gme(cs,g,gme_flux_h=str_xx_gme)
1227 call smooth_gme(cs,g,gme_flux_q=str_xy_gme)
1229 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1230 str_xx(i,j) = (str_xx(i,j) + str_xx_gme(i,j)) * (h(i,j,k) * cs%reduction_xx(i,j))
1233 do j=js-1,jeq ;
do i=is-1,ieq
1235 if (cs%no_slip)
then
1236 str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * cs%reduction_xy(i,j))
1238 str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
1242 if (
associated(meke%GME_snk))
then
1243 do j=js,je ;
do i=is,ie
1244 frictwork_gme(i,j,k) = us%s_to_T*gme_coeff_h(i,j,k) * h(i,j,k) * gv%H_to_kg_m2 * grad_vel_mag_bt_h(i,j)
1249 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1250 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
1253 do j=js-1,jeq ;
do i=is-1,ieq
1254 if (cs%no_slip)
then
1255 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * cs%reduction_xy(i,j))
1257 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
1264 do j=js,je ;
do i=isq,ieq
1265 diffu(i,j,k) = ((g%IdyCu(i,j)*(cs%DY2h(i,j) *str_xx(i,j) - &
1266 cs%DY2h(i+1,j)*str_xx(i+1,j)) + &
1267 g%IdxCu(i,j)*(cs%DX2q(i,j-1)*str_xy(i,j-1) - &
1268 cs%DX2q(i,j) *str_xy(i,j))) * &
1269 g%IareaCu(i,j)) / (h_u(i,j) + h_neglect)
1275 do n=1,obc%number_of_segments
1276 if (obc%segment(n)%is_E_or_W)
then
1277 i = obc%segment(n)%HI%IsdB
1278 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
1286 do j=jsq,jeq ;
do i=is,ie
1287 diffv(i,j,k) = ((g%IdyCv(i,j)*(cs%DY2q(i-1,j)*str_xy(i-1,j) - &
1288 cs%DY2q(i,j) *str_xy(i,j)) - &
1289 g%IdxCv(i,j)*(cs%DX2h(i,j) *str_xx(i,j) - &
1290 cs%DX2h(i,j+1)*str_xx(i,j+1))) * &
1291 g%IareaCv(i,j)) / (h_v(i,j) + h_neglect)
1296 do n=1,obc%number_of_segments
1297 if (obc%segment(n)%is_N_or_S)
then
1298 j = obc%segment(n)%HI%JsdB
1299 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
1306 if (find_frictwork)
then ;
do j=js,je ;
do i=is,ie
1309 frictwork(i,j,k) = us%s_to_T*gv%H_to_kg_m2 * ( &
1310 (str_xx(i,j)*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
1311 -str_xx(i,j)*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
1312 +0.25*((str_xy(i,j)*( &
1313 (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
1314 +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
1315 +str_xy(i-1,j-1)*( &
1316 (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
1317 +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
1319 (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
1320 +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
1322 (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
1323 +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
1324 enddo ;
enddo ;
endif
1329 if (find_frictwork .and.
associated(meke))
then ;
if (
associated(meke%mom_src))
then
1331 do j=js,je ;
do i=is,ie
1332 meke%mom_src(i,j) = 0.
1333 meke%GME_snk(i,j) = 0.
1336 if (meke%backscatter_Ro_c /= 0.)
then
1337 do j=js,je ;
do i=is,ie
1338 fath = 0.25*( (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1339 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))) )
1340 shear_mag = us%T_to_s * sqrt(sh_xx(i,j)*sh_xx(i,j) + &
1341 0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
1342 (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
1343 if (cs%answers_2018)
then
1344 fath = (us%s_to_T*fath)**meke%backscatter_Ro_pow
1347 shear_mag = ( ( (us%s_to_T*shear_mag)**meke%backscatter_Ro_pow ) + 1.e-30 ) &
1348 * meke%backscatter_Ro_c
1351 roscl = shear_mag / ( fath + shear_mag )
1353 if (fath <= backscat_subround*shear_mag)
then
1356 sh_f_pow = meke%backscatter_Ro_c * (shear_mag / fath)**meke%backscatter_Ro_pow
1357 roscl = sh_f_pow / (1.0 + sh_f_pow)
1360 meke%mom_src(i,j) = meke%mom_src(i,j) + us%s_to_T*gv%H_to_kg_m2 * ( &
1361 ((str_xx(i,j)-roscl*bhstr_xx(i,j))*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
1362 -(str_xx(i,j)-roscl*bhstr_xx(i,j))*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
1363 +0.25*(((str_xy(i,j)-roscl*bhstr_xy(i,j))*( &
1364 (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
1365 +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
1366 +(str_xy(i-1,j-1)-roscl*bhstr_xy(i-1,j-1))*( &
1367 (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
1368 +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
1369 +((str_xy(i-1,j)-roscl*bhstr_xy(i-1,j))*( &
1370 (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
1371 +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
1372 +(str_xy(i,j-1)-roscl*bhstr_xy(i,j-1))*( &
1373 (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
1374 +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
1377 do j=js,je ;
do i=is,ie
1379 meke%mom_src(i,j) = meke%mom_src(i,j) + max(frictwork_diss(i,j,k), frictworkmax(i,j,k))
1383 if (cs%use_GME .and.
associated(meke))
then
1384 if (
associated(meke%GME_snk))
then
1385 do j=js,je ;
do i=is,ie
1386 meke%GME_snk(i,j) = meke%GME_snk(i,j) + frictwork_gme(i,j,k)
1396 if (cs%id_diffu>0)
call post_data(cs%id_diffu, diffu, cs%diag)
1397 if (cs%id_diffv>0)
call post_data(cs%id_diffv, diffv, cs%diag)
1398 if (cs%id_FrictWork>0)
call post_data(cs%id_FrictWork, frictwork, cs%diag)
1399 if (cs%id_FrictWorkMax>0)
call post_data(cs%id_FrictWorkMax, frictworkmax, cs%diag)
1400 if (cs%id_FrictWork_diss>0)
call post_data(cs%id_FrictWork_diss, frictwork_diss, cs%diag)
1401 if (cs%id_FrictWork_GME>0)
call post_data(cs%id_FrictWork_GME, frictwork_gme, cs%diag)
1402 if (cs%id_Ah_h>0)
call post_data(cs%id_Ah_h, ah_h, cs%diag)
1403 if (cs%id_div_xx_h>0)
call post_data(cs%id_div_xx_h, div_xx_h, cs%diag)
1404 if (cs%id_vort_xy_q>0)
call post_data(cs%id_vort_xy_q, vort_xy_q, cs%diag)
1405 if (cs%id_Ah_q>0)
call post_data(cs%id_Ah_q, ah_q, cs%diag)
1406 if (cs%id_Kh_h>0)
call post_data(cs%id_Kh_h, kh_h, cs%diag)
1407 if (cs%id_Kh_q>0)
call post_data(cs%id_Kh_q, kh_q, cs%diag)
1408 if (cs%id_GME_coeff_h > 0)
call post_data(cs%id_GME_coeff_h, gme_coeff_h, cs%diag)
1409 if (cs%id_GME_coeff_q > 0)
call post_data(cs%id_GME_coeff_q, gme_coeff_q, cs%diag)
1411 if (cs%id_FrictWorkIntz > 0)
then
1413 do i=is,ie ; frictworkintz(i,j) = frictwork(i,j,1) ;
enddo
1414 do k=2,nz ;
do i=is,ie
1415 frictworkintz(i,j) = frictworkintz(i,j) + frictwork(i,j,k)
1418 call post_data(cs%id_FrictWorkIntz, frictworkintz, cs%diag)