Skip to content

Commit

Permalink
Adds GME and smoothing function
Browse files Browse the repository at this point in the history
  • Loading branch information
@gustavo-marques
gustavo-marques committed Jan 18, 2019
1 parent a38fc07 commit d5f5d32
Showing 1 changed file with 211 additions and 13 deletions.
224 changes: 211 additions & 13 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,13 @@ module MOM_hor_visc

use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_domains, only : pass_var
use MOM_domains, only : pass_var, CORNER
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_Leith_viscosity
use MOM_barotropic, only : barotropic_CS, barotropic_get_tav
use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH
use MOM_MEKE_types, only : MEKE_type
use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE
Expand Down Expand Up @@ -73,6 +75,7 @@ module MOM_hor_visc
real :: Kh_aniso !< The anisotropic viscosity in m2 s-1.
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 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 @@ -159,6 +162,7 @@ module MOM_hor_visc
integer :: id_diffu = -1, id_diffv = -1
integer :: id_Ah_h = -1, id_Ah_q = -1
integer :: id_Kh_h = -1, id_Kh_q = -1
integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
integer :: id_vort_xy_q = -1, id_div_xx_h = -1
integer :: id_FrictWork = -1, id_FrictWorkIntz = -1
!!@}
Expand All @@ -179,7 +183,9 @@ module MOM_hor_visc
!! u[is-2:ie+2,js-2:je+2]
!! v[is-2:ie+2,js-2:je+2]
!! h[is-1:ie+1,js-1:je+1]
subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, OBC)

subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, &
thickness_diffuse, G, GV, CS, OBC)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
Expand All @@ -199,24 +205,34 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
!! related to Mesoscale Eddy Kinetic Energy.
type(VarMix_CS), pointer :: VarMix !< Pointer to a structure with fields that
!! specify the spatially variable viscosities
type(hor_visc_CS), pointer :: CS !< Pontrol structure returned by a previous
type(barotropic_CS), pointer :: Barotropic !< Pointer to a structure containing
!! barotropic velocities
type(thickness_diffuse_CS), pointer :: thickness_diffuse !< Pointer to a structure containing
!! interface height diffusivities
type(hor_visc_CS), pointer :: CS !< Control structure returned by a previous
!! call to hor_visc_init.
type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type

! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: &
u0, & ! Laplacian of u (m-1 s-1)
h_u, & ! Thickness interpolated to u points, in H.
vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1)
div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1)
div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1)
ubtav ! zonal barotropic vel. ave. over baroclinic time-step (m s-1)
real, dimension(SZI_(G),SZJB_(G)) :: &
v0, & ! Laplacian of v (m-1 s-1)
h_v, & ! Thickness interpolated to v points, in H.
vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1)
div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1)
div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1)
vbtav ! meridional barotropic vel. ave. over baroclinic time-step (m s-1)
real, dimension(SZI_(G),SZJ_(G)) :: &
dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension (s-1)
div_xx, & ! Estimate of horizontal divergence at h-points (s-1)
sh_xx, & ! horizontal tension (du/dx - dv/dy) (1/sec) including metric terms
sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) (1/sec) including metric terms
str_xx,& ! str_xx is the diagonal term in the stress tensor (H m2 s-2)
str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME (H m2 s-2)
bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2)
FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2)
Leith_Kh_h, & ! Leith Laplacian viscosity at h-points (m2 s-1)
Expand All @@ -227,8 +243,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain (s-1)
dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain (s-1)
sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms
sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) (1/sec) inc. metric terms
str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2)
str_xy_GME, & ! smoothed cross term in the stress tensor from GME (H m2 s-2)
bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2)
vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1)
Leith_Kh_q, & ! Leith Laplacian viscosity at q-points (m2 s-1)
Expand All @@ -238,15 +257,22 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
grad_div_mag_q ! Magnitude of divergence gradient at q-points (m-1 s-1)

real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
Ah_q, & ! biharmonic viscosity at corner points (m4/s)
Kh_q, & ! Laplacian viscosity at corner points (m2/s)
vort_xy_q ! vertical vorticity at corner points (s-1)

Ah_q, & ! biharmonic viscosity at corner points (m4/s)
Kh_q, & ! Laplacian viscosity at corner points (m2/s)
vort_xy_q, & ! vertical vorticity at corner points (s-1)
GME_coeff_q !< GME coeff. at q-points (m2 s-1)

real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
KH_u_GME !< interface height diffusivities in u-columns (m2 s-1)
real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
KH_v_GME !< interface height diffusivities in v-columns (m2 s-1)
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
Ah_h, & ! biharmonic viscosity at thickness points (m4/s)
Kh_h, & ! Laplacian viscosity at thickness points (m2/s)
FrictWork, & ! energy dissipated by lateral friction (W/m2)
div_xx_h ! horizontal divergence (s-1)
div_xx_h, & ! horizontal divergence (s-1)
KH_t_GME, & !< interface height diffusivities in t-columns (m2 s-1)
GME_coeff_h !< GME coeff. at h-points (m2 s-1)

real :: Ah ! biharmonic viscosity (m4/s)
real :: Kh ! Laplacian viscosity (m2/s)
Expand Down Expand Up @@ -278,6 +304,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
real :: FatH ! abs(f) at h-point for MEKE source term (s-1)
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 :: DY_dxBu, DX_dyBu
logical :: rescale_Kh, legacy_bound
logical :: find_FrictWork
Expand Down Expand Up @@ -330,12 +357,46 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
!$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, &
!$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) &
!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
!$OMP sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, &
!$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
!$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
!$OMP div_xx, div_xx_dx, div_xx_dy, &
!$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min)


