22 implicit none ;
private
24 #include <MOM_memory.h>
26 public horizontal_viscosity, hor_visc_init, hor_visc_end
41 logical :: better_bound_kh
45 logical :: better_bound_ah
51 logical :: smagorinsky_kh
53 logical :: smagorinsky_ah
57 logical :: modified_leith
59 logical :: use_beta_in_leith
62 logical :: use_qg_leith_visc
64 logical :: bound_coriolis
67 logical :: use_kh_bg_2d
70 logical :: use_land_mask
75 logical :: anisotropic
77 logical :: dynamic_aniso
79 logical :: res_scale_meke
82 logical :: answers_2018
86 real allocable_,
dimension(NIMEM_,NJMEM_) :: kh_bg_xx
90 real allocable_,
dimension(NIMEM_,NJMEM_) :: kh_bg_2d
94 real allocable_,
dimension(NIMEM_,NJMEM_) :: ah_bg_xx
98 real allocable_,
dimension(NIMEM_,NJMEM_) :: biharm5_const2_xx
103 real allocable_,
dimension(NIMEM_,NJMEM_) :: reduction_xx
106 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
107 kh_max_xx, & !< The maximum permitted Laplacian viscosity [m2 T-1 ~> m2 s-1].
108 ah_max_xx, & !< The maximum permitted biharmonic viscosity [m4 T-1 ~> m4 s-1].
109 n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points
112 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: kh_bg_xy
116 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: ah_bg_xy
120 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: biharm5_const2_xy
125 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: reduction_xy
128 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
129 kh_max_xy, & !< The maximum permitted Laplacian viscosity [m2 T-1 ~> m2 s-1].
130 ah_max_xy, & !< The maximum permitted biharmonic viscosity [m4 T-1 ~> m4 s-1].
131 n1n2_q, & !< Factor n1*n2 in the anisotropic direction tensor at q-points
134 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
135 dx2h, & !< Pre-calculated dx^2 at h points [m2]
136 dy2h, & !< Pre-calculated dy^2 at h points [m2]
137 dx_dyt, & !< Pre-calculated dx/dy at h points [nondim]
139 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
140 dx2q, & !< Pre-calculated dx^2 at q points [m2]
141 dy2q, & !< Pre-calculated dy^2 at q points [m2]
142 dx_dybu, & !< Pre-calculated dx/dy at q points [nondim]
144 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: &
145 idx2dycu, & !< 1/(dx^2 dy) at u points [m-3]
147 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: &
148 idx2dycv, & !< 1/(dx^2 dy) at v points [m-3]
153 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
154 laplac2_const_xx, & !< Laplacian metric-dependent constants [m2]
155 biharm5_const_xx, & !< Biharmonic metric-dependent constants [m5]
156 laplac3_const_xx, & !< Laplacian metric-dependent constants [m3]
157 biharm_const_xx, & !< Biharmonic metric-dependent constants [m4]
160 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
161 laplac2_const_xy, & !< Laplacian metric-dependent constants [m2]
162 biharm5_const_xy, & !< Biharmonic metric-dependent constants [m5]
163 laplac3_const_xy, & !< Laplacian metric-dependent constants [m3]
164 biharm_const_xy, & !< Biharmonic metric-dependent constants [m4]
171 integer :: id_diffu = -1, id_diffv = -1
172 integer :: id_ah_h = -1, id_ah_q = -1
173 integer :: id_kh_h = -1, id_kh_q = -1
174 integer :: id_gme_coeff_h = -1, id_gme_coeff_q = -1
175 integer :: id_vort_xy_q = -1, id_div_xx_h = -1
176 integer :: id_frictwork = -1, id_frictworkintz = -1
177 integer :: id_frictworkmax = -1
178 integer :: id_frictwork_diss = -1, id_frictwork_gme = -1
198 subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
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)), &
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)
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))
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))
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)
1421 end subroutine horizontal_viscosity
1426 subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
1427 type(time_type),
intent(in) :: time
1432 type(
diag_ctrl),
target,
intent(inout) :: diag
1436 real,
dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v
1437 real,
dimension(SZI_(G),SZJB_(G)) :: v0u, v0v
1448 real :: boundcorconst
1454 real :: kh_vel_scale
1455 real :: ah_vel_scale
1456 real :: ah_time_scale
1457 real :: smag_lap_const
1458 real :: smag_bi_const
1459 real :: leith_lap_const
1460 real :: leith_bi_const
1465 real :: bound_cor_vel
1469 real :: kh_pwr_of_sine
1470 logical :: bound_cor_def
1476 logical :: default_2018_answers
1477 character(len=64) :: inputdir, filename
1480 real :: aniso_grid_dir(2)
1481 integer :: aniso_mode
1482 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
1483 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1487 #include "version_variable.h"
1488 character(len=40) :: mdl =
"MOM_hor_visc"
1490 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1491 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1492 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1493 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1495 if (
associated(cs))
then
1496 call mom_error(warning,
"hor_visc_init called with an associated "// &
1497 "control structure.")
1509 cs%bound_Kh = .false. ; cs%better_bound_Kh = .false. ; cs%Smagorinsky_Kh = .false. ; cs%Leith_Kh = .false.
1510 cs%bound_Ah = .false. ; cs%better_bound_Ah = .false. ; cs%Smagorinsky_Ah = .false. ; cs%Leith_Ah = .false.
1511 cs%use_QG_Leith_visc = .false.
1512 cs%bound_Coriolis = .false.
1513 cs%Modified_Leith = .false.
1514 cs%anisotropic = .false.
1515 cs%dynamic_aniso = .false.
1521 call get_param(param_file, mdl,
"GET_ALL_PARAMS", get_all, default=.false.)
1523 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1524 "This sets the default value for the various _2018_ANSWERS parameters.", &
1526 call get_param(param_file, mdl,
"HOR_VISC_2018_ANSWERS", cs%answers_2018, &
1527 "If true, use the order of arithmetic and expressions that recover the "//&
1528 "answers from the end of 2018. Otherwise, use updated and more robust "//&
1529 "forms of the same expressions.", default=default_2018_answers)
1530 call get_param(param_file, mdl,
"LAPLACIAN", cs%Laplacian, &
1531 "If true, use a Laplacian horizontal viscosity.", &
1533 if (cs%Laplacian .or. get_all)
then
1534 call get_param(param_file, mdl,
"KH", kh, &
1535 "The background Laplacian horizontal viscosity.", &
1536 units =
"m2 s-1", default=0.0, scale=us%T_to_s)
1537 call get_param(param_file, mdl,
"KH_BG_MIN", cs%Kh_bg_min, &
1538 "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
1539 units =
"m2 s-1", default=0.0, scale=us%T_to_s)
1540 call get_param(param_file, mdl,
"KH_VEL_SCALE", kh_vel_scale, &
1541 "The velocity scale which is multiplied by the grid "//&
1542 "spacing to calculate the Laplacian viscosity. "//&
1543 "The final viscosity is the largest of this scaled "//&
1544 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
1545 units=
"m s-1", default=0.0, scale=us%T_to_s)
1546 call get_param(param_file, mdl,
"KH_SIN_LAT", kh_sin_lat, &
1547 "The amplitude of a latitudinally-dependent background "//&
1548 "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &
1549 units =
"m2 s-1", default=0.0, scale=us%T_to_s)
1550 if (kh_sin_lat>0. .or. get_all) &
1551 call get_param(param_file, mdl,
"KH_PWR_OF_SINE", kh_pwr_of_sine, &
1552 "The power used to raise SIN(LAT) when using a latitudinally "//&
1553 "dependent background viscosity.", &
1554 units =
"nondim", default=4.0)
1556 call get_param(param_file, mdl,
"SMAGORINSKY_KH", cs%Smagorinsky_Kh, &
1557 "If true, use a Smagorinsky nonlinear eddy viscosity.", &
1559 if (cs%Smagorinsky_Kh .or. get_all) &
1560 call get_param(param_file, mdl,
"SMAG_LAP_CONST", smag_lap_const, &
1561 "The nondimensional Laplacian Smagorinsky constant, "//&
1562 "often 0.15.", units=
"nondim", default=0.0, &
1563 fail_if_missing = cs%Smagorinsky_Kh)
1565 call get_param(param_file, mdl,
"LEITH_KH", cs%Leith_Kh, &
1566 "If true, use a Leith nonlinear eddy viscosity.", &
1569 call get_param(param_file, mdl,
"MODIFIED_LEITH", cs%Modified_Leith, &
1570 "If true, add a term to Leith viscosity which is "//&
1571 "proportional to the gradient of divergence.", &
1573 call get_param(param_file, mdl,
"RES_SCALE_MEKE_VISC", cs%res_scale_MEKE, &
1574 "If true, the viscosity contribution from MEKE is scaled by "//&
1575 "the resolution function.", default=.false.)
1577 if (cs%Leith_Kh .or. get_all)
then
1578 call get_param(param_file, mdl,
"LEITH_LAP_CONST", leith_lap_const, &
1579 "The nondimensional Laplacian Leith constant, "//&
1580 "often set to 1.0", units=
"nondim", default=0.0, &
1581 fail_if_missing = cs%Leith_Kh)
1582 call get_param(param_file, mdl,
"USE_QG_LEITH_VISC", cs%use_QG_Leith_visc, &
1583 "If true, use QG Leith nonlinear eddy viscosity.", &
1585 if (cs%use_QG_Leith_visc .and. .not. cs%Leith_Kh)
call mom_error(fatal, &
1586 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1587 "LEITH_KH must be True when USE_QG_LEITH_VISC=True.")
1589 if (cs%Leith_Kh .or. cs%Leith_Ah .or. get_all)
then
1590 call get_param(param_file, mdl,
"USE_BETA_IN_LEITH", cs%use_beta_in_Leith, &
1591 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1592 default=cs%Leith_Kh)
1593 call get_param(param_file, mdl,
"MODIFIED_LEITH", cs%modified_Leith, &
1594 "If true, add a term to Leith viscosity which is \n"//&
1595 "proportional to the gradient of divergence.", &
1598 call get_param(param_file, mdl,
"BOUND_KH", cs%bound_Kh, &
1599 "If true, the Laplacian coefficient is locally limited "//&
1600 "to be stable.", default=.true.)
1601 call get_param(param_file, mdl,
"BETTER_BOUND_KH", cs%better_bound_Kh, &
1602 "If true, the Laplacian coefficient is locally limited "//&
1603 "to be stable with a better bounding than just BOUND_KH.", &
1604 default=cs%bound_Kh)
1605 call get_param(param_file, mdl,
"ANISOTROPIC_VISCOSITY", cs%anisotropic, &
1606 "If true, allow anistropic viscosity in the Laplacian "//&
1607 "horizontal viscosity.", default=.false.)
1609 if (cs%anisotropic .or. get_all)
then
1610 call get_param(param_file, mdl,
"KH_ANISO", cs%Kh_aniso, &
1611 "The background Laplacian anisotropic horizontal viscosity.", &
1612 units =
"m2 s-1", default=0.0, scale=us%T_to_s)
1613 call get_param(param_file, mdl,
"ANISOTROPIC_MODE", aniso_mode, &
1614 "Selects the mode for setting the direction of anistropy.\n"//&
1615 "\t 0 - Points along the grid i-direction.\n"//&
1616 "\t 1 - Points towards East.\n"//&
1617 "\t 2 - Points along the flow direction, U/|U|.", &
1619 select case (aniso_mode)
1621 call get_param(param_file, mdl,
"ANISO_GRID_DIR", aniso_grid_dir, &
1622 "The vector pointing in the direction of anistropy for "//&
1623 "horizont viscosity. n1,n2 are the i,j components relative "//&
1624 "to the grid.", units =
"nondim", fail_if_missing=.true.)
1626 call get_param(param_file, mdl,
"ANISO_GRID_DIR", aniso_grid_dir, &
1627 "The vector pointing in the direction of anistropy for "//&
1628 "horizont viscosity. n1,n2 are the i,j components relative "//&
1629 "to the spherical coordinates.", units =
"nondim", fail_if_missing=.true.)
1633 call get_param(param_file, mdl,
"BIHARMONIC", cs%biharmonic, &
1634 "If true, use a biharmonic horizontal viscosity. "//&
1635 "BIHARMONIC may be used with LAPLACIAN.", &
1637 if (cs%biharmonic .or. get_all)
then
1638 call get_param(param_file, mdl,
"AH", ah, &
1639 "The background biharmonic horizontal viscosity.", &
1640 units =
"m4 s-1", default=0.0, scale=us%T_to_s)
1641 call get_param(param_file, mdl,
"AH_VEL_SCALE", ah_vel_scale, &
1642 "The velocity scale which is multiplied by the cube of "//&
1643 "the grid spacing to calculate the biharmonic viscosity. "//&
1644 "The final viscosity is the largest of this scaled "//&
1645 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
1646 units=
"m s-1", default=0.0, scale=us%T_to_s)
1647 call get_param(param_file, mdl,
"AH_TIME_SCALE", ah_time_scale, &
1648 "A time scale whose inverse is multiplied by the fourth "//&
1649 "power of the grid spacing to calculate biharmonic viscosity. "//&
1650 "The final viscosity is the largest of all viscosity "//&
1651 "formulations in use. 0.0 means that it's not used.", &
1652 units=
"s", default=0.0, scale=us%s_to_T)
1653 call get_param(param_file, mdl,
"SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
1654 "If true, use a biharmonic Smagorinsky nonlinear eddy "//&
1655 "viscosity.", default=.false.)
1656 call get_param(param_file, mdl,
"LEITH_AH", cs%Leith_Ah, &
1657 "If true, use a biharmonic Leith nonlinear eddy "//&
1658 "viscosity.", default=.false.)
1660 call get_param(param_file, mdl,
"BOUND_AH", cs%bound_Ah, &
1661 "If true, the biharmonic coefficient is locally limited "//&
1662 "to be stable.", default=.true.)
1663 call get_param(param_file, mdl,
"BETTER_BOUND_AH", cs%better_bound_Ah, &
1664 "If true, the biharmonic coefficient is locally limited "//&
1665 "to be stable with a better bounding than just BOUND_AH.", &
1666 default=cs%bound_Ah)
1668 if (cs%Smagorinsky_Ah .or. get_all)
then
1669 call get_param(param_file, mdl,
"SMAG_BI_CONST",smag_bi_const, &
1670 "The nondimensional biharmonic Smagorinsky constant, "//&
1671 "typically 0.015 - 0.06.", units=
"nondim", default=0.0, &
1672 fail_if_missing = cs%Smagorinsky_Ah)
1674 call get_param(param_file, mdl,
"BOUND_CORIOLIS", bound_cor_def, default=.false.)
1675 call get_param(param_file, mdl,
"BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &
1676 "If true use a viscosity that increases with the square "//&
1677 "of the velocity shears, so that the resulting viscous "//&
1678 "drag is of comparable magnitude to the Coriolis terms "//&
1679 "when the velocity differences between adjacent grid "//&
1680 "points is 0.5*BOUND_CORIOLIS_VEL. The default is the "//&
1681 "value of BOUND_CORIOLIS (or false).", default=bound_cor_def)
1682 if (cs%bound_Coriolis .or. get_all)
then
1683 call get_param(param_file, mdl,
"MAXVEL", maxvel, default=3.0e8)
1684 bound_cor_vel = maxvel
1685 call get_param(param_file, mdl,
"BOUND_CORIOLIS_VEL", bound_cor_vel, &
1686 "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&
1687 "the biharmonic drag to have comparable magnitude to the "//&
1688 "Coriolis acceleration. The default is set by MAXVEL.", &
1689 units=
"m s-1", default=maxvel, scale=us%m_s_to_L_T)
1693 if (cs%Leith_Ah .or. get_all) &
1694 call get_param(param_file, mdl,
"LEITH_BI_CONST", leith_bi_const, &
1695 "The nondimensional biharmonic Leith constant, "//&
1696 "typical values are thus far undetermined.", units=
"nondim", default=0.0, &
1697 fail_if_missing = cs%Leith_Ah)
1701 call get_param(param_file, mdl,
"USE_LAND_MASK_FOR_HVISC", cs%use_land_mask, &
1702 "If true, use Use the land mask for the computation of thicknesses "//&
1703 "at velocity locations. This eliminates the dependence on arbitrary "//&
1704 "values over land or outside of the domain. Default is False in order to "//&
1705 "maintain answers with legacy experiments but should be changed to True "//&
1706 "for new experiments.", default=.false.)
1708 if (cs%better_bound_Ah .or. cs%better_bound_Kh .or. get_all) &
1709 call get_param(param_file, mdl,
"HORVISC_BOUND_COEF", cs%bound_coef, &
1710 "The nondimensional coefficient of the ratio of the "//&
1711 "viscosity bounds to the theoretical maximum for "//&
1712 "stability without considering other terms.", units=
"nondim", &
1715 call get_param(param_file, mdl,
"NOSLIP", cs%no_slip, &
1716 "If true, no slip boundary conditions are used; otherwise "//&
1717 "free slip boundary conditions are assumed. The "//&
1718 "implementation of the free slip BCs on a C-grid is much "//&
1719 "cleaner than the no slip BCs. The use of free slip BCs "//&
1720 "is strongly encouraged, and no slip BCs are not used with "//&
1721 "the biharmonic viscosity.", default=.false.)
1723 call get_param(param_file, mdl,
"USE_KH_BG_2D", cs%use_Kh_bg_2d, &
1724 "If true, read a file containing 2-d background harmonic "//&
1725 "viscosities. The final viscosity is the maximum of the other "//&
1726 "terms and this background value.", default=.false.)
1728 call get_param(param_file, mdl,
"USE_GME", cs%use_GME, &
1729 "If true, use the GM+E backscatter scheme in association \n"//&
1730 "with the Gent and McWilliams parameterization.", default=.false.)
1732 if (cs%use_GME)
then
1733 call get_param(param_file, mdl,
"SPLIT", split, &
1734 "Use the split time stepping if true.", default=.true., &
1736 if (.not. split)
call mom_error(fatal,
"ERROR: Currently, USE_GME = True "// &
1737 "cannot be used with SPLIT=False.")
1740 if (cs%bound_Kh .or. cs%bound_Ah .or. cs%better_bound_Kh .or. cs%better_bound_Ah) &
1741 call get_param(param_file, mdl,
"DT", dt, &
1742 "The (baroclinic) dynamics time step.", units=
"s", scale=us%s_to_T, &
1743 fail_if_missing=.true.)
1745 if (cs%no_slip .and. cs%biharmonic) &
1746 call mom_error(fatal,
"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
1747 "at the same time in MOM.")
1749 if (.not.(cs%Laplacian .or. cs%biharmonic))
then
1751 if ( max(g%domain%niglobal, g%domain%njglobal)>2 )
call mom_error(warning, &
1752 "hor_visc_init: It is usually a very bad idea not to use either "//&
1753 "LAPLACIAN or BIHARMONIC viscosity.")
1757 deg2rad = atan(1.0) / 45.
1759 alloc_(cs%dx2h(isd:ied,jsd:jed)) ; cs%dx2h(:,:) = 0.0
1760 alloc_(cs%dy2h(isd:ied,jsd:jed)) ; cs%dy2h(:,:) = 0.0
1761 alloc_(cs%dx2q(isdb:iedb,jsdb:jedb)) ; cs%dx2q(:,:) = 0.0
1762 alloc_(cs%dy2q(isdb:iedb,jsdb:jedb)) ; cs%dy2q(:,:) = 0.0
1763 alloc_(cs%dx_dyT(isd:ied,jsd:jed)) ; cs%dx_dyT(:,:) = 0.0
1764 alloc_(cs%dy_dxT(isd:ied,jsd:jed)) ; cs%dy_dxT(:,:) = 0.0
1765 alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0
1766 alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0
1768 if (cs%Laplacian)
then
1769 alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed)) ; cs%Kh_bg_xx(:,:) = 0.0
1770 alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0
1771 if (cs%bound_Kh .or. cs%better_bound_Kh)
then
1772 alloc_(cs%Kh_Max_xx(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xx(:,:) = 0.0
1773 alloc_(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xy(:,:) = 0.0
1775 if (cs%Smagorinsky_Kh)
then
1776 alloc_(cs%Laplac2_const_xx(isd:ied,jsd:jed)) ; cs%Laplac2_const_xx(:,:) = 0.0
1777 alloc_(cs%Laplac2_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac2_const_xy(:,:) = 0.0
1779 if (cs%Leith_Kh)
then
1780 alloc_(cs%Laplac3_const_xx(isd:ied,jsd:jed)) ; cs%Laplac3_const_xx(:,:) = 0.0
1781 alloc_(cs%Laplac3_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac3_const_xy(:,:) = 0.0
1784 alloc_(cs%reduction_xx(isd:ied,jsd:jed)) ; cs%reduction_xx(:,:) = 0.0
1785 alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0
1787 if (cs%anisotropic)
then
1788 alloc_(cs%n1n2_h(isd:ied,jsd:jed)) ; cs%n1n2_h(:,:) = 0.0
1789 alloc_(cs%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; cs%n1n1_m_n2n2_h(:,:) = 0.0
1790 alloc_(cs%n1n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n2_q(:,:) = 0.0
1791 alloc_(cs%n1n1_m_n2n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n1_m_n2n2_q(:,:) = 0.0
1792 select case (aniso_mode)
1794 call align_aniso_tensor_to_grid(cs, aniso_grid_dir(1), aniso_grid_dir(2))
1798 cs%dynamic_aniso = .true.
1800 call mom_error(fatal,
"MOM_hor_visc: "//&
1801 "Runtime parameter ANISOTROPIC_MODE is out of range.")
1805 if (cs%use_Kh_bg_2d)
then
1806 alloc_(cs%Kh_bg_2d(isd:ied,jsd:jed)) ; cs%Kh_bg_2d(:,:) = 0.0
1807 call get_param(param_file, mdl,
"KH_BG_2D_FILENAME", filename, &
1808 'The filename containing a 2d map of "Kh".', &
1809 default=
'KH_background_2d.nc')
1810 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1811 inputdir = slasher(inputdir)
1812 call mom_read_data(trim(inputdir)//trim(filename),
'Kh', cs%Kh_bg_2d, &
1813 g%domain, timelevel=1, scale=us%T_to_s)
1814 call pass_var(cs%Kh_bg_2d, g%domain)
1817 if (cs%biharmonic)
then
1818 alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0
1819 alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0
1820 alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed)) ; cs%Idxdy2u(:,:) = 0.0
1821 alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb)) ; cs%Idxdy2v(:,:) = 0.0
1823 alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed)) ; cs%Ah_bg_xx(:,:) = 0.0
1824 alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0
1825 if (cs%bound_Ah .or. cs%better_bound_Ah)
then
1826 alloc_(cs%Ah_Max_xx(isd:ied,jsd:jed)) ; cs%Ah_Max_xx(:,:) = 0.0
1827 alloc_(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_Max_xy(:,:) = 0.0
1829 if (cs%Smagorinsky_Ah)
then
1830 alloc_(cs%Biharm_const_xx(isd:ied,jsd:jed)) ; cs%Biharm_const_xx(:,:) = 0.0
1831 alloc_(cs%Biharm_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const_xy(:,:) = 0.0
1832 if (cs%bound_Coriolis)
then
1833 alloc_(cs%Biharm_const2_xx(isd:ied,jsd:jed)) ; cs%Biharm_const2_xx(:,:) = 0.0
1834 alloc_(cs%Biharm_const2_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const2_xy(:,:) = 0.0
1837 if (cs%Leith_Ah)
then
1838 alloc_(cs%biharm5_const_xx(isd:ied,jsd:jed)) ; cs%biharm5_const_xx(:,:) = 0.0
1839 alloc_(cs%biharm5_const_xy(isdb:iedb,jsdb:jedb)) ; cs%biharm5_const_xy(:,:) = 0.0
1843 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
1844 cs%DX2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%DY2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)
1845 cs%DX_dyBu(i,j) = g%dxBu(i,j)*g%IdyBu(i,j) ; cs%DY_dxBu(i,j) = g%dyBu(i,j)*g%IdxBu(i,j)
1847 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
1848 cs%DX2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%DY2h(i,j) = g%dyT(i,j)*g%dyT(i,j)
1849 cs%DX_dyT(i,j) = g%dxT(i,j)*g%IdyT(i,j) ; cs%DY_dxT(i,j) = g%dyT(i,j)*g%IdxT(i,j)
1852 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1853 cs%reduction_xx(i,j) = 1.0
1854 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1855 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &
1856 cs%reduction_xx(i,j) = g%dy_Cu(i,j) / g%dyCu(i,j)
1857 if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &
1858 (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &
1859 cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / g%dyCu(i-1,j)
1860 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1861 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &
1862 cs%reduction_xx(i,j) = g%dx_Cv(i,j) / g%dxCv(i,j)
1863 if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &
1864 (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &
1865 cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / g%dxCv(i,j-1)
1868 do j=js-1,jeq ;
do i=is-1,ieq
1869 cs%reduction_xy(i,j) = 1.0
1870 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1871 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &
1872 cs%reduction_xy(i,j) = g%dy_Cu(i,j) / g%dyCu(i,j)
1873 if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &
1874 (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &
1875 cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / g%dyCu(i,j+1)
1876 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1877 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &
1878 cs%reduction_xy(i,j) = g%dx_Cv(i,j) / g%dxCv(i,j)
1879 if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &
1880 (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &
1881 cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / g%dxCv(i+1,j)
1884 if (cs%Laplacian)
then
1887 if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)
1890 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1892 grid_sp_h2 = (2.0*cs%DX2h(i,j)*cs%DY2h(i,j)) / (cs%DX2h(i,j) + cs%DY2h(i,j))
1893 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1894 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xx(i,j) = smag_lap_const * grid_sp_h2
1895 if (cs%Leith_Kh) cs%Laplac3_const_xx(i,j) = leith_lap_const * grid_sp_h3
1897 cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))
1900 if (cs%use_Kh_bg_2d) cs%Kh_bg_xx(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xx(i,j))
1903 if (kh_sin_lat>0.)
then
1904 slat_fn = abs( sin( deg2rad * g%geoLatT(i,j) ) ) ** kh_pwr_of_sine
1905 cs%Kh_bg_xx(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xx(i,j))
1908 if (cs%bound_Kh .and. .not.cs%better_bound_Kh)
then
1910 cs%Kh_Max_xx(i,j) = kh_limit * grid_sp_h2
1911 cs%Kh_bg_xx(i,j) = min(cs%Kh_bg_xx(i,j), cs%Kh_Max_xx(i,j))
1916 do j=js-1,jeq ;
do i=is-1,ieq
1918 grid_sp_q2 = (2.0*cs%DX2q(i,j)*cs%DY2q(i,j)) / (cs%DX2q(i,j) + cs%DY2q(i,j))
1919 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1920 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xy(i,j) = smag_lap_const * grid_sp_q2
1921 if (cs%Leith_Kh) cs%Laplac3_const_xy(i,j) = leith_lap_const * grid_sp_q3
1923 cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))
1927 if (cs%use_Kh_bg_2d) cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xy(i,j))
1930 if (kh_sin_lat>0.)
then
1931 slat_fn = abs( sin( deg2rad * g%geoLatBu(i,j) ) ) ** kh_pwr_of_sine
1932 cs%Kh_bg_xy(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xy(i,j))
1935 if (cs%bound_Kh .and. .not.cs%better_bound_Kh)
then
1937 cs%Kh_Max_xy(i,j) = kh_limit * grid_sp_q2
1938 cs%Kh_bg_xy(i,j) = min(cs%Kh_bg_xy(i,j), cs%Kh_Max_xy(i,j))
1943 if (cs%biharmonic)
then
1945 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
1946 cs%IDX2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)
1947 cs%IDXDY2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))
1949 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
1950 cs%IDX2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)
1951 cs%IDXDY2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))
1954 cs%Ah_bg_xy(:,:) = 0.0
1957 if (cs%better_bound_Ah .or. cs%bound_Ah) ah_limit = 0.3 / (dt*64.0)
1958 if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &
1959 boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))
1960 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1961 grid_sp_h2 = (2.0*cs%DX2h(i,j)*cs%DY2h(i,j)) / (cs%DX2h(i,j)+cs%DY2h(i,j))
1962 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1964 if (cs%Smagorinsky_Ah)
then
1965 cs%Biharm_const_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)
1966 if (cs%bound_Coriolis)
then
1967 fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &
1968 abs(g%CoriolisBu(i-1,j)), abs(g%CoriolisBu(i,j)))
1969 cs%Biharm_const2_xx(i,j) = us%m_to_L**2*(grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
1970 (fmax * boundcorconst)
1973 if (cs%Leith_Ah)
then
1974 cs%biharm5_const_xx(i,j) = leith_bi_const * (grid_sp_h3 * grid_sp_h2)
1976 cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
1977 if (ah_time_scale > 0.) cs%Ah_bg_xx(i,j) = &
1978 max(cs%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / ah_time_scale)
1979 if (cs%bound_Ah .and. .not.cs%better_bound_Ah)
then
1980 cs%Ah_Max_xx(i,j) = ah_limit * (grid_sp_h2 * grid_sp_h2)
1981 cs%Ah_bg_xx(i,j) = min(cs%Ah_bg_xx(i,j), cs%Ah_Max_xx(i,j))
1984 do j=js-1,jeq ;
do i=is-1,ieq
1985 grid_sp_q2 = (2.0*cs%DX2q(i,j)*cs%DY2q(i,j)) / (cs%DX2q(i,j)+cs%DY2q(i,j))
1986 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1988 if (cs%Smagorinsky_Ah)
then
1989 cs%Biharm_const_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)
1990 if (cs%bound_Coriolis)
then
1991 cs%Biharm_const2_xy(i,j) = us%m_to_L**2*(grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
1992 (abs(g%CoriolisBu(i,j)) * boundcorconst)
1995 if (cs%Leith_Ah)
then
1996 cs%biharm5_const_xy(i,j) = leith_bi_const * (grid_sp_q3 * grid_sp_q2)
1999 cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
2000 if (ah_time_scale > 0.) cs%Ah_bg_xy(i,j) = &
2001 max(cs%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / ah_time_scale)
2002 if (cs%bound_Ah .and. .not.cs%better_bound_Ah)
then
2003 cs%Ah_Max_xy(i,j) = ah_limit * (grid_sp_q2 * grid_sp_q2)
2004 cs%Ah_bg_xy(i,j) = min(cs%Ah_bg_xy(i,j), cs%Ah_Max_xy(i,j))
2010 if (cs%Laplacian .and. cs%better_bound_Kh)
then
2012 do j=jsq,jeq+1 ;
do i=isq,ieq+1
2014 (cs%DY2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &
2015 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
2016 (cs%DX2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &
2017 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
2018 cs%Kh_Max_xx(i,j) = 0.0
2020 cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom
2022 do j=js-1,jeq ;
do i=is-1,ieq
2024 (cs%DX2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &
2025 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
2026 (cs%DY2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &
2027 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
2028 cs%Kh_Max_xy(i,j) = 0.0
2030 cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom
2036 if (cs%biharmonic .and. cs%better_bound_Ah)
then
2038 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
2039 u0u(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*cs%DY_dxT(i+1,j)*(g%IdyCu(i+1,j) + g%IdyCu(i,j)) + &
2040 cs%DY2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) + &
2041 cs%IDX2dyCu(i,j)*(cs%DX2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
2042 cs%DX2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1)) )
2044 u0v(i,j) = cs%IDXDY2u(i,j)*(cs%DY2h(i+1,j)*cs%DX_dyT(i+1,j)*(g%IdxCv(i+1,j) + g%IdxCv(i+1,j-1)) + &
2045 cs%DY2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) + &
2046 cs%IDX2dyCu(i,j)*(cs%DX2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) + &
2047 cs%DX2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1)) )
2049 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
2050 v0u(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
2051 cs%DY2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j)) ) + &
2052 cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*cs%DY_dxT(i,j+1)*(g%IdyCu(i,j+1) + g%IdyCu(i-1,j+1)) + &
2053 cs%DX2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) )
2055 v0v(i,j) = cs%IDXDY2v(i,j)*(cs%DY2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) + &
2056 cs%DY2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j)) ) + &
2057 cs%IDX2dyCv(i,j)*(cs%DX2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j)) + &
2058 cs%DX2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) )
2061 do j=jsq,jeq+1 ;
do i=isq,ieq+1
2064 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0u(i,j) + g%IdyCu(i-1,j)*u0u(i-1,j)) + &
2065 cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0u(i,j) + g%IdxCv(i,j-1)*v0u(i,j-1))) * &
2066 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
2068 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0v(i,j) + g%IdyCu(i-1,j)*u0v(i-1,j)) + &
2069 cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0v(i,j) + g%IdxCv(i,j-1)*v0v(i,j-1))) * &
2070 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
2071 cs%Ah_Max_xx(i,j) = 0.0
2073 cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom
2076 do j=js-1,jeq ;
do i=is-1,ieq
2079 (cs%DX_dyBu(i,j)*(u0u(i,j+1)*g%IdxCu(i,j+1) + u0u(i,j)*g%IdxCu(i,j)) + &
2080 cs%DY_dxBu(i,j)*(v0u(i+1,j)*g%IdyCv(i+1,j) + v0u(i,j)*g%IdyCv(i,j))) * &
2081 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
2083 (cs%DX_dyBu(i,j)*(u0v(i,j+1)*g%IdxCu(i,j+1) + u0v(i,j)*g%IdxCu(i,j)) + &
2084 cs%DY_dxBu(i,j)*(v0v(i+1,j)*g%IdyCv(i+1,j) + v0v(i,j)*g%IdyCv(i,j))) * &
2085 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
2086 cs%Ah_Max_xy(i,j) = 0.0
2088 cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom
2094 cs%id_diffu = register_diag_field(
'ocean_model',
'diffu', diag%axesCuL, time, &
2095 'Zonal Acceleration from Horizontal Viscosity',
'm s-2', conversion=us%s_to_T)
2097 cs%id_diffv = register_diag_field(
'ocean_model',
'diffv', diag%axesCvL, time, &
2098 'Meridional Acceleration from Horizontal Viscosity',
'm s-2', conversion=us%s_to_T)
2100 if (cs%biharmonic)
then
2101 cs%id_Ah_h = register_diag_field(
'ocean_model',
'Ahh', diag%axesTL, time, &
2102 'Biharmonic Horizontal Viscosity at h Points',
'm4 s-1', conversion=us%s_to_T, &
2103 cmor_field_name=
'difmxybo', &
2104 cmor_long_name=
'Ocean lateral biharmonic viscosity', &
2105 cmor_standard_name=
'ocean_momentum_xy_biharmonic_diffusivity')
2107 cs%id_Ah_q = register_diag_field(
'ocean_model',
'Ahq', diag%axesBL, time, &
2108 'Biharmonic Horizontal Viscosity at q Points',
'm4 s-1', conversion=us%s_to_T)
2111 if (cs%Laplacian)
then
2112 cs%id_Kh_h = register_diag_field(
'ocean_model',
'Khh', diag%axesTL, time, &
2113 'Laplacian Horizontal Viscosity at h Points',
'm2 s-1', conversion=us%s_to_T, &
2114 cmor_field_name=
'difmxylo', &
2115 cmor_long_name=
'Ocean lateral Laplacian viscosity', &
2116 cmor_standard_name=
'ocean_momentum_xy_laplacian_diffusivity')
2118 cs%id_Kh_q = register_diag_field(
'ocean_model',
'Khq', diag%axesBL, time, &
2119 'Laplacian Horizontal Viscosity at q Points',
'm2 s-1', conversion=us%s_to_T)
2121 if (cs%Leith_Kh)
then
2122 cs%id_vort_xy_q = register_diag_field(
'ocean_model',
'vort_xy_q', diag%axesBL, time, &
2123 'Vertical vorticity at q Points',
's-1')
2125 cs%id_div_xx_h = register_diag_field(
'ocean_model',
'div_xx_h', diag%axesTL, time, &
2126 'Horizontal divergence at h Points',
's-1')
2131 if (cs%use_GME)
then
2132 cs%id_GME_coeff_h = register_diag_field(
'ocean_model',
'GME_coeff_h', diag%axesTL, time, &
2133 'GME coefficient at h Points',
'm2 s-1', conversion=us%s_to_T)
2135 cs%id_GME_coeff_q = register_diag_field(
'ocean_model',
'GME_coeff_q', diag%axesBL, time, &
2136 'GME coefficient at q Points',
'm2 s-1', conversion=us%s_to_T)
2138 cs%id_FrictWork_GME = register_diag_field(
'ocean_model',
'FrictWork_GME',diag%axesTL,time,&
2139 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)',
'W m-2')
2142 cs%id_FrictWork = register_diag_field(
'ocean_model',
'FrictWork',diag%axesTL,time,&
2143 'Integral work done by lateral friction terms',
'W m-2')
2145 cs%id_FrictWork_diss = register_diag_field(
'ocean_model',
'FrictWork_diss',diag%axesTL,time,&
2146 'Integral work done by lateral friction terms (excluding diffusion of energy)',
'W m-2')
2148 if (
associated(meke))
then
2149 if (
associated(meke%mom_src))
then
2150 cs%id_FrictWorkMax = register_diag_field(
'ocean_model',
'FrictWorkMax', diag%axesTL, time,&
2151 'Maximum possible integral work done by lateral friction terms',
'W m-2')
2155 cs%id_FrictWorkIntz = register_diag_field(
'ocean_model',
'FrictWorkIntz',diag%axesT1,time, &
2156 'Depth integrated work done by lateral friction',
'W m-2', &
2157 cmor_field_name=
'dispkexyfo', &
2158 cmor_long_name=
'Depth integrated ocean kinetic energy dissipation due to lateral friction',&
2159 cmor_standard_name=
'ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
2161 if (cs%Laplacian .or. get_all)
then
2164 end subroutine hor_visc_init
2168 subroutine align_aniso_tensor_to_grid(CS, n1, n2)
2170 real,
intent(in) :: n1
2171 real,
intent(in) :: n2
2173 real :: recip_n2_norm
2176 recip_n2_norm = n1**2 + n2**2
2177 if (recip_n2_norm > 0.) recip_n2_norm = 1./recip_n2_norm
2179 cs%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
2180 cs%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
2181 cs%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
2182 cs%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
2184 end subroutine align_aniso_tensor_to_grid
2188 subroutine smooth_gme(CS,G,GME_flux_h,GME_flux_q)
2192 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(inout) :: GME_flux_h
2194 real,
dimension(SZIB_(G),SZJB_(G)),
optional,
intent(inout) :: GME_flux_q
2198 real,
dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original
2199 real,
dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original
2200 real :: wc, ww, we, wn, ws
2201 integer :: i, j, k, s
2207 if (
present(gme_flux_h))
then
2208 call pass_var(gme_flux_h, g%Domain)
2209 gme_flux_h_original = gme_flux_h
2214 if (g%mask2dT(i,j)==0.) cycle
2217 ww = 0.125 * g%mask2dT(i-1,j)
2218 we = 0.125 * g%mask2dT(i+1,j)
2219 ws = 0.125 * g%mask2dT(i,j-1)
2220 wn = 0.125 * g%mask2dT(i,j+1)
2221 wc = 1.0 - (ww+we+wn+ws)
2223 gme_flux_h(i,j) = wc * gme_flux_h_original(i,j) &
2224 + ww * gme_flux_h_original(i-1,j) &
2225 + we * gme_flux_h_original(i+1,j) &
2226 + ws * gme_flux_h_original(i,j-1) &
2227 + wn * gme_flux_h_original(i,j+1)
2232 if (
present(gme_flux_q))
then
2233 call pass_var(gme_flux_q, g%Domain, position=corner, complete=.true.)
2234 gme_flux_q_original = gme_flux_q
2236 do j = g%JscB, g%JecB
2237 do i = g%IscB, g%IecB
2239 if (g%mask2dBu(i,j)==0.) cycle
2242 ww = 0.125 * g%mask2dBu(i-1,j)
2243 we = 0.125 * g%mask2dBu(i+1,j)
2244 ws = 0.125 * g%mask2dBu(i,j-1)
2245 wn = 0.125 * g%mask2dBu(i,j+1)
2246 wc = 1.0 - (ww+we+wn+ws)
2248 gme_flux_q(i,j) = wc * gme_flux_q_original(i,j) &
2249 + ww * gme_flux_q_original(i-1,j) &
2250 + we * gme_flux_q_original(i+1,j) &
2251 + ws * gme_flux_q_original(i,j-1) &
2252 + wn * gme_flux_q_original(i,j+1)
2258 end subroutine smooth_gme
2261 subroutine hor_visc_end(CS)
2265 if (cs%Laplacian .or. cs%biharmonic)
then
2266 dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)
2267 dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)
2268 dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)
2271 if (cs%Laplacian)
then
2272 dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)
2273 if (cs%bound_Kh)
then
2274 dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
2276 if (cs%Smagorinsky_Kh)
then
2277 dealloc_(cs%Laplac2_const_xx) ; dealloc_(cs%Laplac2_const_xy)
2279 if (cs%Leith_Kh)
then
2280 dealloc_(cs%Laplac3_const_xx) ; dealloc_(cs%Laplac3_const_xy)
2284 if (cs%biharmonic)
then
2285 dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)
2286 dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)
2287 dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)
2288 if (cs%bound_Ah)
then
2289 dealloc_(cs%Ah_Max_xx) ; dealloc_(cs%Ah_Max_xy)
2291 if (cs%Smagorinsky_Ah)
then
2292 dealloc_(cs%Biharm5_const_xx) ; dealloc_(cs%Biharm5_const_xy)
2293 if (cs%bound_Coriolis)
then
2294 dealloc_(cs%Biharm5_const2_xx) ; dealloc_(cs%Biharm5_const2_xy)
2297 if (cs%Leith_Ah)
then
2298 dealloc_(cs%Biharm_const_xx) ; dealloc_(cs%Biharm_const_xy)
2301 if (cs%anisotropic)
then
2304 dealloc_(cs%n1n1_m_n2n2_h)
2305 dealloc_(cs%n1n1_m_n2n2_q)
2309 end subroutine hor_visc_end