Skip to content

Commit

Permalink
Replaced calc_vert_vort_mag() with calc_Leith_viscosity()
Browse files Browse the repository at this point in the history 
- The subroutine that was calculating terms used in the Leith viscosities
  now returns the viscosities themselves.
  • Loading branch information
@adcroft
adcroft committed Dec 5, 2017
1 parent c59f87e commit c546d3d
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 32 deletions.
27 changes: 12 additions & 15 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ module MOM_hor_visc
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
use MOM_lateral_mixing_coeffs, only : calc_vert_vort_mag
use MOM_lateral_mixing_coeffs, only : calc_Leith_viscosity
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 @@ -282,14 +282,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
str_xx,& ! str_xx is the diagonal term in the stress tensor (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)
vert_vort_mag_h ! Magnitude of vertical vorticity at h-points |dv/dx - du/dy| (1/sec)
Leith_Kh_h, & ! Leith Laplacian viscosity at h-points (m2 s-1)
Leith_Ah_h ! Leith bi-harmonic viscosity at h-points (m4 s-1)

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain (s-1)
sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms
str_xy, & ! str_xy is the cross term in the stress tensor (H m2 s-2)
bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2)
vert_vort_mag_q ! Magnitude of vertical vorticity at q-points |dv/dx - du/dy| (1/sec)
Leith_Kh_q, & ! Leith Laplacian viscosity at q-points (m2 s-1)
Leith_Ah_q ! Leith bi-harmonic viscosity at q-points (m4 s-1)

real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
Ah_q, & ! biharmonic viscosity at corner points (m4/s)
Expand All @@ -307,7 +309,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
real :: AhLth ! 2D Leith biharmonic viscosity (m4/s)
real :: KhLth ! 2D Leith Laplacian viscosity (m2/s)
real :: Shear_mag ! magnitude of the shear (1/s)
! real :: Vort_mag ! magnitude of the vorticity (1/s)
real :: h2uq, h2vq ! temporary variables in units of H^2 (i.e. m2 or kg2 m-4).
real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner
! points; these are first interpolated to u or v velocity
Expand Down Expand Up @@ -370,12 +371,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
!$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, &
!$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
!$OMP find_FrictWork,FrictWork,use_MEKE_Ku,MEKE, &
!$OMP vert_vort_mag_h, vert_vort_mag_q) &
!$OMP Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q) &
!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
!$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
!$OMP AhLth,KhLth, &
!$OMP mod_Leith, &
!$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min)
do k=1,nz

Expand Down Expand Up @@ -545,7 +545,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
endif

