Skip to content

Commit

Permalink
Refactored MOM_hor_visc. Moved GME out of loops where downgradient
Browse files Browse the repository at this point in the history 
viscosity is calculated. Moved calculation of FrictWork ahead of GME
so I can separate the work done by downgradient stuff from that of GME.
  • Loading branch information
@sdbachman
sdbachman committed Apr 12, 2019
1 parent bd98d83 commit 807072b
Showing 1 changed file with 136 additions and 62 deletions.
198 changes: 136 additions & 62 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ module MOM_hor_visc
logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function
!! of state. This is set depending on ANISOTROPIC_MODE.
logical :: use_GME !< If true, use GME backscatter scheme.
real :: GME_dt
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx
!< The background Laplacian viscosity at h points, in units
!! of m2 s-1. The actual viscosity may be the larger of this
Expand Down Expand Up @@ -240,7 +241,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
beta_h, & ! Gradient of planetary vorticity at h-points (m-1 s-1)
grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points (m-1 s-1)
grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points (m-1 s-1)
grad_div_mag_h ! Magnitude of divergence gradient at h-points (m-1 s-1)
grad_div_mag_h, & ! Magnitude of divergence gradient at h-points (m-1 s-1)
dudx, &
dvdy

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain (s-1)
Expand All @@ -257,8 +260,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points (m-1 s-1)
grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points (m-1 s-1)
grad_div_mag_q, & ! Magnitude of divergence gradient at q-points (m-1 s-1)
hq ! harmonic mean of the harmonic means of the u- & v-
! point thicknesses, in H; This form guarantees that hq/hu < 4.
hq ! harmonic mean of the harmonic means of the u- & v point thicknesses, in H; This form guarantees that hq/hu < 4.

real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
Ah_q, & ! biharmonic viscosity at corner points (m4/s)
Expand Down Expand Up @@ -310,8 +312,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
real :: local_strain ! Local variable for interpolating computed strain rates (s-1).
real :: epsilon
real :: GME_coeff ! The GME (negative) viscosity coefficient (m2 s-1)
real :: GME_coeff_limiter
real :: DY_dxBu, DX_dyBu
real :: H0
real :: tmp
logical :: rescale_Kh, legacy_bound
logical :: find_FrictWork
logical :: apply_OBC = .false.
Expand All @@ -326,7 +330,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
h_neglect3 = h_neglect**3
inv_PI3 = 1.0/((4.0*atan(1.0))**3)
inv_PI6 = inv_PI3**2
epsilon = 1.e-15
epsilon = 1.e-8

if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then
apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally
Expand Down Expand Up @@ -373,7 +377,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
if (CS%use_GME) then
! GME tapers off above this depth
H0 = 1000.0

! initialize diag. array with zeros
GME_coeff_h(:,:,:) = 0.0
GME_coeff_q(:,:,:) = 0.0
Expand Down Expand Up @@ -631,10 +635,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
if (CS%no_slip) then
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
dudy(I,J) = (2.0-G%mask2dBu(I,J)) * dudy(I,J)
dvdx(I,J) = (2.0-G%mask2dBu(I,J)) * dvdx(I,J)
enddo ; enddo
else
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
dudy(I,J) = G%mask2dBu(I,J) * dudy(I,J)
dvdx(I,J) = G%mask2dBu(I,J) * dvdx(I,J)
enddo ; enddo
endif

Expand All @@ -661,6 +669,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
(G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - &
G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k))))*G%IareaT(i,j)/ &
(h(i,j,k) + GV%H_subroundoff)
dudx(i,j) = 0.5*(G%dyCu(I,j) * u(I,j,k) * (h(i+1,j,k)+h(i,j,k)) - &
G%dyCu(I-1,j) * u(I-1,j,k) * (h(i-1,j,k)+h(i,j,k)) )*G%IareaT(i,j)/ &
(h(i,j,k) + GV%H_subroundoff) * G%mask2dcu(I,j)
dvdy(i,j) = 0.5*(G%dxCv(i,J) * v(i,J,k) * (h(i,j,k)+h(i,j+1,k)) - &
G%dxCv(i,J-1)*v(i,J-1,k)*(h(i,j,k)+h(i,j-1,k)))*G%IareaT(i,j)/ &
(h(i,j,k) + GV%H_subroundoff) * G%mask2dcv(i,J)
enddo ; enddo

call pass_var(div_xx, G%Domain, complete=.true.)
Expand Down Expand Up @@ -801,28 +815,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
str_xx(i,j) = 0.0
endif ! Laplacian

if (CS%use_GME) then
GME_coeff = MIN(G%bathyT(i,j)/H0,1.0)*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)) * &
0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0) * &
( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + &
(0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 ) / &
( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
(0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
(0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
epsilon)
! apply mask
GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

! simple way to limit this coeff
GME_coeff = MIN(GME_coeff,1.0E5)

if (CS%id_GME_coeff_h>0) GME_coeff_h(i,j,k) = GME_coeff

str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j)

endif ! CS%use_GME

if (CS%anisotropic) then
! Shearing-strain averaged to h-points
local_strain = 0.25 * ( (sh_xy(I,J) + sh_xy(I-1,J-1)) + (sh_xy(I-1,J) + sh_xy(I,J-1)) )
Expand Down Expand Up @@ -871,18 +863,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

enddo ; enddo

! applying GME diagonal term
if (CS%use_GME) then
call smooth_GME(CS,G,GME_flux_h=str_xx_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))
enddo ; enddo
else
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
endif

if (CS%biharmonic) then
! Gradient of Laplacian, for use in bi-harmonic term
do J=js-1,Jeq ; do I=is-1,Ieq
Expand Down Expand Up @@ -998,28 +978,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
str_xy(I,J) = 0.0
endif ! Laplacian

if (CS%use_GME) then
GME_coeff = MIN(G%bathyT(i,j)/H0,1.0)*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)) * &
0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0) * &
( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I,j+1,k)) )**2 + &
(0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i+1,J,k)) )**2 ) / &
( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
(0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
(0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
epsilon)
! apply mask
GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

! simple way to limit this coeff
GME_coeff = MIN(GME_coeff,1.0E5)

if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J)

endif

if (CS%anisotropic) then
! Horizontal-tension averaged to q-points
local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )
Expand Down Expand Up @@ -1063,6 +1021,115 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

