Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(*)Fix discretization issues with USE_GME=True #1481

Merged
merged 2 commits into from
Aug 30, 2021
Merged

(*)Fix discretization issues with USE_GME=True #1481

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dc205c3
(*)Fix discretization issues with USE_GME=True
Hallberg-NOAA Aug 26, 2021
7ab8f12
Merge branch 'dev/gfdl' into hor_visc_GME_fix
marshallward Aug 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 44 additions & 24 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -276,6 +276,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2]
grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2]
grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2]
GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim]
htot, & ! The total thickness of all layers [Z ~> m]
boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]

real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2]
Expand All @@ -301,6 +303,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2]
! This form guarantees that hq/hu < 4.
grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2]
GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim]
boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim]

real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
Expand Down Expand Up @@ -338,6 +341,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner
! points; these are first interpolated to u or v velocity
! points where masks are applied [H ~> m or kg m-2].
real :: h_arith_q ! The arithmetic mean total thickness at q points [Z ~> m]
real :: h_harm_q ! The harmonic mean total thickness at q points [Z ~> m]
real :: I_hq ! The inverse of the arithmetic mean total thickness at q points [Z-1 ~> m-1]
real :: I_GME_h0 ! The inverse of GME tapering scale [Z-1 ~> m-1]
real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2]
real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6]
real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m]
Expand Down Expand Up @@ -501,6 +508,35 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
(0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)
enddo ; enddo

do j=js-1,je+1 ; do i=is-1,ie+1
htot(i,j) = 0.0
enddo ; enddo
do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
htot(i,j) = htot(i,j) + GV%H_to_Z*h(i,j,k)
enddo ; enddo ; enddo

I_GME_h0 = 1.0 / CS%GME_h0
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if (grad_vel_mag_bt_h(i,j)>0) then
GME_effic_h(i,j) = CS%GME_efficiency * boundary_mask_h(i,j) * &
(MIN(htot(i,j) * I_GME_h0, 1.0)**2)
else
GME_effic_h(i,j) = 0.0
endif
enddo ; enddo

do J=js-1,Jeq ; do I=is-1,Ieq
if (grad_vel_mag_bt_q(I,J)>0) then
h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1)))
I_hq = 1.0 / h_arith_q
h_harm_q = 0.25 * h_arith_q * ((htot(i,j)*I_hq + htot(i+1,j+1)*I_hq) + &
(htot(i+1,j)*I_hq + htot(i,j+1)*I_hq))
GME_effic_q(I,J) = CS%GME_efficiency * boundary_mask_q(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2)
else
GME_effic_q(I,J) = 0.0
endif
enddo ; enddo

endif ! use_GME

!$OMP parallel do default(none) &
Expand All @@ -509,7 +545,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
!$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
!$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
!$OMP use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, &
!$OMP backscat_subround, GME_coeff_limiter, &
!$OMP backscat_subround, GME_coeff_limiter, GME_effic_h, GME_effic_q, &
!$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, &
!$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, &
!$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, &
Expand Down Expand Up @@ -1388,15 +1424,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
call pass_vector(KH_u_GME, KH_v_GME, G%Domain)

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if (grad_vel_mag_bt_h(i,j)>0) then
GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * &
(0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k)))
else
GME_coeff = 0.0
endif

! apply mask
GME_coeff = GME_coeff * boundary_mask_h(i,j)
GME_coeff = GME_effic_h(i,j) * 0.25 * &
((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k)))
GME_coeff = MIN(GME_coeff, CS%GME_limiter)

if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff
Expand All @@ -1405,17 +1434,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo

do J=js-1,Jeq ; do I=is-1,Ieq
if (grad_vel_mag_bt_q(I,J)>0) then
!### This expression is not rotationally invariant - bathyT is to the SW of q points,
! and it needs parentheses in the sum of the 4 diffusivities.
GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * &
(0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k)))
else
GME_coeff = 0.0
endif

! apply mask
GME_coeff = GME_coeff * boundary_mask_q(I,J)
GME_coeff = GME_effic_q(I,J) * 0.25 * &
((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k)))
GME_coeff = MIN(GME_coeff, CS%GME_limiter)

if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
Expand All @@ -1424,8 +1444,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo

! Applying GME diagonal term. This is linear and the arguments can be rescaled.
call smooth_GME(CS,G,GME_flux_h=str_xx_GME)
call smooth_GME(CS,G,GME_flux_q=str_xy_GME)
call smooth_GME(CS, G, GME_flux_h=str_xx_GME)
call smooth_GME(CS, G, GME_flux_q=str_xy_GME)

do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j))
Expand All @@ -1446,7 +1466,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo
endif

else ! use_GME
else ! .not. use_GME
do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
enddo ; enddo
Expand Down