23 implicit none ;
private
25 #include <MOM_memory.h>
43 logical :: better_bound_kh
47 logical :: better_bound_ah
53 logical :: smagorinsky_kh
55 logical :: smagorinsky_ah
59 logical :: modified_leith
61 logical :: use_beta_in_leith
64 logical :: use_qg_leith_visc
66 logical :: bound_coriolis
69 logical :: use_kh_bg_2d
72 logical :: use_land_mask
77 logical :: anisotropic
78 logical :: add_les_viscosity
81 logical :: dynamic_aniso
83 logical :: res_scale_meke
86 logical :: answers_2018
91 real :: gme_efficiency
94 real allocable_,
dimension(NIMEM_,NJMEM_) :: kh_bg_xx
98 real allocable_,
dimension(NIMEM_,NJMEM_) :: kh_bg_2d
102 real allocable_,
dimension(NIMEM_,NJMEM_) :: ah_bg_xx
111 real allocable_,
dimension(NIMEM_,NJMEM_) :: reduction_xx
114 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
115 kh_max_xx, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
116 ah_max_xx, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
117 n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points
120 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: kh_bg_xy
124 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: ah_bg_xy
133 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: reduction_xy
136 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
137 kh_max_xy, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
138 ah_max_xy, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
139 n1n2_q, & !< Factor n1*n2 in the anisotropic direction tensor at q-points
142 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
143 dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2]
144 dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2]
145 dx_dyt, & !< Pre-calculated dx/dy at h points [nondim]
147 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
148 dx2q, & !< Pre-calculated dx^2 at q points [L2 ~> m2]
149 dy2q, & !< Pre-calculated dy^2 at q points [L2 ~> m2]
150 dx_dybu, & !< Pre-calculated dx/dy at q points [nondim]
152 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: &
153 idx2dycu, & !< 1/(dx^2 dy) at u points [L-3 ~> m-3]
155 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: &
156 idx2dycv, & !< 1/(dx^2 dy) at v points [L-3 ~> m-3]
161 real allocable_,
dimension(NIMEM_,NJMEM_) :: &
162 laplac2_const_xx, & !< Laplacian metric-dependent constants [L2 ~> m2]
163 biharm5_const_xx, & !< Biharmonic metric-dependent constants [L5 ~> m5]
164 laplac3_const_xx, & !< Laplacian metric-dependent constants [L3 ~> m3]
165 biharm_const_xx, & !< Biharmonic metric-dependent constants [L4 ~> m4]
168 real allocable_,
dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
169 laplac2_const_xy, & !< Laplacian metric-dependent constants [L2 ~> m2]
170 biharm5_const_xy, & !< Biharmonic metric-dependent constants [L5 ~> m5]
171 laplac3_const_xy, & !< Laplacian metric-dependent constants [L3 ~> m3]
172 biharm_const_xy, & !< Biharmonic metric-dependent constants [L4 ~> m4]
179 integer :: id_diffu = -1, id_diffv = -1
180 integer :: id_ah_h = -1, id_ah_q = -1
181 integer :: id_kh_h = -1, id_kh_q = -1
182 integer :: id_gme_coeff_h = -1, id_gme_coeff_q = -1
183 integer :: id_vort_xy_q = -1, id_div_xx_h = -1
184 integer :: id_frictwork = -1, id_frictworkintz = -1
185 integer :: id_frictwork_gme = -1
205 subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
209 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
211 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
213 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
215 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
218 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
233 real,
dimension(SZIB_(G),SZJ_(G)) :: &
239 real,
dimension(SZI_(G),SZJB_(G)) :: &
245 real,
dimension(SZI_(G),SZJ_(G)) :: &
257 grad_vort_mag_h_2d, &
265 real,
dimension(SZIB_(G),SZJB_(G)) :: &
267 ddel2vdx, ddel2udy, &
278 grad_vort_mag_q_2d, &
286 real,
dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
290 gme_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]
293 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: &
295 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: &
297 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
304 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
314 real :: vert_vort_mag
324 real :: visc_bound_rem
333 real :: gme_coeff_limiter
338 real :: backscat_subround
341 logical :: rescale_kh, legacy_bound
342 logical :: find_frictwork
343 logical :: apply_obc = .false.
344 logical :: use_meke_ku
345 logical :: use_meke_au
346 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
347 integer :: i, j, k, n
348 real :: inv_pi3, inv_pi2, inv_pi5
349 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
350 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
352 h_neglect = gv%H_subroundoff
353 h_neglect3 = h_neglect**3
354 inv_pi3 = 1.0/((4.0*atan(1.0))**3)
355 inv_pi2 = 1.0/((4.0*atan(1.0))**2)
356 inv_pi5 = inv_pi3 * inv_pi2
361 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then
362 apply_obc = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally
364 endif ;
endif ;
endif
366 if (.not.
associated(cs))
call mom_error(fatal, &
367 "MOM_hor_visc: Module must be initialized before it is used.")
368 if (.not.(cs%Laplacian .or. cs%biharmonic))
return
370 find_frictwork = (cs%id_FrictWork > 0)
371 if (cs%id_FrictWorkIntz > 0) find_frictwork = .true.
372 if (
associated(meke))
then
373 if (
associated(meke%mom_src)) find_frictwork = .true.
374 backscat_subround = 0.0
375 if (find_frictwork .and.
associated(meke%mom_src) .and. (meke%backscatter_Ro_c > 0.0) .and. &
376 (meke%backscatter_Ro_Pow /= 0.0)) &
377 backscat_subround = (1.0e-16/meke%backscatter_Ro_c)**(1.0/meke%backscatter_Ro_Pow)
381 if (
associated(varmix))
then
382 rescale_kh = varmix%Resoln_scaled_Kh
383 if (rescale_kh .and. &
384 (.not.
associated(varmix%Res_fn_h) .or. .not.
associated(varmix%Res_fn_q))) &
385 call mom_error(fatal,
"MOM_hor_visc: VarMix%Res_fn_h and " //&
386 "VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
388 legacy_bound = (cs%Smagorinsky_Kh .or. cs%Leith_Kh) .and. &
389 (cs%bound_Kh .and. .not.cs%better_bound_Kh)
392 use_meke_ku =
associated(meke%Ku)
393 use_meke_au =
associated(meke%Au)
396 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
397 boundary_mask_h(i,j) = (g%mask2dCu(i,j) * g%mask2dCv(i,j) * g%mask2dCu(i-1,j) * g%mask2dCv(i,j-1))
400 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
401 boundary_mask_q(i,j) = (g%mask2dCv(i,j) * g%mask2dCv(i+1,j) * g%mask2dCu(i,j) * g%mask2dCu(i,j-1))
405 gme_coeff_h(:,:,:) = 0.0
406 gme_coeff_q(:,:,:) = 0.0
407 str_xx_gme(:,:) = 0.0
408 str_xy_gme(:,:) = 0.0
414 do j=js-1,je+1 ;
do i=is-1,ie+1
415 dudx_bt(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * ubtav(i,j) - &
416 g%IdyCu(i-1,j) * ubtav(i-1,j))
417 dvdy_bt(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * vbtav(i,j) - &
418 g%IdxCv(i,j-1) * vbtav(i,j-1))
421 call pass_vector(dudx_bt, dvdy_bt, g%Domain, stagger=bgrid_ne)
423 do j=jsq,jeq+1 ;
do i=isq,ieq+1
424 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
428 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
429 dvdx_bt(i,j) = cs%DY_dxBu(i,j)*(vbtav(i+1,j)*g%IdyCv(i+1,j) &
430 - vbtav(i,j)*g%IdyCv(i,j))
431 dudy_bt(i,j) = cs%DX_dyBu(i,j)*(ubtav(i,j+1)*g%IdxCu(i,j+1) &
432 - ubtav(i,j)*g%IdxCu(i,j))
435 call pass_vector(dvdx_bt, dudy_bt, g%Domain, stagger=agrid)
438 do j=js-1,jeq ;
do i=is-1,ieq
439 sh_xy_bt(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
442 do j=js-1,jeq ;
do i=is-1,ieq
443 sh_xy_bt(i,j) = g%mask2dBu(i,j) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
447 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
448 grad_vel_mag_bt_h(i,j) = boundary_mask_h(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
449 (0.25*((dvdx_bt(i,j)+dvdx_bt(i-1,j-1))+(dvdx_bt(i,j-1)+dvdx_bt(i-1,j))))**2 + &
450 (0.25*((dudy_bt(i,j)+dudy_bt(i-1,j-1))+(dudy_bt(i,j-1)+dudy_bt(i-1,j))))**2)
453 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
454 grad_vel_mag_bt_q(i,j) = boundary_mask_q(i,j) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
455 (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + &
456 (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)
495 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
496 dudx(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * u(i,j,k) - &
497 g%IdyCu(i-1,j) * u(i-1,j,k))
498 dvdy(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * v(i,j,k) - &
499 g%IdxCv(i,j-1) * v(i,j-1,k))
500 sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
504 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
505 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))
506 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))
514 if (cs%use_land_mask)
then
515 do j=js-2,je+2 ;
do i=isq-1,ieq+1
516 h_u(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k))
518 do j=jsq-1,jeq+1 ;
do i=is-2,ie+2
519 h_v(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k))
522 do j=js-2,je+2 ;
do i=isq-1,ieq+1
523 h_u(i,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
525 do j=jsq-1,jeq+1 ;
do i=is-2,ie+2
526 h_v(i,j) = 0.5 * (h(i,j,k) + h(i,j+1,k))
532 if (apply_obc)
then ;
do n=1,obc%number_of_segments
533 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
534 if (obc%zero_strain .or. obc%freeslip_strain .or. obc%computed_strain)
then
535 if (obc%segment(n)%is_N_or_S .and. (j >= js-2) .and. (j <= jeq+1))
then
536 do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
537 if (obc%zero_strain)
then
538 dvdx(i,j) = 0. ; dudy(i,j) = 0.
539 elseif (obc%freeslip_strain)
then
541 elseif (obc%computed_strain)
then
543 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &
544 (obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%IdxCu(i,j)
546 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &
547 (u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdxCu(i,j+1)
549 elseif (obc%specified_strain)
then
551 dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j)*g%dxBu(i,j)
553 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)
557 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-2) .and. (i <= ieq+1))
then
558 do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
559 if (obc%zero_strain)
then
560 dvdx(i,j) = 0. ; dudy(i,j) = 0.
561 elseif (obc%freeslip_strain)
then
563 elseif (obc%computed_strain)
then
564 if (obc%segment(n)%direction == obc_direction_e)
then
565 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &
566 (obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%IdyCv(i,j)
568 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &
569 (v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdyCv(i+1,j)
571 elseif (obc%specified_strain)
then
572 if (obc%segment(n)%direction == obc_direction_e)
then
573 dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i,j)*g%dxBu(i,j)
575 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)
588 if ((j >= jsq-1) .and. (j <= jeq+1))
then
589 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
594 if ((j >= jsq-1) .and. (j <= jeq+1))
then
595 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
596 h_v(i,j) = h(i,j+1,k)
599 elseif (obc%segment(n)%direction == obc_direction_e)
then
600 if ((i >= isq-1) .and. (i <= ieq+1))
then
601 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
605 elseif (obc%segment(n)%direction == obc_direction_w)
then
606 if ((i >= isq-1) .and. (i <= ieq+1))
then
607 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
608 h_u(i,j) = h(i+1,j,k)
614 if (apply_obc)
then ;
do n=1,obc%number_of_segments
615 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
617 if ((j >= js-2) .and. (j <= je))
then
618 do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
619 h_u(i,j+1) = h_u(i,j)
623 if ((j >= js-1) .and. (j <= je+1))
then
624 do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+1,obc%segment(n)%HI%ied)
625 h_u(i,j) = h_u(i,j+1)
628 elseif (obc%segment(n)%direction == obc_direction_e)
then
629 if ((i >= is-2) .and. (i <= ie))
then
630 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
631 h_v(i+1,j) = h_v(i,j)
634 elseif (obc%segment(n)%direction == obc_direction_w)
then
635 if ((i >= is-1) .and. (i <= ie+1))
then
636 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
637 h_v(i,j) = h_v(i+1,j)
646 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
647 sh_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) + dudy(i,j) )
650 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
651 sh_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) + dudy(i,j) )
656 if (cs%biharmonic)
then
657 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
658 del2u(i,j) = cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*sh_xx(i+1,j) - cs%dy2h(i,j)*sh_xx(i,j)) + &
659 cs%Idx2dyCu(i,j)*(cs%dx2q(i,j)*sh_xy(i,j) - cs%dx2q(i,j-1)*sh_xy(i,j-1))
661 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
662 del2v(i,j) = cs%Idxdy2v(i,j)*(cs%dy2q(i,j)*sh_xy(i,j) - cs%dy2q(i-1,j)*sh_xy(i-1,j)) - &
663 cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*sh_xx(i,j+1) - cs%dx2h(i,j)*sh_xx(i,j))
665 if (apply_obc) then;
if (obc%zero_biharmonic)
then
666 do n=1,obc%number_of_segments
667 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
668 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1))
then
669 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
672 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1))
then
673 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
681 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
685 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
686 vort_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) - dudy(i,j) )
689 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
690 vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )
695 do j=js-2,jeq+1 ;
do i=is-1,ieq+1
696 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)
697 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))
700 do j=js-1,jeq+1 ;
do i=is-2,ieq+1
701 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)
702 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))
705 if (cs%modified_Leith)
then
707 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
708 div_xx(i,j) = dudx(i,j) + dvdy(i,j)
712 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+1
713 div_xx_dx(i,j) = g%IdxCu(i,j)*(div_xx(i+1,j) - div_xx(i,j))
715 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+2
716 div_xx_dy(i,j) = g%IdyCv(i,j)*(div_xx(i,j+1) - div_xx(i,j))
720 do j=jsq,jeq+1 ;
do i=isq,ieq+1
721 grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i-1,j)))**2 + &
722 (0.5 * (div_xx_dy(i,j) + div_xx_dy(i,j-1)))**2)
725 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+1
726 grad_div_mag_q(i,j) =sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i,j+1)))**2 + &
727 (0.5 * (div_xx_dy(i,j) + div_xx_dy(i+1,j)))**2)
732 do j=jsq-1,jeq+2 ;
do i=is-2,ieq+1
735 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+2
738 do j=jsq,jeq+1 ;
do i=isq,ieq+1
739 grad_div_mag_h(i,j) = 0.0
741 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+1
742 grad_div_mag_q(i,j) = 0.0
748 if (cs%use_beta_in_Leith)
then
749 do j=js-2,jeq+1 ;
do i=is-1,ieq+1
750 vort_xy_dx(i,j) = vort_xy_dx(i,j) + 0.5 * ( g%dF_dx(i,j) + g%dF_dx(i,j+1))
752 do j=js-1,jeq+1 ;
do i=is-2,ieq+1
753 vort_xy_dy(i,j) = vort_xy_dy(i,j) + 0.5 * ( g%dF_dy(i,j) + g%dF_dy(i+1,j))
757 if (cs%use_QG_Leith_visc)
then
759 do j=jsq,jeq+1 ;
do i=isq,ieq+1
760 grad_vort_mag_h_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
761 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
763 do j=js-1,jeq ;
do i=is-1,ieq
764 grad_vort_mag_q_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
765 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
770 vort_xy_dx, vort_xy_dy)
774 do j=jsq,jeq+1 ;
do i=isq,ieq+1
775 grad_vort_mag_h(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
776 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
778 do j=js-1,jeq ;
do i=is-1,ieq
779 grad_vort_mag_q(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
780 (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
787 do j=jsq,jeq+1 ;
do i=isq,ieq+1
788 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah))
then
789 shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
790 0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
791 (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
793 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
794 if (cs%use_QG_Leith_visc)
then
795 vert_vort_mag = min(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3.*grad_vort_mag_h_2d(i,j))
797 vert_vort_mag = (grad_vort_mag_h(i,j) + grad_div_mag_h(i,j))
800 if (cs%better_bound_Ah .or. cs%better_bound_Kh)
then
801 hrat_min = min(1.0, min(h_u(i,j), h_u(i-1,j), h_v(i,j), h_v(i,j-1)) / &
802 (h(i,j,k) + h_neglect) )
806 if (cs%Laplacian)
then
809 kh = cs%Kh_bg_xx(i,j)
810 if (cs%add_LES_viscosity)
then
811 if (cs%Smagorinsky_Kh) kh = kh + cs%Laplac2_const_xx(i,j) * shear_mag
812 if (cs%Leith_Kh) kh = kh + cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3
814 if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xx(i,j) * shear_mag )
815 if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3)
818 if (rescale_kh) kh = varmix%Res_fn_h(i,j) * kh
819 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_h(i,j)
821 if (legacy_bound) kh = min(kh, cs%Kh_Max_xx(i,j))
822 kh = max( kh, cs%Kh_bg_min )
824 kh = kh + meke%Ku(i,j) * meke_res_fn
825 if (cs%anisotropic) kh = kh + cs%Kh_aniso * ( 1. - cs%n1n2_h(i,j)**2 )
829 if (cs%better_bound_Kh)
then
830 if (kh >= hrat_min*cs%Kh_Max_xx(i,j))
then
832 kh = hrat_min*cs%Kh_Max_xx(i,j)
834 visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xx(i,j))
838 if ((cs%id_Kh_h>0) .or. find_frictwork .or. cs%debug) kh_h(i,j,k) = kh
839 if (cs%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)
841 str_xx(i,j) = -kh * sh_xx(i,j)
846 if (cs%anisotropic)
then
848 local_strain = 0.25 * ( (sh_xy(i,j) + sh_xy(i-1,j-1)) + (sh_xy(i-1,j) + sh_xy(i,j-1)) )
850 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
853 if (cs%biharmonic)
then
856 ahsm = 0.0; ahlth = 0.0
857 if ((cs%Smagorinsky_Ah) .or. (cs%Leith_Ah))
then
858 if (cs%Smagorinsky_Ah)
then
859 if (cs%bound_Coriolis)
then
860 ahsm = shear_mag * (cs%Biharm_const_xx(i,j) + &
861 cs%Biharm_const2_xx(i,j)*shear_mag)
863 ahsm = cs%Biharm_const_xx(i,j) * shear_mag
866 if (cs%Leith_Ah) ahlth = cs%Biharm5_const_xx(i,j) * vert_vort_mag * inv_pi5
867 ah = max(max(cs%Ah_bg_xx(i,j), ahsm), ahlth)
868 if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
869 ah = min(ah, cs%Ah_Max_xx(i,j))
871 ah = cs%Ah_bg_xx(i,j)
874 if (use_meke_au) ah = ah + meke%Au(i,j)
876 if (cs%better_bound_Ah)
then
877 ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xx(i,j))
880 if ((cs%id_Ah_h>0) .or. find_frictwork .or. cs%debug) ah_h(i,j,k) = ah
882 str_xx(i,j) = str_xx(i,j) + ah * &
883 (cs%DY_dxT(i,j) * (g%IdyCu(i,j)*del2u(i,j) - g%IdyCu(i-1,j)*del2u(i-1,j)) - &
884 cs%DX_dyT(i,j) * (g%IdxCv(i,j)*del2v(i,j) - g%IdxCv(i,j-1)*del2v(i,j-1)))
887 bhstr_xx(i,j) = ah * &
888 (cs%DY_dxT(i,j) * (g%IdyCu(i,j)*del2u(i,j) - g%IdyCu(i-1,j)*del2u(i-1,j)) - &
889 cs%DX_dyT(i,j) * (g%IdxCv(i,j)*del2v(i,j) - g%IdxCv(i,j-1)*del2v(i,j-1)))
890 bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
896 if (cs%biharmonic)
then
898 do j=js-1,jeq ;
do i=is-1,ieq
899 ddel2vdx(i,j) = cs%DY_dxBu(i,j)*(del2v(i+1,j)*g%IdyCv(i+1,j) - del2v(i,j)*g%IdyCv(i,j))
900 ddel2udy(i,j) = cs%DX_dyBu(i,j)*(del2u(i,j+1)*g%IdxCu(i,j+1) - del2u(i,j)*g%IdxCu(i,j))
903 if (apply_obc)
then ;
if (obc%zero_strain .or. obc%freeslip_strain)
then
904 do n=1,obc%number_of_segments
905 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
906 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= jeq))
then
907 do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
908 if (obc%zero_strain)
then
909 ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.
910 elseif (obc%freeslip_strain)
then
914 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ieq))
then
915 do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
916 if (obc%zero_strain)
then
917 ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.
918 elseif (obc%freeslip_strain)
then
929 do j=js-1,jeq ;
do i=is-1,ieq
930 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah))
then
931 shear_mag = sqrt(sh_xy(i,j)*sh_xy(i,j) + &
932 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + &
933 (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j))))
935 if ((cs%Leith_Kh) .or. (cs%Leith_Ah))
then
936 if (cs%use_QG_Leith_visc)
then
937 vert_vort_mag = min(grad_vort_mag_q(i,j) + grad_div_mag_q(i,j), 3.*grad_vort_mag_q_2d(i,j))
939 vert_vort_mag = (grad_vort_mag_q(i,j) + grad_div_mag_q(i,j))
942 h2uq = 4.0 * h_u(i,j) * h_u(i,j+1)
943 h2vq = 4.0 * h_v(i,j) * h_v(i+1,j)
944 hq(i,j) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
945 ((h_u(i,j) + h_u(i,j+1)) + (h_v(i,j) + h_v(i+1,j))))
947 if (cs%better_bound_Ah .or. cs%better_bound_Kh)
then
948 hrat_min = min(1.0, min(h_u(i,j), h_u(i,j+1), h_v(i,j), h_v(i+1,j)) / &
949 (hq(i,j) + h_neglect) )
953 if (cs%no_slip .and. (g%mask2dBu(i,j) < 0.5))
then
954 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
955 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0)
then
958 hu = g%mask2dCu(i,j) * h_u(i,j) + g%mask2dCu(i,j+1) * h_u(i,j+1)
959 hv = g%mask2dCv(i,j) * h_v(i,j) + g%mask2dCv(i+1,j) * h_v(i+1,j)
960 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) * &
961 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) == 0.0)
then
967 hq(i,j) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
968 hrat_min = min(1.0, min(hu, hv) / (hq(i,j) + h_neglect) )
973 if (cs%Laplacian)
then
976 kh = cs%Kh_bg_xy(i,j)
977 if (cs%add_LES_viscosity)
then
978 if (cs%Smagorinsky_Kh) kh = kh + cs%Laplac2_const_xx(i,j) * shear_mag
979 if (cs%Leith_Kh) kh = kh + cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3
981 if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xy(i,j) * shear_mag )
982 if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xy(i,j) * vert_vort_mag*inv_pi3)
985 if (rescale_kh) kh = varmix%Res_fn_q(i,j) * kh
986 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_q(i,j)
988 if (legacy_bound) kh = min(kh, cs%Kh_Max_xy(i,j))
989 kh = max( kh, cs%Kh_bg_min )
990 if (use_meke_ku)
then
991 kh = kh + 0.25*( (meke%Ku(i,j) + meke%Ku(i+1,j+1)) + &
992 (meke%Ku(i+1,j) + meke%Ku(i,j+1)) ) * meke_res_fn
995 if (cs%anisotropic) kh = kh + cs%Kh_aniso * cs%n1n2_q(i,j)**2
999 if (cs%better_bound_Kh)
then
1000 if (kh >= hrat_min*cs%Kh_Max_xy(i,j))
then
1001 visc_bound_rem = 0.0
1002 kh = hrat_min*cs%Kh_Max_xy(i,j)
1003 elseif (cs%Kh_Max_xy(i,j)>0.)
then
1004 visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xy(i,j))
1008 if (cs%id_Kh_q>0 .or. cs%debug) kh_q(i,j,k) = kh
1009 if (cs%id_vort_xy_q>0) vort_xy_q(i,j,k) = vort_xy(i,j)
1011 str_xy(i,j) = -kh * sh_xy(i,j)
1016 if (cs%anisotropic)
then
1018 local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )
1020 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
1023 if (cs%biharmonic)
then
1026 ahsm = 0.0 ; ahlth = 0.0
1027 if (cs%Smagorinsky_Ah .or. cs%Leith_Ah)
then
1028 if (cs%Smagorinsky_Ah)
then
1029 if (cs%bound_Coriolis)
then
1030 ahsm = shear_mag * (cs%Biharm_const_xy(i,j) + &
1031 cs%Biharm_const2_xy(i,j)*shear_mag)
1033 ahsm = cs%Biharm_const_xy(i,j) * shear_mag
1036 if (cs%Leith_Ah) ahlth = cs%Biharm5_const_xy(i,j) * vert_vort_mag * inv_pi5
1037 ah = max(max(cs%Ah_bg_xy(i,j), ahsm), ahlth)
1038 if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
1039 ah = min(ah, cs%Ah_Max_xy(i,j))
1041 ah = cs%Ah_bg_xy(i,j)
1044 if (use_meke_au)
then
1045 ah = ah + 0.25*( (meke%Au(i,j) + meke%Au(i+1,j+1)) + &
1046 (meke%Au(i+1,j) + meke%Au(i,j+1)) )
1049 if (cs%better_bound_Ah)
then
1050 ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xy(i,j))
1053 if (cs%id_Ah_q>0 .or. cs%debug) ah_q(i,j,k) = ah
1055 str_xy(i,j) = str_xy(i,j) + ah * ( ddel2vdx(i,j) + ddel2udy(i,j) )
1058 bhstr_xy(i,j) = ah * ( ddel2vdx(i,j) + ddel2udy(i,j) ) * &
1059 (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
1065 if (cs%use_GME)
then
1066 if (cs%answers_2018)
then
1067 do j=js,je ;
do i=is,ie
1068 grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
1069 (0.25*((dvdx(i,j) + dvdx(i-1,j-1)) + (dvdx(i,j-1) + dvdx(i-1,j))) )**2 + &
1070 (0.25*((dudy(i,j) + dudy(i-1,j-1)) + (dudy(i,j-1) + dudy(i-1,j))) )**2)
1071 max_diss_rate_h(i,j,k) = 2.0 * meke%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
1074 do j=js,je ;
do i=is,ie
1075 grad_vel_mag_h(i,j) = boundary_mask_h(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + &
1076 ((0.25*((dvdx(i,j) + dvdx(i-1,j-1)) + (dvdx(i,j-1) + dvdx(i-1,j))) )**2 + &
1077 (0.25*((dudy(i,j) + dudy(i-1,j-1)) + (dudy(i,j-1) + dudy(i-1,j))) )**2))
1078 max_diss_rate_h(i,j,k) = 2.0 * meke%MEKE(i,j) * sqrt(grad_vel_mag_h(i,j))
1082 if (cs%answers_2018)
then
1083 do j = g%JscB, g%JecB ;
do i = g%IscB, g%IecB
1084 grad_vel_mag_q(i,j) = boundary_mask_q(i,j) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
1085 (0.25*((dvdx(i,j)+dvdx(i-1,j-1)) + (dvdx(i,j-1)+dvdx(i-1,j))) )**2 + &
1086 (0.25*((dudy(i,j)+dudy(i-1,j-1)) + (dudy(i,j-1)+dudy(i-1,j))) )**2)
1088 max_diss_rate_q(i,j,k) = 0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)+ &
1089 meke%MEKE(i,j+1)+meke%MEKE(i+1,j+1)) * sqrt(grad_vel_mag_q(i,j))
1092 do j = g%JscB, g%JecB ;
do i = g%IscB, g%IecB
1093 grad_vel_mag_q(i,j) = boundary_mask_q(i,j) * ((dudx(i,j)**2 + dvdy(i,j)**2) + &
1094 ((0.25*((dvdx(i,j)+dvdx(i-1,j-1)) + (dvdx(i,j-1)+dvdx(i-1,j))) )**2 + &
1095 (0.25*((dudy(i,j)+dudy(i-1,j-1)) + (dudy(i,j-1)+dudy(i-1,j))) )**2))
1097 max_diss_rate_q(i,j,k) = 0.5*((meke%MEKE(i,j) + meke%MEKE(i+1,j+1)) + &
1098 (meke%MEKE(i+1,j) + meke%MEKE(i,j+1))) * sqrt(grad_vel_mag_q(i,j))
1102 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1103 if ((grad_vel_mag_bt_h(i,j)>0) .and. (max_diss_rate_h(i,j,k)>0))
then
1104 gme_coeff = (min(g%bathyT(i,j)/cs%GME_h0,1.0)**2) * cs%GME_efficiency*max_diss_rate_h(i,j,k) / &
1105 grad_vel_mag_bt_h(i,j)
1111 gme_coeff = gme_coeff * boundary_mask_h(i,j)
1112 gme_coeff = min(gme_coeff, cs%GME_limiter)
1114 if ((cs%id_GME_coeff_h>0) .or. find_frictwork) gme_coeff_h(i,j,k) = gme_coeff
1115 str_xx_gme(i,j) = gme_coeff * sh_xx_bt(i,j)
1119 do j=js-1,jeq ;
do i=is-1,ieq
1121 if ((grad_vel_mag_bt_q(i,j)>0) .and. (max_diss_rate_q(i,j,k)>0))
then
1122 gme_coeff = (min(g%bathyT(i,j)/cs%GME_h0,1.0)**2) * cs%GME_efficiency*max_diss_rate_q(i,j,k) / &
1123 grad_vel_mag_bt_q(i,j)
1129 gme_coeff = gme_coeff * boundary_mask_q(i,j)
1130 gme_coeff = min(gme_coeff, cs%GME_limiter)
1132 if (cs%id_GME_coeff_q>0) gme_coeff_q(i,j,k) = gme_coeff
1133 str_xy_gme(i,j) = gme_coeff * sh_xy_bt(i,j)
1141 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1142 str_xx(i,j) = (str_xx(i,j) + str_xx_gme(i,j)) * (h(i,j,k) * cs%reduction_xx(i,j))
1145 do j=js-1,jeq ;
do i=is-1,ieq
1147 if (cs%no_slip)
then
1148 str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * cs%reduction_xy(i,j))
1150 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))
1154 if (
associated(meke%GME_snk))
then
1155 do j=js,je ;
do i=is,ie
1156 frictwork_gme(i,j,k) = gme_coeff_h(i,j,k) * h(i,j,k) * gv%H_to_RZ * grad_vel_mag_bt_h(i,j)
1161 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1162 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
1165 do j=js-1,jeq ;
do i=is-1,ieq
1166 if (cs%no_slip)
then
1167 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * cs%reduction_xy(i,j))
1169 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
1176 do j=js,je ;
do i=isq,ieq
1177 diffu(i,j,k) = ((g%IdyCu(i,j)*(cs%dy2h(i,j) *str_xx(i,j) - &
1178 cs%dy2h(i+1,j)*str_xx(i+1,j)) + &
1179 g%IdxCu(i,j)*(cs%dx2q(i,j-1)*str_xy(i,j-1) - &
1180 cs%dx2q(i,j) *str_xy(i,j))) * &
1181 g%IareaCu(i,j)) / (h_u(i,j) + h_neglect)
1187 do n=1,obc%number_of_segments
1188 if (obc%segment(n)%is_E_or_W)
then
1189 i = obc%segment(n)%HI%IsdB
1190 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
1198 do j=jsq,jeq ;
do i=is,ie
1199 diffv(i,j,k) = ((g%IdyCv(i,j)*(cs%dy2q(i-1,j)*str_xy(i-1,j) - &
1200 cs%dy2q(i,j) *str_xy(i,j)) - &
1201 g%IdxCv(i,j)*(cs%dx2h(i,j) *str_xx(i,j) - &
1202 cs%dx2h(i,j+1)*str_xx(i,j+1))) * &
1203 g%IareaCv(i,j)) / (h_v(i,j) + h_neglect)
1208 do n=1,obc%number_of_segments
1209 if (obc%segment(n)%is_N_or_S)
then
1210 j = obc%segment(n)%HI%JsdB
1211 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
1218 if (find_frictwork)
then ;
do j=js,je ;
do i=is,ie
1221 frictwork(i,j,k) = gv%H_to_RZ * ( &
1222 (str_xx(i,j)*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
1223 -str_xx(i,j)*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
1224 +0.25*((str_xy(i,j)*( &
1225 (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
1226 +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
1227 +str_xy(i-1,j-1)*( &
1228 (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
1229 +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
1231 (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
1232 +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
1234 (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
1235 +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
1236 enddo ;
enddo ;
endif
1241 if (find_frictwork .and.
associated(meke))
then ;
if (
associated(meke%mom_src))
then
1243 do j=js,je ;
do i=is,ie
1244 meke%mom_src(i,j) = 0.
1246 if (
associated(meke%GME_snk))
then
1247 do j=js,je ;
do i=is,ie
1248 meke%GME_snk(i,j) = 0.
1252 if (meke%backscatter_Ro_c /= 0.)
then
1253 do j=js,je ;
do i=is,ie
1254 fath = 0.25*( (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1255 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))) )
1256 shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
1257 0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
1258 (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
1259 if (cs%answers_2018)
then
1260 fath = (us%s_to_T*fath)**meke%backscatter_Ro_pow
1263 shear_mag = ( ( (us%s_to_T*shear_mag)**meke%backscatter_Ro_pow ) + 1.e-30 ) &
1264 * meke%backscatter_Ro_c
1267 roscl = shear_mag / ( fath + shear_mag )
1269 if (fath <= backscat_subround*shear_mag)
then
1272 sh_f_pow = meke%backscatter_Ro_c * (shear_mag / fath)**meke%backscatter_Ro_pow
1273 roscl = sh_f_pow / (1.0 + sh_f_pow)
1278 meke%mom_src(i,j) = meke%mom_src(i,j) + gv%H_to_RZ * ( &
1279 ((str_xx(i,j)-roscl*bhstr_xx(i,j))*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j) &
1280 -(str_xx(i,j)-roscl*bhstr_xx(i,j))*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
1281 +0.25*(((str_xy(i,j)-roscl*bhstr_xy(i,j))*( &
1282 (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j) &
1283 +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) ) &
1284 +(str_xy(i-1,j-1)-roscl*bhstr_xy(i-1,j-1))*( &
1285 (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1) &
1286 +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
1287 +((str_xy(i-1,j)-roscl*bhstr_xy(i-1,j))*( &
1288 (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j) &
1289 +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) ) &
1290 +(str_xy(i,j-1)-roscl*bhstr_xy(i,j-1))*( &
1291 (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1) &
1292 +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
1296 do j=js,je ;
do i=is,ie
1297 meke%mom_src(i,j) = meke%mom_src(i,j) + frictwork(i,j,k)
1300 if (cs%use_GME .and.
associated(meke))
then ;
if (
associated(meke%GME_snk))
then
1301 do j=js,je ;
do i=is,ie
1302 meke%GME_snk(i,j) = meke%GME_snk(i,j) + frictwork_gme(i,j,k)
1311 if (cs%id_diffu>0)
call post_data(cs%id_diffu, diffu, cs%diag)
1312 if (cs%id_diffv>0)
call post_data(cs%id_diffv, diffv, cs%diag)
1313 if (cs%id_FrictWork>0)
call post_data(cs%id_FrictWork, frictwork, cs%diag)
1314 if (cs%id_FrictWork_GME>0)
call post_data(cs%id_FrictWork_GME, frictwork_gme, cs%diag)
1315 if (cs%id_Ah_h>0)
call post_data(cs%id_Ah_h, ah_h, cs%diag)
1316 if (cs%id_div_xx_h>0)
call post_data(cs%id_div_xx_h, div_xx_h, cs%diag)
1317 if (cs%id_vort_xy_q>0)
call post_data(cs%id_vort_xy_q, vort_xy_q, cs%diag)
1318 if (cs%id_Ah_q>0)
call post_data(cs%id_Ah_q, ah_q, cs%diag)
1319 if (cs%id_Kh_h>0)
call post_data(cs%id_Kh_h, kh_h, cs%diag)
1320 if (cs%id_Kh_q>0)
call post_data(cs%id_Kh_q, kh_q, cs%diag)
1321 if (cs%id_GME_coeff_h > 0)
call post_data(cs%id_GME_coeff_h, gme_coeff_h, cs%diag)
1322 if (cs%id_GME_coeff_q > 0)
call post_data(cs%id_GME_coeff_q, gme_coeff_q, cs%diag)
1325 if (cs%Laplacian)
then
1326 call hchksum(kh_h,
"Kh_h", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
1327 call bchksum(kh_q,
"Kh_q", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
1329 if (cs%biharmonic)
call hchksum(ah_h,
"Ah_h", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
1330 if (cs%biharmonic)
call bchksum(ah_q,
"Ah_q", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
1333 if (cs%id_FrictWorkIntz > 0)
then
1335 do i=is,ie ; frictworkintz(i,j) = frictwork(i,j,1) ;
enddo
1336 do k=2,nz ;
do i=is,ie
1337 frictworkintz(i,j) = frictworkintz(i,j) + frictwork(i,j,k)
1340 call post_data(cs%id_FrictWorkIntz, frictworkintz, cs%diag)
1348 subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
1349 type(time_type),
intent(in) :: time
1354 type(
diag_ctrl),
target,
intent(inout) :: diag
1358 real,
dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v
1359 real,
dimension(SZI_(G),SZJB_(G)) :: v0u, v0v
1370 real :: boundcorconst
1376 real :: kh_vel_scale
1377 real :: ah_vel_scale
1378 real :: ah_time_scale
1379 real :: smag_lap_const
1380 real :: smag_bi_const
1381 real :: leith_lap_const
1382 real :: leith_bi_const
1387 real :: bound_cor_vel
1391 real :: kh_pwr_of_sine
1392 logical :: bound_cor_def
1400 logical :: default_2018_answers
1402 character(len=64) :: inputdir, filename
1405 real :: aniso_grid_dir(2)
1406 integer :: aniso_mode
1407 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
1408 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
1412 #include "version_variable.h"
1413 character(len=40) :: mdl =
"MOM_hor_visc"
1415 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1416 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1417 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1418 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1420 if (
associated(cs))
then
1421 call mom_error(warning,
"hor_visc_init called with an associated "// &
1422 "control structure.")
1434 cs%bound_Kh = .false. ; cs%better_bound_Kh = .false. ; cs%Smagorinsky_Kh = .false. ; cs%Leith_Kh = .false.
1435 cs%bound_Ah = .false. ; cs%better_bound_Ah = .false. ; cs%Smagorinsky_Ah = .false. ; cs%Leith_Ah = .false.
1436 cs%use_QG_Leith_visc = .false.
1437 cs%bound_Coriolis = .false.
1438 cs%Modified_Leith = .false.
1439 cs%anisotropic = .false.
1440 cs%dynamic_aniso = .false.
1446 call get_param(param_file, mdl,
"GET_ALL_PARAMS", get_all, default=.false.)
1448 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1449 "This sets the default value for the various _2018_ANSWERS parameters.", &
1451 call get_param(param_file, mdl,
"HOR_VISC_2018_ANSWERS", cs%answers_2018, &
1452 "If true, use the order of arithmetic and expressions that recover the "//&
1453 "answers from the end of 2018. Otherwise, use updated and more robust "//&
1454 "forms of the same expressions.", default=default_2018_answers)
1456 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false.)
1458 call get_param(param_file, mdl,
"LAPLACIAN", cs%Laplacian, &
1459 "If true, use a Laplacian horizontal viscosity.", &
1461 if (cs%Laplacian .or. get_all)
then
1462 call get_param(param_file, mdl,
"KH", kh, &
1463 "The background Laplacian horizontal viscosity.", &
1464 units =
"m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1465 call get_param(param_file, mdl,
"KH_BG_MIN", cs%Kh_bg_min, &
1466 "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
1467 units =
"m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1468 call get_param(param_file, mdl,
"KH_VEL_SCALE", kh_vel_scale, &
1469 "The velocity scale which is multiplied by the grid "//&
1470 "spacing to calculate the Laplacian viscosity. "//&
1471 "The final viscosity is the largest of this scaled "//&
1472 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
1473 units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
1474 call get_param(param_file, mdl,
"KH_SIN_LAT", kh_sin_lat, &
1475 "The amplitude of a latitudinally-dependent background "//&
1476 "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &
1477 units =
"m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1478 if (kh_sin_lat>0. .or. get_all) &
1479 call get_param(param_file, mdl,
"KH_PWR_OF_SINE", kh_pwr_of_sine, &
1480 "The power used to raise SIN(LAT) when using a latitudinally "//&
1481 "dependent background viscosity.", &
1482 units =
"nondim", default=4.0)
1484 call get_param(param_file, mdl,
"SMAGORINSKY_KH", cs%Smagorinsky_Kh, &
1485 "If true, use a Smagorinsky nonlinear eddy viscosity.", &
1487 if (cs%Smagorinsky_Kh .or. get_all) &
1488 call get_param(param_file, mdl,
"SMAG_LAP_CONST", smag_lap_const, &
1489 "The nondimensional Laplacian Smagorinsky constant, "//&
1490 "often 0.15.", units=
"nondim", default=0.0, &
1491 fail_if_missing = cs%Smagorinsky_Kh)
1493 call get_param(param_file, mdl,
"LEITH_KH", cs%Leith_Kh, &
1494 "If true, use a Leith nonlinear eddy viscosity.", &
1497 call get_param(param_file, mdl,
"MODIFIED_LEITH", cs%Modified_Leith, &
1498 "If true, add a term to Leith viscosity which is "//&
1499 "proportional to the gradient of divergence.", &
1501 call get_param(param_file, mdl,
"RES_SCALE_MEKE_VISC", cs%res_scale_MEKE, &
1502 "If true, the viscosity contribution from MEKE is scaled by "//&
1503 "the resolution function.", default=.false.)
1505 if (cs%Leith_Kh .or. get_all)
then
1506 call get_param(param_file, mdl,
"LEITH_LAP_CONST", leith_lap_const, &
1507 "The nondimensional Laplacian Leith constant, "//&
1508 "often set to 1.0", units=
"nondim", default=0.0, &
1509 fail_if_missing = cs%Leith_Kh)
1510 call get_param(param_file, mdl,
"USE_QG_LEITH_VISC", cs%use_QG_Leith_visc, &
1511 "If true, use QG Leith nonlinear eddy viscosity.", &
1513 if (cs%use_QG_Leith_visc .and. .not. cs%Leith_Kh)
call mom_error(fatal, &
1514 "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1515 "LEITH_KH must be True when USE_QG_LEITH_VISC=True.")
1517 if (cs%Leith_Kh .or. cs%Leith_Ah .or. get_all)
then
1518 call get_param(param_file, mdl,
"USE_BETA_IN_LEITH", cs%use_beta_in_Leith, &
1519 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1520 default=cs%Leith_Kh)
1521 call get_param(param_file, mdl,
"MODIFIED_LEITH", cs%modified_Leith, &
1522 "If true, add a term to Leith viscosity which is \n"//&
1523 "proportional to the gradient of divergence.", &
1526 call get_param(param_file, mdl,
"BOUND_KH", cs%bound_Kh, &
1527 "If true, the Laplacian coefficient is locally limited "//&
1528 "to be stable.", default=.true.)
1529 call get_param(param_file, mdl,
"BETTER_BOUND_KH", cs%better_bound_Kh, &
1530 "If true, the Laplacian coefficient is locally limited "//&
1531 "to be stable with a better bounding than just BOUND_KH.", &
1532 default=cs%bound_Kh)
1533 call get_param(param_file, mdl,
"ANISOTROPIC_VISCOSITY", cs%anisotropic, &
1534 "If true, allow anistropic viscosity in the Laplacian "//&
1535 "horizontal viscosity.", default=.false.)
1536 call get_param(param_file, mdl,
"ADD_LES_VISCOSITY", cs%add_LES_viscosity, &
1537 "If true, adds the viscosity from Smagorinsky and Leith to the "//&
1538 "background viscosity instead of taking the maximum.", default=.false.)
1540 if (cs%anisotropic .or. get_all)
then
1541 call get_param(param_file, mdl,
"KH_ANISO", cs%Kh_aniso, &
1542 "The background Laplacian anisotropic horizontal viscosity.", &
1543 units =
"m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1544 call get_param(param_file, mdl,
"ANISOTROPIC_MODE", aniso_mode, &
1545 "Selects the mode for setting the direction of anistropy.\n"//&
1546 "\t 0 - Points along the grid i-direction.\n"//&
1547 "\t 1 - Points towards East.\n"//&
1548 "\t 2 - Points along the flow direction, U/|U|.", &
1550 select case (aniso_mode)
1552 call get_param(param_file, mdl,
"ANISO_GRID_DIR", aniso_grid_dir, &
1553 "The vector pointing in the direction of anistropy for "//&
1554 "horizont viscosity. n1,n2 are the i,j components relative "//&
1555 "to the grid.", units =
"nondim", fail_if_missing=.true.)
1557 call get_param(param_file, mdl,
"ANISO_GRID_DIR", aniso_grid_dir, &
1558 "The vector pointing in the direction of anistropy for "//&
1559 "horizont viscosity. n1,n2 are the i,j components relative "//&
1560 "to the spherical coordinates.", units =
"nondim", fail_if_missing=.true.)
1564 call get_param(param_file, mdl,
"BIHARMONIC", cs%biharmonic, &
1565 "If true, use a biharmonic horizontal viscosity. "//&
1566 "BIHARMONIC may be used with LAPLACIAN.", &
1568 if (cs%biharmonic .or. get_all)
then
1569 call get_param(param_file, mdl,
"AH", ah, &
1570 "The background biharmonic horizontal viscosity.", &
1571 units =
"m4 s-1", default=0.0, scale=us%m_to_L**4*us%T_to_s)
1572 call get_param(param_file, mdl,
"AH_VEL_SCALE", ah_vel_scale, &
1573 "The velocity scale which is multiplied by the cube of "//&
1574 "the grid spacing to calculate the biharmonic viscosity. "//&
1575 "The final viscosity is the largest of this scaled "//&
1576 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
1577 units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
1578 call get_param(param_file, mdl,
"AH_TIME_SCALE", ah_time_scale, &
1579 "A time scale whose inverse is multiplied by the fourth "//&
1580 "power of the grid spacing to calculate biharmonic viscosity. "//&
1581 "The final viscosity is the largest of all viscosity "//&
1582 "formulations in use. 0.0 means that it's not used.", &
1583 units=
"s", default=0.0, scale=us%s_to_T)
1584 call get_param(param_file, mdl,
"SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
1585 "If true, use a biharmonic Smagorinsky nonlinear eddy "//&
1586 "viscosity.", default=.false.)
1587 call get_param(param_file, mdl,
"LEITH_AH", cs%Leith_Ah, &
1588 "If true, use a biharmonic Leith nonlinear eddy "//&
1589 "viscosity.", default=.false.)
1591 call get_param(param_file, mdl,
"BOUND_AH", cs%bound_Ah, &
1592 "If true, the biharmonic coefficient is locally limited "//&
1593 "to be stable.", default=.true.)
1594 call get_param(param_file, mdl,
"BETTER_BOUND_AH", cs%better_bound_Ah, &
1595 "If true, the biharmonic coefficient is locally limited "//&
1596 "to be stable with a better bounding than just BOUND_AH.", &
1597 default=cs%bound_Ah)
1599 if (cs%Smagorinsky_Ah .or. get_all)
then
1600 call get_param(param_file, mdl,
"SMAG_BI_CONST",smag_bi_const, &
1601 "The nondimensional biharmonic Smagorinsky constant, "//&
1602 "typically 0.015 - 0.06.", units=
"nondim", default=0.0, &
1603 fail_if_missing = cs%Smagorinsky_Ah)
1605 call get_param(param_file, mdl,
"BOUND_CORIOLIS", bound_cor_def, default=.false.)
1606 call get_param(param_file, mdl,
"BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &
1607 "If true use a viscosity that increases with the square "//&
1608 "of the velocity shears, so that the resulting viscous "//&
1609 "drag is of comparable magnitude to the Coriolis terms "//&
1610 "when the velocity differences between adjacent grid "//&
1611 "points is 0.5*BOUND_CORIOLIS_VEL. The default is the "//&
1612 "value of BOUND_CORIOLIS (or false).", default=bound_cor_def)
1613 if (cs%bound_Coriolis .or. get_all)
then
1614 call get_param(param_file, mdl,
"MAXVEL", maxvel, default=3.0e8)
1615 bound_cor_vel = maxvel
1616 call get_param(param_file, mdl,
"BOUND_CORIOLIS_VEL", bound_cor_vel, &
1617 "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&
1618 "the biharmonic drag to have comparable magnitude to the "//&
1619 "Coriolis acceleration. The default is set by MAXVEL.", &
1620 units=
"m s-1", default=maxvel, scale=us%m_s_to_L_T)
1624 if (cs%Leith_Ah .or. get_all) &
1625 call get_param(param_file, mdl,
"LEITH_BI_CONST", leith_bi_const, &
1626 "The nondimensional biharmonic Leith constant, "//&
1627 "typical values are thus far undetermined.", units=
"nondim", default=0.0, &
1628 fail_if_missing = cs%Leith_Ah)
1632 call get_param(param_file, mdl,
"USE_LAND_MASK_FOR_HVISC", cs%use_land_mask, &
1633 "If true, use Use the land mask for the computation of thicknesses "//&
1634 "at velocity locations. This eliminates the dependence on arbitrary "//&
1635 "values over land or outside of the domain. Default is False in order to "//&
1636 "maintain answers with legacy experiments but should be changed to True "//&
1637 "for new experiments.", default=.false.)
1639 if (cs%better_bound_Ah .or. cs%better_bound_Kh .or. get_all) &
1640 call get_param(param_file, mdl,
"HORVISC_BOUND_COEF", cs%bound_coef, &
1641 "The nondimensional coefficient of the ratio of the "//&
1642 "viscosity bounds to the theoretical maximum for "//&
1643 "stability without considering other terms.", units=
"nondim", &
1646 call get_param(param_file, mdl,
"NOSLIP", cs%no_slip, &
1647 "If true, no slip boundary conditions are used; otherwise "//&
1648 "free slip boundary conditions are assumed. The "//&
1649 "implementation of the free slip BCs on a C-grid is much "//&
1650 "cleaner than the no slip BCs. The use of free slip BCs "//&
1651 "is strongly encouraged, and no slip BCs are not used with "//&
1652 "the biharmonic viscosity.", default=.false.)
1654 call get_param(param_file, mdl,
"USE_KH_BG_2D", cs%use_Kh_bg_2d, &
1655 "If true, read a file containing 2-d background harmonic "//&
1656 "viscosities. The final viscosity is the maximum of the other "//&
1657 "terms and this background value.", default=.false.)
1659 call get_param(param_file, mdl,
"USE_GME", cs%use_GME, &
1660 "If true, use the GM+E backscatter scheme in association \n"//&
1661 "with the Gent and McWilliams parameterization.", default=.false.)
1663 if (cs%use_GME)
then
1664 call get_param(param_file, mdl,
"SPLIT", split, &
1665 "Use the split time stepping if true.", default=.true., &
1667 if (.not. split)
call mom_error(fatal,
"ERROR: Currently, USE_GME = True "// &
1668 "cannot be used with SPLIT=False.")
1670 call get_param(param_file, mdl,
"USE_MEKE", use_meke, &
1671 "If true, turns on the MEKE scheme which calculates\n"// &
1672 "a sub-grid mesoscale eddy kinetic energy budget.", &
1675 if (.not. use_meke)
call mom_error(fatal,
"ERROR: Currently, USE_GME = True "// &
1676 "cannot be used with USE_MEKE=False.")
1678 call get_param(param_file, mdl,
"GME_H0", cs%GME_h0, &
1679 "The strength of GME tapers quadratically to zero when the bathymetric "//&
1680 "depth is shallower than GME_H0.", units=
"m", scale=us%m_to_Z, &
1683 call get_param(param_file, mdl,
"GME_EFFICIENCY", cs%GME_efficiency, &
1684 "The nondimensional prefactor multiplying the GME coefficient.", &
1685 units=
"nondim", default=1.0)
1687 call get_param(param_file, mdl,
"GME_LIMITER", cs%GME_limiter, &
1688 "The absolute maximum value the GME coefficient is allowed to take.", &
1689 units=
"m2 s-1", scale=us%m_to_L**2*us%T_to_s, default=1.0e7)
1693 if (cs%bound_Kh .or. cs%bound_Ah .or. cs%better_bound_Kh .or. cs%better_bound_Ah) &
1694 call get_param(param_file, mdl,
"DT", dt, &
1695 "The (baroclinic) dynamics time step.", units=
"s", scale=us%s_to_T, &
1696 fail_if_missing=.true.)
1698 if (cs%no_slip .and. cs%biharmonic) &
1699 call mom_error(fatal,
"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
1700 "at the same time in MOM.")
1702 if (.not.(cs%Laplacian .or. cs%biharmonic))
then
1704 if ( max(g%domain%niglobal, g%domain%njglobal)>2 )
call mom_error(warning, &
1705 "hor_visc_init: It is usually a very bad idea not to use either "//&
1706 "LAPLACIAN or BIHARMONIC viscosity.")
1710 deg2rad = atan(1.0) / 45.
1712 alloc_(cs%dx2h(isd:ied,jsd:jed)) ; cs%dx2h(:,:) = 0.0
1713 alloc_(cs%dy2h(isd:ied,jsd:jed)) ; cs%dy2h(:,:) = 0.0
1714 alloc_(cs%dx2q(isdb:iedb,jsdb:jedb)) ; cs%dx2q(:,:) = 0.0
1715 alloc_(cs%dy2q(isdb:iedb,jsdb:jedb)) ; cs%dy2q(:,:) = 0.0
1716 alloc_(cs%dx_dyT(isd:ied,jsd:jed)) ; cs%dx_dyT(:,:) = 0.0
1717 alloc_(cs%dy_dxT(isd:ied,jsd:jed)) ; cs%dy_dxT(:,:) = 0.0
1718 alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0
1719 alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0
1721 if (cs%Laplacian)
then
1722 alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed)) ; cs%Kh_bg_xx(:,:) = 0.0
1723 alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0
1724 if (cs%bound_Kh .or. cs%better_bound_Kh)
then
1725 alloc_(cs%Kh_Max_xx(isd:ied,jsd:jed)) ; cs%Kh_Max_xx(:,:) = 0.0
1726 alloc_(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xy(:,:) = 0.0
1728 if (cs%Smagorinsky_Kh)
then
1729 alloc_(cs%Laplac2_const_xx(isd:ied,jsd:jed)) ; cs%Laplac2_const_xx(:,:) = 0.0
1730 alloc_(cs%Laplac2_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac2_const_xy(:,:) = 0.0
1732 if (cs%Leith_Kh)
then
1733 alloc_(cs%Laplac3_const_xx(isd:ied,jsd:jed)) ; cs%Laplac3_const_xx(:,:) = 0.0
1734 alloc_(cs%Laplac3_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac3_const_xy(:,:) = 0.0
1737 alloc_(cs%reduction_xx(isd:ied,jsd:jed)) ; cs%reduction_xx(:,:) = 0.0
1738 alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0
1740 if (cs%anisotropic)
then
1741 alloc_(cs%n1n2_h(isd:ied,jsd:jed)) ; cs%n1n2_h(:,:) = 0.0
1742 alloc_(cs%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; cs%n1n1_m_n2n2_h(:,:) = 0.0
1743 alloc_(cs%n1n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n2_q(:,:) = 0.0
1744 alloc_(cs%n1n1_m_n2n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n1_m_n2n2_q(:,:) = 0.0
1745 select case (aniso_mode)
1751 cs%dynamic_aniso = .true.
1753 call mom_error(fatal,
"MOM_hor_visc: "//&
1754 "Runtime parameter ANISOTROPIC_MODE is out of range.")
1758 if (cs%use_Kh_bg_2d)
then
1759 alloc_(cs%Kh_bg_2d(isd:ied,jsd:jed)) ; cs%Kh_bg_2d(:,:) = 0.0
1760 call get_param(param_file, mdl,
"KH_BG_2D_FILENAME", filename, &
1761 'The filename containing a 2d map of "Kh".', &
1762 default=
'KH_background_2d.nc')
1763 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1764 inputdir = slasher(inputdir)
1765 call mom_read_data(trim(inputdir)//trim(filename),
'Kh', cs%Kh_bg_2d, &
1766 g%domain, timelevel=1, scale=us%m_to_L**2*us%T_to_s)
1767 call pass_var(cs%Kh_bg_2d, g%domain)
1770 if (cs%biharmonic)
then
1771 alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0
1772 alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0
1773 alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed)) ; cs%Idxdy2u(:,:) = 0.0
1774 alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb)) ; cs%Idxdy2v(:,:) = 0.0
1776 alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed)) ; cs%Ah_bg_xx(:,:) = 0.0
1777 alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0
1778 if (cs%bound_Ah .or. cs%better_bound_Ah)
then
1779 alloc_(cs%Ah_Max_xx(isd:ied,jsd:jed)) ; cs%Ah_Max_xx(:,:) = 0.0
1780 alloc_(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_Max_xy(:,:) = 0.0
1782 if (cs%Smagorinsky_Ah)
then
1783 alloc_(cs%Biharm_const_xx(isd:ied,jsd:jed)) ; cs%Biharm_const_xx(:,:) = 0.0
1784 alloc_(cs%Biharm_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const_xy(:,:) = 0.0
1785 if (cs%bound_Coriolis)
then
1786 alloc_(cs%Biharm_const2_xx(isd:ied,jsd:jed)) ; cs%Biharm_const2_xx(:,:) = 0.0
1787 alloc_(cs%Biharm_const2_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const2_xy(:,:) = 0.0
1790 if (cs%Leith_Ah)
then
1791 alloc_(cs%biharm5_const_xx(isd:ied,jsd:jed)) ; cs%biharm5_const_xx(:,:) = 0.0
1792 alloc_(cs%biharm5_const_xy(isdb:iedb,jsdb:jedb)) ; cs%biharm5_const_xy(:,:) = 0.0
1796 do j=js-2,jeq+1 ;
do i=is-2,ieq+1
1797 cs%dx2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%dy2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)
1798 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)
1800 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
1801 cs%dx2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%dy2h(i,j) = g%dyT(i,j)*g%dyT(i,j)
1802 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)
1805 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1806 cs%reduction_xx(i,j) = 1.0
1807 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1808 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &
1809 cs%reduction_xx(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
1810 if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &
1811 (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &
1812 cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / (g%dyCu(i-1,j))
1813 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1814 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &
1815 cs%reduction_xx(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
1816 if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &
1817 (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &
1818 cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / (g%dxCv(i,j-1))
1821 do j=js-1,jeq ;
do i=is-1,ieq
1822 cs%reduction_xy(i,j) = 1.0
1823 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
1824 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &
1825 cs%reduction_xy(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
1826 if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &
1827 (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &
1828 cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / (g%dyCu(i,j+1))
1829 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
1830 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &
1831 cs%reduction_xy(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
1832 if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &
1833 (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &
1834 cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / (g%dxCv(i+1,j))
1837 if (cs%Laplacian)
then
1840 if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)
1843 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1845 grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j) + cs%dy2h(i,j))
1846 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1847 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xx(i,j) = smag_lap_const * grid_sp_h2
1848 if (cs%Leith_Kh) cs%Laplac3_const_xx(i,j) = leith_lap_const * grid_sp_h3
1850 cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))
1853 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))
1856 if (kh_sin_lat>0.)
then
1857 slat_fn = abs( sin( deg2rad * g%geoLatT(i,j) ) ) ** kh_pwr_of_sine
1858 cs%Kh_bg_xx(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xx(i,j))
1861 if (cs%bound_Kh .and. .not.cs%better_bound_Kh)
then
1863 cs%Kh_Max_xx(i,j) = kh_limit * grid_sp_h2
1864 cs%Kh_bg_xx(i,j) = min(cs%Kh_bg_xx(i,j), cs%Kh_Max_xx(i,j))
1869 do j=js-1,jeq ;
do i=is-1,ieq
1871 grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j) + cs%dy2q(i,j))
1872 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1873 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xy(i,j) = smag_lap_const * grid_sp_q2
1874 if (cs%Leith_Kh) cs%Laplac3_const_xy(i,j) = leith_lap_const * grid_sp_q3
1876 cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))
1880 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))
1883 if (kh_sin_lat>0.)
then
1884 slat_fn = abs( sin( deg2rad * g%geoLatBu(i,j) ) ) ** kh_pwr_of_sine
1885 cs%Kh_bg_xy(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xy(i,j))
1888 if (cs%bound_Kh .and. .not.cs%better_bound_Kh)
then
1890 cs%Kh_Max_xy(i,j) = kh_limit * grid_sp_q2
1891 cs%Kh_bg_xy(i,j) = min(cs%Kh_bg_xy(i,j), cs%Kh_Max_xy(i,j))
1896 if (cs%biharmonic)
then
1898 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
1899 cs%Idx2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)
1900 cs%Idxdy2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))
1902 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
1903 cs%Idx2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)
1904 cs%Idxdy2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))
1907 cs%Ah_bg_xy(:,:) = 0.0
1910 if (cs%better_bound_Ah .or. cs%bound_Ah) ah_limit = 0.3 / (dt*64.0)
1911 if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &
1912 boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))
1913 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1914 grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j)+cs%dy2h(i,j))
1915 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
1917 if (cs%Smagorinsky_Ah)
then
1918 cs%Biharm_const_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)
1919 if (cs%bound_Coriolis)
then
1920 fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &
1921 abs(g%CoriolisBu(i-1,j)), abs(g%CoriolisBu(i,j)))
1922 cs%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
1923 (fmax * boundcorconst)
1926 if (cs%Leith_Ah)
then
1927 cs%biharm5_const_xx(i,j) = leith_bi_const * (grid_sp_h3 * grid_sp_h2)
1929 cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
1930 if (ah_time_scale > 0.) cs%Ah_bg_xx(i,j) = &
1931 max(cs%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / ah_time_scale)
1932 if (cs%bound_Ah .and. .not.cs%better_bound_Ah)
then
1933 cs%Ah_Max_xx(i,j) = ah_limit * (grid_sp_h2 * grid_sp_h2)
1934 cs%Ah_bg_xx(i,j) = min(cs%Ah_bg_xx(i,j), cs%Ah_Max_xx(i,j))
1937 do j=js-1,jeq ;
do i=is-1,ieq
1938 grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j)+cs%dy2q(i,j))
1939 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
1941 if (cs%Smagorinsky_Ah)
then
1942 cs%Biharm_const_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)
1943 if (cs%bound_Coriolis)
then
1944 cs%Biharm_const2_xy(i,j) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
1945 (abs(g%CoriolisBu(i,j)) * boundcorconst)
1948 if (cs%Leith_Ah)
then
1949 cs%biharm5_const_xy(i,j) = leith_bi_const * (grid_sp_q3 * grid_sp_q2)
1952 cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
1953 if (ah_time_scale > 0.) cs%Ah_bg_xy(i,j) = &
1954 max(cs%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / ah_time_scale)
1955 if (cs%bound_Ah .and. .not.cs%better_bound_Ah)
then
1956 cs%Ah_Max_xy(i,j) = ah_limit * (grid_sp_q2 * grid_sp_q2)
1957 cs%Ah_bg_xy(i,j) = min(cs%Ah_bg_xy(i,j), cs%Ah_Max_xy(i,j))
1963 if (cs%Laplacian .and. cs%better_bound_Kh)
then
1965 do j=jsq,jeq+1 ;
do i=isq,ieq+1
1967 (cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &
1968 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
1969 (cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &
1970 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
1971 cs%Kh_Max_xx(i,j) = 0.0
1973 cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom
1975 do j=js-1,jeq ;
do i=is-1,ieq
1977 (cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &
1978 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
1979 (cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &
1980 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
1981 cs%Kh_Max_xy(i,j) = 0.0
1983 cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom
1986 call hchksum(cs%Kh_Max_xx,
"Kh_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
1987 call bchksum(cs%Kh_Max_xy,
"Kh_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
1993 if (cs%biharmonic .and. cs%better_bound_Ah)
then
1995 do j=js-1,jeq+1 ;
do i=isq-1,ieq+1
1996 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)) + &
1997 cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) + &
1998 cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
1999 cs%dx2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1)) ) )
2001 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)) + &
2002 cs%dy2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) + &
2003 cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) + &
2004 cs%dx2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1)) ) )
2006 do j=jsq-1,jeq+1 ;
do i=is-1,ieq+1
2007 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)) + &
2008 cs%dy2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j)) ) + &
2009 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)) + &
2010 cs%dx2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) )
2012 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)) + &
2013 cs%dy2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j)) ) + &
2014 cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j)) + &
2015 cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) )
2018 do j=jsq,jeq+1 ;
do i=isq,ieq+1
2021 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0u(i,j) + g%IdyCu(i-1,j)*u0u(i-1,j)) + &
2022 cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0u(i,j) + g%IdxCv(i,j-1)*v0u(i,j-1))) * &
2023 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
2025 (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0v(i,j) + g%IdyCu(i-1,j)*u0v(i-1,j)) + &
2026 cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0v(i,j) + g%IdxCv(i,j-1)*v0v(i,j-1))) * &
2027 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
2028 cs%Ah_Max_xx(i,j) = 0.0
2030 cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom
2033 do j=js-1,jeq ;
do i=is-1,ieq
2036 (cs%DX_dyBu(i,j)*(u0u(i,j+1)*g%IdxCu(i,j+1) + u0u(i,j)*g%IdxCu(i,j)) + &
2037 cs%DY_dxBu(i,j)*(v0u(i+1,j)*g%IdyCv(i+1,j) + v0u(i,j)*g%IdyCv(i,j))) * &
2038 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
2040 (cs%DX_dyBu(i,j)*(u0v(i,j+1)*g%IdxCu(i,j+1) + u0v(i,j)*g%IdxCu(i,j)) + &
2041 cs%DY_dxBu(i,j)*(v0v(i+1,j)*g%IdyCv(i+1,j) + v0v(i,j)*g%IdyCv(i,j))) * &
2042 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
2043 cs%Ah_Max_xy(i,j) = 0.0
2045 cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom
2048 call hchksum(cs%Ah_Max_xx,
"Ah_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
2049 call bchksum(cs%Ah_Max_xy,
"Ah_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
2055 cs%id_diffu = register_diag_field(
'ocean_model',
'diffu', diag%axesCuL, time, &
2056 'Zonal Acceleration from Horizontal Viscosity',
'm s-2', conversion=us%L_T2_to_m_s2)
2058 cs%id_diffv = register_diag_field(
'ocean_model',
'diffv', diag%axesCvL, time, &
2059 'Meridional Acceleration from Horizontal Viscosity',
'm s-2', conversion=us%L_T2_to_m_s2)
2061 if (cs%biharmonic)
then
2062 cs%id_Ah_h = register_diag_field(
'ocean_model',
'Ahh', diag%axesTL, time, &
2063 'Biharmonic Horizontal Viscosity at h Points',
'm4 s-1', conversion=us%L_to_m**4*us%s_to_T, &
2064 cmor_field_name=
'difmxybo', &
2065 cmor_long_name=
'Ocean lateral biharmonic viscosity', &
2066 cmor_standard_name=
'ocean_momentum_xy_biharmonic_diffusivity')
2068 cs%id_Ah_q = register_diag_field(
'ocean_model',
'Ahq', diag%axesBL, time, &
2069 'Biharmonic Horizontal Viscosity at q Points',
'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)
2072 if (cs%Laplacian)
then
2073 cs%id_Kh_h = register_diag_field(
'ocean_model',
'Khh', diag%axesTL, time, &
2074 'Laplacian Horizontal Viscosity at h Points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
2075 cmor_field_name=
'difmxylo', &
2076 cmor_long_name=
'Ocean lateral Laplacian viscosity', &
2077 cmor_standard_name=
'ocean_momentum_xy_laplacian_diffusivity')
2079 cs%id_Kh_q = register_diag_field(
'ocean_model',
'Khq', diag%axesBL, time, &
2080 'Laplacian Horizontal Viscosity at q Points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2082 if (cs%Leith_Kh)
then
2083 cs%id_vort_xy_q = register_diag_field(
'ocean_model',
'vort_xy_q', diag%axesBL, time, &
2084 'Vertical vorticity at q Points',
's-1', conversion=us%s_to_T)
2086 cs%id_div_xx_h = register_diag_field(
'ocean_model',
'div_xx_h', diag%axesTL, time, &
2087 'Horizontal divergence at h Points',
's-1', conversion=us%s_to_T)
2092 if (cs%use_GME)
then
2093 cs%id_GME_coeff_h = register_diag_field(
'ocean_model',
'GME_coeff_h', diag%axesTL, time, &
2094 'GME coefficient at h Points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2096 cs%id_GME_coeff_q = register_diag_field(
'ocean_model',
'GME_coeff_q', diag%axesBL, time, &
2097 'GME coefficient at q Points',
'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2099 cs%id_FrictWork_GME = register_diag_field(
'ocean_model',
'FrictWork_GME',diag%axesTL,time,&
2100 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', &
2101 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2)
2104 cs%id_FrictWork = register_diag_field(
'ocean_model',
'FrictWork',diag%axesTL,time,&
2105 'Integral work done by lateral friction terms', &
2106 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2)
2108 cs%id_FrictWorkIntz = register_diag_field(
'ocean_model',
'FrictWorkIntz',diag%axesT1,time, &
2109 'Depth integrated work done by lateral friction', &
2110 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%s_to_T**3*us%L_to_m**2, &
2111 cmor_field_name=
'dispkexyfo', &
2112 cmor_long_name=
'Depth integrated ocean kinetic energy dissipation due to lateral friction',&
2113 cmor_standard_name=
'ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
2115 if (cs%Laplacian .or. get_all)
then
2124 real,
intent(in) :: n1
2125 real,
intent(in) :: n2
2127 real :: recip_n2_norm
2130 recip_n2_norm = n1**2 + n2**2
2131 if (recip_n2_norm > 0.) recip_n2_norm = 1./recip_n2_norm
2133 cs%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
2134 cs%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
2135 cs%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
2136 cs%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
2142 subroutine smooth_gme(CS,G,GME_flux_h,GME_flux_q)
2146 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(inout) :: GME_flux_h
2148 real,
dimension(SZIB_(G),SZJB_(G)),
optional,
intent(inout) :: GME_flux_q
2152 real,
dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original
2153 real,
dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original
2154 real :: wc, ww, we, wn, ws
2155 integer :: i, j, k, s
2160 if (
present(gme_flux_h))
then
2161 call pass_var(gme_flux_h, g%Domain)
2162 gme_flux_h_original = gme_flux_h
2167 if (g%mask2dT(i,j)==0.) cycle
2170 ww = 0.125 * g%mask2dT(i-1,j)
2171 we = 0.125 * g%mask2dT(i+1,j)
2172 ws = 0.125 * g%mask2dT(i,j-1)
2173 wn = 0.125 * g%mask2dT(i,j+1)
2174 wc = 1.0 - (ww+we+wn+ws)
2176 gme_flux_h(i,j) = wc * gme_flux_h_original(i,j) &
2177 + ww * gme_flux_h_original(i-1,j) &
2178 + we * gme_flux_h_original(i+1,j) &
2179 + ws * gme_flux_h_original(i,j-1) &
2180 + wn * gme_flux_h_original(i,j+1)
2185 if (
present(gme_flux_q))
then
2186 call pass_var(gme_flux_q, g%Domain, position=corner, complete=.true.)
2187 gme_flux_q_original = gme_flux_q
2189 do j = g%JscB, g%JecB
2190 do i = g%IscB, g%IecB
2192 if (g%mask2dBu(i,j)==0.) cycle
2195 ww = 0.125 * g%mask2dBu(i-1,j)
2196 we = 0.125 * g%mask2dBu(i+1,j)
2197 ws = 0.125 * g%mask2dBu(i,j-1)
2198 wn = 0.125 * g%mask2dBu(i,j+1)
2199 wc = 1.0 - (ww+we+wn+ws)
2201 gme_flux_q(i,j) = wc * gme_flux_q_original(i,j) &
2202 + ww * gme_flux_q_original(i-1,j) &
2203 + we * gme_flux_q_original(i+1,j) &
2204 + ws * gme_flux_q_original(i,j-1) &
2205 + wn * gme_flux_q_original(i,j+1)
2218 if (cs%Laplacian .or. cs%biharmonic)
then
2219 dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)
2220 dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)
2221 dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)
2224 if (cs%Laplacian)
then
2225 dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)
2226 if (cs%bound_Kh)
then
2227 dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
2229 if (cs%Smagorinsky_Kh)
then
2230 dealloc_(cs%Laplac2_const_xx) ; dealloc_(cs%Laplac2_const_xy)
2232 if (cs%Leith_Kh)
then
2233 dealloc_(cs%Laplac3_const_xx) ; dealloc_(cs%Laplac3_const_xy)
2237 if (cs%biharmonic)
then
2238 dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)
2239 dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)
2240 dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)
2241 if (cs%bound_Ah)
then
2242 dealloc_(cs%Ah_Max_xx) ; dealloc_(cs%Ah_Max_xy)
2244 if (cs%Smagorinsky_Ah)
then
2245 dealloc_(cs%Biharm5_const_xx) ; dealloc_(cs%Biharm5_const_xy)
2250 if (cs%Leith_Ah)
then
2251 dealloc_(cs%Biharm_const_xx) ; dealloc_(cs%Biharm_const_xy)
2254 if (cs%anisotropic)
then
2257 dealloc_(cs%n1n1_m_n2n2_h)
2258 dealloc_(cs%n1n1_m_n2n2_q)