diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 451c366409..a12721a2ea 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -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 @@ -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) @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7f981f7c04..1e5e7fdbea 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -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 @@ -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