enddo ; enddo


if (find_FrictWork) then ; do j=js,je ; do i=is,ie
! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
FrictWork(i,j,k) = GV%H_to_kg_m2 * ( &
(str_xx(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) &
-str_xx(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) &
+0.25*((str_xy(I,J)*( &
(u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) &
+(v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) &
+str_xy(I-1,J-1)*( &
(u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) &
+(v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1) )) &
+(str_xy(I-1,J)*( &
(u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) &
+(v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J) ) &
+str_xy(I,J-1)*( &
(u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) &
+(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) )
enddo ; enddo ; endif


if (CS%use_GME) then

do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1

! GME_coeff = (MIN(G%bathyT(i,j)/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)) * &
! sqrt(0.5*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)),0.0) * &
! ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I-1,j,k)) )**2 + &
! (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i,J-1,k)) )**2 ) / &
! ( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
! (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
! (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
! epsilon))

GME_coeff = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / &
SQRT( dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
(0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
(0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + &
epsilon)

! apply mask
GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

! GME_coeff_limiter = 1e5
! GME_coeff_limiter = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / sqrt(dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
! (0.25*(dvdx_bt(I,J)+dvdx_bt(I-1,J)+dvdx_bt(I,J-1)+dvdx_bt(I-1,J-1)) )**2 + &
! (0.25*(dudy_bt(I,J)+dudy_bt(I-1,J)+dudy_bt(I,J-1)+dudy_bt(I-1,J-1)) )**2 + epsilon)

! simple way to limit this coeff
! GME_coeff = MIN(GME_coeff,GME_coeff_limiter)

if (CS%id_GME_coeff_h>0) GME_coeff_h(i,j,k) = GME_coeff

str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j)

enddo ; enddo

endif ! CS%use_GME

! applying GME diagonal term
if (CS%use_GME) then
call smooth_GME(CS,G,GME_flux_h=str_xx_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))
enddo ; enddo
else
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
endif

if (CS%use_GME) then
do J=js-1,Jeq ; do I=is-1,Ieq

! GME_coeff = (MIN(G%bathyT(i,j)/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)) * &
! sqrt( 0.25*MAX((VarMix%N2_u(I,j,k)+VarMix%N2_u(I,j+1,k) + &
! VarMix%N2_v(i,J,k)+VarMix%N2_v(i+1,J,k)),0.0) * &
! ( (0.5*(VarMix%slope_x(I,j,k)+VarMix%slope_x(I,j+1,k)) )**2 + &
! (0.5*(VarMix%slope_y(i,J,k)+VarMix%slope_y(i+1,J,k)) )**2 ) / &
! ( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
! (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
! (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
! epsilon))

GME_coeff = 2.0* MAX(0.0,MEKE%MEKE(i,j)) / &
SQRT( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
(0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
(0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + &
epsilon)

! apply mask
GME_coeff = GME_coeff * (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))

! GME_coeff_limiter = 1e5
! GME_coeff_limiter = 2.0 * MAX(0.0,MEKE%MEKE(i,j)) / sqrt( dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
! (0.25*(dudx_bt(i,j)+dudx_bt(i+1,j)+dudx_bt(i,j+1)+dudx_bt(i+1,j+1)))**2 + &
! (0.25*(dvdy_bt(i,j)+dvdy_bt(i+1,j)+dvdy_bt(i,j+1)+dvdy_bt(i+1,j+1)) )**2 + epsilon)

! simple way to limit this coeff
! GME_coeff = MIN(GME_coeff,GME_coeff_limiter)

if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J)

enddo ; enddo
endif

! applying GME diagonal term
if (CS%use_GME) then
call smooth_GME(CS,G,GME_flux_q=str_xy_GME)
Expand All @@ -1087,6 +1154,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
enddo ; enddo
endif






! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
do j=js,je ; do I=Isq,Ieq
diffu(I,j,k) = ((G%IdyCu(I,j)*(CS%DY2h(i,j) *str_xx(i,j) - &
Expand Down Expand Up @@ -1644,6 +1716,8 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS)
CS%reduction_xy(I,J) = G%dx_Cv(i+1,J) / G%dxCv(i+1,J)
enddo ; enddo

CS%GME_dt = dt

if (CS%Laplacian) then
! The 0.3 below was 0.4 in MOM1.10. The change in hq requires
! this to be less than 1/3, rather than 1/2 as before.
Expand Down Expand Up @@ -1946,7 +2020,7 @@ subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q)
integer :: i, j, k, s

!do s=1,CS%n_smooth
do s=1,5
do s=1,1

! Update halos
if (present(GME_flux_h)) then
Expand Down

0 comments on commit 807072b

Please sign in to comment.