if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
call calc_vert_vort_mag(VarMix, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q)
call calc_Leith_viscosity(VarMix, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q)
endif

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand All @@ -568,8 +568,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
if ((CS%Smagorinsky_Kh) .or. (CS%Leith_Kh)) then
if (CS%Smagorinsky_Kh) &
KhSm = CS%LAPLAC_CONST_xx(i,j) * Shear_mag
if (CS%Leith_Kh) &
KhLth = CS%LAPLAC3_CONST_xx(i,j) * vert_vort_mag_h(i,j)
if (CS%Leith_Kh) KhLth = Leith_Kh_h(i,j)
! Note: move Leith outside of resolution function
Kh = Kh_scale * MAX(KhLth, MAX(CS%Kh_bg_xx(i,j), KhSm))
if (CS%bound_Kh .and. .not.CS%better_bound_Kh) &
Kh = MIN(Kh, CS%Kh_Max_xx(i,j))
Expand Down Expand Up @@ -610,8 +610,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
AhSm = CS%BIHARM_CONST_xx(i,j) * Shear_mag
endif
endif
if (CS%Leith_Ah) &
AhLth = vert_vort_mag_h(i,j) * (CS%BIHARM_CONST_xx(i,j))
if (CS%Leith_Ah) AhLth = Leith_Ah_h(i,j)
Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm),AhLth)
if (CS%bound_Ah .and. .not.CS%better_bound_Ah) &
Ah = MIN(Ah, CS%Ah_Max_xx(i,j))
Expand Down Expand Up @@ -719,8 +718,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
if ((CS%Smagorinsky_Kh) .or. (CS%Leith_Kh)) then
if (CS%Smagorinsky_Kh) &
KhSm = CS%LAPLAC_CONST_xy(I,J) * Shear_mag
if (CS%Leith_Kh) &
KhLth = CS%LAPLAC3_CONST_xy(I,J) * vert_vort_mag_q(I,J)
if (CS%Leith_Kh) KhLth = Leith_Kh_q(I,J)
Kh = Kh_scale * MAX(MAX(CS%Kh_bg_xy(I,J), KhSm), KhLth)
if (CS%bound_Kh .and. .not.CS%better_bound_Kh) &
Kh = MIN(Kh, CS%Kh_Max_xy(I,J))
Expand Down Expand Up @@ -764,8 +762,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
AhSm = CS%BIHARM_CONST_xy(I,J) * Shear_mag
endif
endif
if (CS%Leith_Ah) &
AhLth = vert_vort_mag_q(I,J) * (CS%BIHARM5_CONST_xy(I,J))
if (CS%Leith_Ah) AhLth = Leith_Ah_q(I,J)
Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm),AhLth)
if (CS%bound_Ah .and. .not.CS%better_bound_Ah) &
Ah = MIN(Ah, CS%Ah_Max_xy(I,J))
Expand Down
81 changes: 64 additions & 17 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,16 +113,20 @@ module MOM_lateral_mixing_coeffs
logical :: use_QG_Leith !< If true, enables the QG Leith scheme
logical :: Leith_Kh !< If true, enables the Leith scheme
logical :: modified_Leith !< if true, include the divergence contribution to Leith viscosity
real :: Leith_Lap_const !< The non-dimensional coefficient in the Leith viscosity
logical :: Leith_Ah !< If true, enables the bi-harmonic Leith scheme
real :: Leith_bi_const !< The non-dimensional coefficient in the bi-harmonic Leith viscosity
logical :: no_slip !< If true, no slip boundary conditions are used.
!! Otherwise free slip boundary conditions are assumed.
!! The implementation of the free slip boundary
!! conditions on a C-grid is much cleaner than the
!! no slip boundary conditions. The use of free slip
!! b.c.s is strongly encouraged. The no slip b.c.s
!! are not implemented with the biharmonic viscosity.
real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
Laplac3_const_xx, & ! Laplacian metric-dependent constants (nondim)
biharm5_const_xx ! Biharmonic metric-dependent constants (nondim)
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
Laplac3_const_xy, & ! Laplacian metric-dependent constants (nondim)
biharm5_const_xy ! Biharmonic metric-dependent constants (nondim)

! Diagnostics
!>@{
Expand All @@ -140,7 +144,7 @@ module MOM_lateral_mixing_coeffs
end type VarMix_CS

public VarMix_init, calc_slope_functions, calc_resoln_function
public calc_vert_vort_mag
public calc_Leith_viscosity

contains

Expand Down Expand Up @@ -735,19 +739,19 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes)

end subroutine calc_slope_functions_using_just_e

!> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes
subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_mag_q)
!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
subroutine calc_Leith_viscosity(CS, G, GV, u, v, h, k, Leith_Kh_h, Leith_Kh_q, Leith_Ah_h, Leith_Ah_q)
type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow (m s-1)
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow (m s-1)
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: h !< Layer thickness (m or kg m-2)
integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: vert_vort_mag_h !< Magnitude of vertical component
!! of vorticity at h-ponts (s-1)
real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: vert_vort_mag_q !< Magnitude of vertical component
!! of vorticity at q-ponts (s-1)
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity at h-points (m2 s-1)
real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity at q-points (m2 s-1)
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity at h-points (m4 s-1)
real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity at q-points (m4 s-1)
! Local variables
real, dimension(SZIB_(G),SZJB_(G)) :: vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1)
dudy, & ! Meridional shear of zonal velocity (s-1)
Expand All @@ -760,7 +764,7 @@ subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_
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)
real, dimension(SZI_(G),SZJ_(G)) :: div_xx ! Estimate of horizontal divergence at h-points (s-1)
real :: mod_Leith, DY_dxBu, DX_dyBu
real :: mod_Leith, DY_dxBu, DX_dyBu, vert_vort_mag
integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand Down Expand Up @@ -828,25 +832,29 @@ subroutine calc_vert_vort_mag(CS, G, GV, u, v, h, k, vert_vort_mag_h, vert_vort_

mod_Leith = 0.; if (CS%modified_Leith) mod_Leith = 1.0

! Magnitude of vorticity at h-points
! h-point viscosities
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
vert_vort_mag_h(i,j) = sqrt( &
vert_vort_mag = sqrt( &
0.5*((vort_xy_dx(i,J-1)*vort_xy_dx(i,J-1) + vort_xy_dx(i,J)*vort_xy_dx(i,J)) + &
(vort_xy_dy(I-1,j)*vort_xy_dy(I-1,j) + vort_xy_dy(I,j)*vort_xy_dy(I,j))) + &
mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)*div_xx_dx(I-1,j)) + &
(div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1))))
if (CS%Leith_Kh) Leith_Kh_h(i,j) = CS%Laplac3_const_xx(i,j) * vert_vort_mag
if (CS%Leith_Ah) Leith_Ah_h(i,j) = CS%biharm5_const_xx(i,j) * vert_vort_mag
enddo ; enddo

