From dc205c3a29a1d9304010582d329f03c87f5fec33 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 26 Aug 2021 16:46:43 -0400 Subject: [PATCH] (*)Fix discretization issues with USE_GME=True Corrected the inconsistently staggered expressions for the total depths used to scale away the GM+E backscatter scheme that is enabled with USE_GME. As a part of these changes, the nominal bathymetric depths have been replaced by the sum of the layer thicknesses, and parentheses were added to some of the GME expressions so that they will be rotationally invariant. Although there are no changes to the answers in the MOM6-examples test suite, answers will change in any cases where USE_GME is true, and GME_H0 is larger than the minimum ocean depth. This commit will address MOM6 issue #1466, which can now be closed. --- .../lateral/MOM_hor_visc.F90 | 68 ++++++++++++------- 1 file changed, 44 insertions(+), 24 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index b4f857dec4..cf1712afc6 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -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] @@ -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)) :: & @@ -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] @@ -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) & @@ -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, & @@ -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 @@ -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 @@ -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)) @@ -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