Skip to content

Commit

Permalink
Moved Leith vorticity magnitude calc into MOM_lateral_mixing_coeffs.F90
Browse files Browse the repository at this point in the history 
- Block of calculations for vorticity magnitude used in Leith are now
  in a subroutine in MOM_lateral_mixing_coeffs.F90
  • Loading branch information
@adcroft
adcroft committed Dec 4, 2017
1 parent e5d4bdf commit 273745c
Show file tree
Hide file tree
Showing 2 changed files with 109 additions and 72 deletions.
75 changes: 4 additions & 71 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
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_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 @@ -280,7 +281,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
sh_xx, & ! 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)
bhstr_xx,& ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2)
div_xx, & ! horizontal divergence (du/dx + dv/dy) (1/sec) including metric terms
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)

Expand All @@ -289,17 +289,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
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)
vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (1/sec)
vert_vort_mag_q ! Magnitude of vertical vorticity at q-points |dv/dx - du/dy| (1/sec)

real, dimension(SZI_(G),SZJB_(G)) :: &
vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms
div_xx_dy ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 sec-1) including metric terms

real, dimension(SZIB_(G),SZJ_(G)) :: &
vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms
div_xx_dx ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 sec-1) including metric terms

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)
Expand Down Expand Up @@ -386,8 +377,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
!$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 vort_xy,vort_xy_dx,vort_xy_dy,AhLth,KhLth, &
!$OMP div_xx, div_xx_dx, div_xx_dy, mod_Leith, &
!$OMP AhLth,KhLth, &
!$OMP mod_Leith, &
!$OMP Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min)
do k=1,nz

Expand All @@ -409,11 +400,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
G%IdyCu(I-1,j) * u(I-1,j,k)) - &
CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - &
G%IdxCv(i,J-1)*v(i,J-1,k)))
div_xx(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%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) + h_neglect)
enddo ; enddo
! Components for the shearing strain
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
Expand Down Expand Up @@ -535,15 +521,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
enddo ; enddo
endif

! Divergence gradient
do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j))
enddo ; enddo

do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j))
enddo ; enddo

! Coefficient for modified Leith
if (CS%Modified_Leith) then
mod_Leith = 1.0
Expand Down Expand Up @@ -578,51 +555,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
endif

if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
! Components for the vertical vorticity
! Note this a simple re-calculation of shearing components using the same discretization.
! We will consider using a circulation based calculation of vorticity later.
! Also note this will need OBC boundary conditions re-applied...
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J))
dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j))
enddo ; enddo
! Vorticity
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) )
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) )
enddo ; enddo
endif

! Vorticity gradient
do J=js-2,Jeq+1 ; do I=is-1,Ieq+1
vort_xy_dx(i,J) = CS%DY_dxBu(I,J)*(vort_xy(I,J)*G%IdyCu(I,j) - vort_xy(I-1,J)*G%IdyCu(I-1,j))
enddo ; enddo

do J=js-1,Jeq+1 ; do I=is-2,Ieq+1
vort_xy_dy(I,j) = CS%DX_dyBu(I,J)*(vort_xy(I,J)*G%IdxCv(i,J) - vort_xy(I,J-1)*G%IdxCv(i,J-1))
enddo ; enddo

! Magnitude of vorticity at h-points
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
vert_vort_mag_h(i,j) = 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))))
enddo ; enddo

! Magnitude of vorticity at q-points
do J=js-1,Jeq ; do I=is-1,Ieq
vert_vort_mag_q(I,J) = 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))))
enddo ; enddo
call calc_vert_vort_mag(G, GV, u, v, h, k, CS%no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q)
endif

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down
106 changes: 105 additions & 1 deletion src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,7 @@ module MOM_lateral_mixing_coeffs
end type VarMix_CS

public VarMix_init, calc_slope_functions, calc_resoln_function
public calc_vert_vort_mag

contains

Expand Down Expand Up @@ -719,7 +720,110 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes)

end subroutine calc_slope_functions_using_just_e

!> Initializes the variables mixing coefficients container
!> Calculates the magnitude of the vertical component of vorticity for use in the Leith-like schemes
subroutine calc_vert_vort_mag(G, GV, u, v, h, k, no_slip, mod_Leith, vert_vort_mag_h, vert_vort_mag_q)
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
logical, intent(in) :: no_slip !< True if vorticity should have no-slip BCs
real, intent(in) :: mod_Leith !< Non-dimensional coefficient multiplying the
!! divergence contribution to the Leith viscosity.
!! Set to zero for conventional Leith.
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)
! 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)
dvdx ! Zonal shear of meridional velocity (s-1)
real, dimension(SZI_(G),SZJB_(G)) :: &
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)

real, dimension(SZIB_(G),SZJ_(G)) :: &
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 :: DY_dxBu, DX_dyBu
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

! Divergence
do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
div_xx(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%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)
enddo ; enddo

! Divergence gradient
do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j))
enddo ; enddo

do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j))
enddo ; enddo

! Components for the vertical vorticity
! Note this a simple re-calculation of shearing components using the same discretization.
! We will consider using a circulation based calculation of vorticity later.
! Also note this will need OBC boundary conditions re-applied...
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
dvdx(I,J) = DY_dxBu * (v(i+1,J,k) * G%IdyCv(i+1,J) - v(i,J,k) * G%IdyCv(i,J))
DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
dudy(I,J) = DX_dyBu * (u(I,j+1,k) * G%IdxCu(I,j+1) - u(I,j,k) * G%IdxCu(I,j))
enddo ; enddo

! Vorticity
if (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) )
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) )
enddo ; enddo
endif

! Vorticity gradient
do J=js-2,Jeq+1 ; do I=is-1,Ieq+1
DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j))
enddo ; enddo

do J=js-1,Jeq+1 ; do I=is-2,Ieq+1
DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1))
enddo ; enddo

! Magnitude of vorticity at h-points
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
vert_vort_mag_h(i,j) = 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))))
enddo ; enddo

! Magnitude of vorticity at q-points
do J=js-1,Jeq ; do I=is-1,Ieq
vert_vort_mag_q(I,J) = 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))))
enddo ; enddo

end subroutine calc_vert_vort_mag

subroutine VarMix_init(Time, G, param_file, diag, CS)
type(time_type), intent(in) :: Time !< Current model time
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
Expand Down

0 comments on commit 273745c

Please sign in to comment.