! Magnitude of vorticity at q-points
! q-point viscosities
do J=js-1,Jeq ; do I=is-1,Ieq
vert_vort_mag_q(I,J) = sqrt( &
vert_vort_mag = sqrt( &
0.5*((vort_xy_dx(i,J)*vort_xy_dx(i,J) + vort_xy_dx(i+1,J)*vort_xy_dx(i+1,J)) + &
(vort_xy_dy(I,j)*vort_xy_dy(I,j) + vort_xy_dy(I,j+1)*vort_xy_dy(I,j+1))) + &
mod_Leith*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + &
(div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J))))
if (CS%Leith_Kh) Leith_Kh_q(I,J) = CS%Laplac3_const_xy(I,J) * vert_vort_mag
if (CS%Leith_Ah) Leith_Ah_q(I,J) = CS%biharm5_const_xx(I,J) * vert_vort_mag
enddo ; enddo

end subroutine calc_vert_vort_mag
end subroutine calc_Leith_viscosity

subroutine VarMix_init(Time, G, param_file, diag, CS)
type(time_type), intent(in) :: Time !< Current model time
Expand All @@ -862,6 +870,9 @@ subroutine VarMix_init(Time, G, param_file, diag, CS)
! value is roughly (pi / (the age of the universe) )^2.
logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
real :: MLE_front_length
real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity
real :: Leith_bi_const ! The non-dimensional coefficient in the bi-harmonic Leith viscosity
real :: DX2, DY2, grid_sp_2, grid_sp_3 ! Intermediate quantities for Leith metrics
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
Expand Down Expand Up @@ -1158,11 +1169,11 @@ subroutine VarMix_init(Time, G, param_file, diag, CS)
"If true, add a term to Leith viscosity which is \n"//&
"proportional to the gradient of divergence.", &
default=.false.)
call get_param(param_file, mdl, "LEITH_LAP_CONST", CS%Leith_Lap_const, &
call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
"The nondimensional Laplacian Leith constant, \n"//&
"often set to 1.0", units="nondim", default=0.0, &
fail_if_missing = CS%Leith_Kh)
call get_param(param_file, mdl, "LEITH_BI_CONST", CS%Leith_bi_const, &
call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, &
"The nondimensional biharmonic Leith constant, \n"//&
"typical values are thus far undetermined.", units="nondim", default=0.0, &
fail_if_missing = CS%Leith_Ah)
Expand All @@ -1176,6 +1187,42 @@ subroutine VarMix_init(Time, G, param_file, diag, CS)
"is strongly encouraged, and no slip BCs are not used with \n"//&
"the biharmonic viscosity.", default=.false.)
endif
if (CS%Leith_Kh) then
allocate(CS%Laplac3_const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_const_xx(:,:) = 0.0
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
DX2 = G%dxT(i,j)*G%dxT(i,j)
DY2 = G%dyT(i,j)*G%dyT(i,j)
grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2)
grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_3
enddo ; enddo
allocate(CS%Laplac3_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_const_xy(:,:) = 0.0
do J=js-1,Jeq ; do I=is-1,Ieq
DX2 = G%dxBu(I,J)*G%dxBu(I,J)
DY2 = G%dyBu(I,J)*G%dyBu(I,J)
grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2)
grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_3
enddo ; enddo
endif
if (CS%Leith_Ah) then
allocate(CS%biharm5_const_xx(isd:ied,jsd:jed)) ; CS%biharm5_const_xx(:,:) = 0.0
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
DX2 = G%dxT(i,j)*G%dxT(i,j)
DY2 = G%dyT(i,j)*G%dyT(i,j)
grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2)
grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
CS%biharm5_const_xx(i,j) = Leith_bi_const * (grid_sp_2 * grid_sp_3)
enddo ; enddo
allocate(CS%biharm5_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm5_const_xy(:,:) = 0.0
do J=js-1,Jeq ; do I=is-1,Ieq
DX2 = G%dxBu(I,J)*G%dxBu(I,J)
DY2 = G%dyBu(I,J)*G%dyBu(I,J)
grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2)
grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
CS%biharm5_const_xy(I,J) = Leith_bi_const * (grid_sp_2 * grid_sp_3)
enddo ; enddo
endif

! If nothing is being stored in this class then deallocate
if (in_use) then
Expand Down

0 comments on commit c546d3d

Please sign in to comment.