Skip to content

Commit

Permalink
Add diags for Lapl. and Bihar grid Reynolds #s
Browse files Browse the repository at this point in the history 
grid_Re_Kh = (U sqtr(dx2))/Kh

grid_Re_Kh = (U dx3)/Ah

where dx2 is the harmonic mean of the squares of the grid [L2], and
dx3 is the harmonic mean of the squares of the grid^(3/2) [L3]
  • Loading branch information
@gustavo-marques
gustavo-marques committed Apr 16, 2020
1 parent 6cf28bf commit 65f36d7
Showing 1 changed file with 30 additions and 3 deletions.
33 changes: 30 additions & 3 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,9 @@ module MOM_hor_visc
Kh_Max_xx, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
Ah_Max_xx, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points
n1n1_m_n2n2_h !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points

n1n1_m_n2n2_h, & !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points
grid_sp_h2, & !< Harmonic mean of the squares of the grid [L2 ~> m2]
grid_sp_h3 !< Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Kh_bg_xy
!< The background Laplacian viscosity at q points [L2 T-1 ~> m2 s-1].
!! The actual viscosity may be the larger of this
Expand Down Expand Up @@ -180,6 +181,7 @@ module MOM_hor_visc

!>@{
!! Diagnostic id
integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh = -1
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
Expand Down Expand Up @@ -304,9 +306,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1]
max_diss_rate_h, & ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3]
FrictWork, & ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2]
FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2]
FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2]
div_xx_h ! horizontal divergence [T-1 ~> s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
grid_Re_Kh, & !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim]
grid_Re_Ah, & !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim]
GME_coeff_h !< GME coeff. at h-points [L2 T-1 ~> m2 s-1]
real :: Ah ! biharmonic viscosity [L4 T-1 ~> m4 s-1]
real :: Kh ! Laplacian viscosity [L2 T-1 ~> m2 s-1]
Expand Down Expand Up @@ -842,6 +846,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
endif

if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh

if (CS%id_grid_Re_Kh>0) then
KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2)
grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j)))/Kh
endif

if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)

str_xx(i,j) = -Kh * sh_xx(i,j)
Expand Down Expand Up @@ -890,6 +900,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

if ((CS%id_Ah_h>0) .or. find_FrictWork .or. CS%debug) Ah_h(i,j,k) = Ah

if (CS%id_grid_Re_Ah>0) then
KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2)
grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j))/Ah
endif

str_xx(i,j) = str_xx(i,j) + Ah * &
(CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - &
CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1)))
Expand Down Expand Up @@ -1295,10 +1310,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag)
if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag)
if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag)
if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag)
if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag)
if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag)
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_grid_Re_Kh>0) call post_data(CS%id_grid_Re_Kh, grid_Re_Kh, 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)
Expand Down Expand Up @@ -1658,6 +1675,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
ALLOC_(CS%dx_dyBu(IsdB:IedB,JsdB:JedB)) ; CS%dx_dyBu(:,:) = 0.0
ALLOC_(CS%dy_dxBu(IsdB:IedB,JsdB:JedB)) ; CS%dy_dxBu(:,:) = 0.0
if (CS%Laplacian) then
ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0
ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0
ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0
if (CS%bound_Kh .or. CS%better_bound_Kh) then
Expand Down Expand Up @@ -1710,6 +1728,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
ALLOC_(CS%Idxdy2v(isd:ied,JsdB:JedB)) ; CS%Idxdy2v(:,:) = 0.0
ALLOC_(CS%Ah_bg_xx(isd:ied,jsd:jed)) ; CS%Ah_bg_xx(:,:) = 0.0
ALLOC_(CS%Ah_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_bg_xy(:,:) = 0.0
ALLOC_(CS%grid_sp_h3(IsdB:IedB,JsdB:JedB)); CS%grid_sp_h3(:,:) = 0.0
if (CS%bound_Ah .or. CS%better_bound_Ah) then
ALLOC_(CS%Ah_Max_xx(isd:ied,jsd:jed)) ; CS%Ah_Max_xx(:,:) = 0.0
ALLOC_(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy(:,:) = 0.0
Expand Down Expand Up @@ -1777,6 +1796,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
! Static factors in the Smagorinsky and Leith schemes
grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j) + CS%dy2h(i,j))
CS%grid_sp_h2(i,j) = grid_sp_h2
grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
if (CS%Smagorinsky_Kh) CS%Laplac2_const_xx(i,j) = Smag_Lap_const * grid_sp_h2
if (CS%Leith_Kh) CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_h3
Expand Down Expand Up @@ -1837,6 +1857,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j)+CS%dy2h(i,j))
grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
CS%grid_sp_h3(i,j) = grid_sp_h3
if (CS%Smagorinsky_Ah) then
CS%Biharm_const_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2)
if (CS%bound_Coriolis) then
Expand Down Expand Up @@ -1979,6 +2000,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')
CS%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, Time, &
'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T)
CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, &
'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim')
endif
if (CS%Laplacian) then
CS%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, Time, &
Expand All @@ -1988,6 +2011,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')
CS%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, Time, &
'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
CS%id_grid_Re_Kh = register_diag_field('ocean_model', 'grid_Re_Kh', diag%axesTL, Time, &
'Grid Reynolds number for the Laplacian horizontal viscosity at h points', 'nondim')
if (CS%Leith_Kh) then
CS%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, Time, &
'Vertical vorticity at q Points', 's-1', conversion=US%s_to_T)
Expand Down Expand Up @@ -2105,6 +2130,7 @@ subroutine hor_visc_end(CS)
endif
if (CS%Laplacian) then
DEALLOC_(CS%Kh_bg_xx) ; DEALLOC_(CS%Kh_bg_xy)
DEALLOC_(CS%grid_sp_h2)
if (CS%bound_Kh) then
DEALLOC_(CS%Kh_Max_xx) ; DEALLOC_(CS%Kh_Max_xy)
endif
Expand All @@ -2116,6 +2142,7 @@ subroutine hor_visc_end(CS)
endif
endif
if (CS%biharmonic) then
DEALLOC_(CS%grid_sp_h3)
DEALLOC_(CS%Idx2dyCu) ; DEALLOC_(CS%Idx2dyCv)
DEALLOC_(CS%Idxdy2u) ; DEALLOC_(CS%Idxdy2v)
DEALLOC_(CS%Ah_bg_xx) ; DEALLOC_(CS%Ah_bg_xy)
Expand Down

0 comments on commit 65f36d7

Please sign in to comment.