if (CS%use_GME) then
call barotropic_get_tav(Barotropic, ubtav, vbtav, G)

do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - &
G%IdyCu(I-1,j) * ubtav(I-1,j))
dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - &
G%IdxCv(i,J-1) * vbtav(i,J-1))
sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
enddo ; enddo

! Components for the barotropic shearing strain
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
dvdx_bt(I,J) = CS%DY_dxBu(I,J)*(vbtav(i+1,J)*G%IdyCv(i+1,J) &
- vbtav(i,J)*G%IdyCv(i,J))
dudy_bt(I,J) = CS%DX_dyBu(I,J)*(ubtav(I,j+1)*G%IdxCu(I,j+1) &
- ubtav(I,j)*G%IdxCu(I,j))
enddo ; enddo

if (CS%no_slip) then
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
enddo ; enddo
else
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
enddo ; enddo
endif

endif ! use_GME

do k=1,nz

! The following are the forms of the horizontal tension and horizontal
Expand Down Expand Up @@ -636,8 +697,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,

endif ! CS%Leith_Kh



do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
Expand Down Expand Up @@ -691,6 +750,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
str_xx(i,j) = 0.0
endif ! Laplacian

if (CS%use_GME) then
call thickness_diffuse_get_KH(thickness_diffuse, KH_t_GME, KH_u_GME, KH_v_GME, G)
GME_coeff = KH_t_GME(i,j,k) * &
0.5*(VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)) * &
( (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)

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

endif

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

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 @@ -740,6 +816,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
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)
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 @@ -854,6 +938,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
str_xy(I,J) = 0.0
endif ! Laplacian

if (CS%use_GME) then
GME_coeff = ( 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*(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.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(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))

str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J)
endif

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

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 @@ -902,6 +1003,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
endif
enddo ; enddo

! applying GME diagonal term
if (CS%use_GME) then
call smooth_GME(CS,G,GME_flux_q=str_xy_GME)
do J=js-1,Jeq ; do I=is-1,Ieq
if (CS%no_slip) then
str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq * CS%reduction_xy(I,J))
else
str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
endif
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 @@ -1021,6 +1134,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag)
if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag)
if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag)
if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag)
if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag)

if (CS%id_FrictWorkIntz > 0) then
do j=js,je
Expand Down Expand Up @@ -1309,6 +1424,9 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS)
"viscosities. The final viscosity is the maximum of the other "//&
"terms and this background value.", default=.false.)

call get_param(param_file, mdl, "USE_GME", CS%use_GME, &
"If true, use the GM+E backscatter scheme in association \n"//&
"with the Gent and McWilliams parameterization.", default=.false.)

if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) &
call get_param(param_file, mdl, "DT", dt, &
Expand Down Expand Up @@ -1696,6 +1814,14 @@ subroutine hor_visc_init(Time, G, param_file, diag, CS)

endif

if (CS%use_GME) then
CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, &
'GME coefficient at h Points', 'm^2 s-1')

CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, &
'GME coefficient at q Points', 'm^2 s-1')
endif

CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,&
'Integral work done by lateral friction terms', 'W m-2')

Expand Down Expand Up @@ -1730,6 +1856,78 @@ subroutine align_aniso_tensor_to_grid(CS, n1, n2)

end subroutine align_aniso_tensor_to_grid

!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any
!! horizontal two-grid-point noise
subroutine smooth_GME(CS,G,GME_flux_h,GME_flux_q)
! Arguments
type(hor_visc_CS), pointer :: CS !< Control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid
real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!<
real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!<

! local variables
real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original
real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original
real :: wc, ww, we, wn, ws ! averaging weights for smoothing
integer :: i, j, k, s

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

! Update halos
if (present(GME_flux_h)) then
call pass_var(GME_flux_h, G%Domain)
GME_flux_h_original = GME_flux_h
! apply smoothing on GME
do j = G%jsc, G%jec
do i = G%isc, G%iec
! skip land points
if (G%mask2dT(i,j)==0.) cycle

! compute weights
ww = 0.125 * G%mask2dT(i-1,j)
we = 0.125 * G%mask2dT(i+1,j)
ws = 0.125 * G%mask2dT(i,j-1)
wn = 0.125 * G%mask2dT(i,j+1)
wc = 1.0 - (ww+we+wn+ws)

GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) &
+ ww * GME_flux_h_original(i-1,j) &
+ we * GME_flux_h_original(i+1,j) &
+ ws * GME_flux_h_original(i,j-1) &
+ wn * GME_flux_h_original(i,j+1)
enddo; enddo
endif

! Update halos
if (present(GME_flux_q)) then
call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true.)
GME_flux_q_original = GME_flux_q
! apply smoothing on GME
do J = G%JscB, G%JecB
do I = G%IscB, G%IecB
! skip land points
if (G%mask2dBu(I,J)==0.) cycle

! compute weights
ww = 0.125 * G%mask2dBu(I-1,J)
we = 0.125 * G%mask2dBu(I+1,J)
ws = 0.125 * G%mask2dBu(I,J-1)
wn = 0.125 * G%mask2dBu(I,J+1)
wc = 1.0 - (ww+we+wn+ws)

GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) &
+ ww * GME_flux_q_original(I-1,J) &
+ we * GME_flux_q_original(I+1,J) &
+ ws * GME_flux_q_original(I,J-1) &
+ wn * GME_flux_q_original(I,J+1)
enddo; enddo
endif

enddo ! s-loop

end subroutine smooth_GME

!> Deallocates any variables allocated in hor_visc_init.
subroutine hor_visc_end(CS)
type(hor_visc_CS), pointer :: CS !< The control structure returned by a
Expand Down

0 comments on commit d5f5d32

Please sign in to comment.