Skip to content

Commit

Permalink
Fixed problems with FrictWorkMax and FrictWorkDiss in MOM_hor_visc.F.…
Browse files Browse the repository at this point in the history 
… More specifically,

I moved the calculation of dudx and dvdy outside of the Leith loop so now it is
being calculated with sh_xx.
  • Loading branch information
@sdbachman
sdbachman committed Apr 17, 2019
1 parent 8df6efd commit 172ae02
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 43 deletions.
30 changes: 15 additions & 15 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,7 +384,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
if (CS%use_GME) then
! GME tapers off above this depth
H0 = 1000.0
FWfrac = 0.1
FWfrac = 1.0
! initialize diag. array with zeros
GME_coeff_h(:,:,:) = 0.0
GME_coeff_q(:,:,:) = 0.0
Expand Down Expand Up @@ -447,10 +447,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

! Calculate horizontal tension
do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
sh_xx(i,j) = (CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - &
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)))
dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - &
G%IdyCu(I-1,j) * u(I-1,j,k))
dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - &
G%IdxCv(i,J-1) * v(i,J-1,k))
sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
enddo ; enddo

! Components for the shearing strain
Expand All @@ -459,6 +460,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
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


! Interpolate the thicknesses to velocity points.
! The extra wide halos are to accommodate the cross-corner-point projections
! in OBCs, which are not ordinarily be necessary, and might not be necessary
Expand Down Expand Up @@ -676,12 +678,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
(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)
dudx(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%IareaT(i,j)/ &
(h(i,j,k) + GV%H_subroundoff) * G%mask2dcu(I,j)
dvdy(i,j) = 0.5*(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) * G%mask2dcv(i,J)
enddo ; enddo

call pass_var(div_xx, G%Domain, complete=.true.)
Expand Down Expand Up @@ -1031,17 +1027,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,

if (find_FrictWork) then
if (CS%biharmonic) call pass_vector(u0, v0, G%Domain)
call pass_var(dudx, G%Domain, complete=.true.)
call pass_var(dvdy, G%Domain, complete=.true.)
call pass_var(dvdx, G%Domain, position=CORNER, complete=.true.)
call pass_var(dudy, G%Domain, position=CORNER, complete=.true.)

do j=js,je ; do i=is,ie
! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
FrictWork_diss(i,j,k) = Kh_h(i,j,k) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
FrictWork_diss(i,j,k) = -Kh_h(i,j,k) * (dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + &
(0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) + &
(0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2) - &
Ah_h(i,j,k) * ((0.5*(u0(I,j) + u0(I-1,j)))**2 + &
(0.5*(v0(i,J) + v0(i,J-1)))**2)

if (associated(MEKE)) then ; if (associated(MEKE%mom_src)) then
FrictWorkMax(i,j,k) = MEKE%MEKE(i,j) * sqrt(dudx(i,j)**2 + dvdy(i,j)**2 + &
FrictWorkMax(i,j,k) = 2.0*MEKE%MEKE(i,j) * sqrt(dudx(i,j)**2 + dvdy(i,j)**2 + &
(0.25*(dvdx(I,J)+dvdx(I-1,J)+dvdx(I,J-1)+dvdx(I-1,J-1)) )**2 + &
(0.25*(dudy(I,J)+dudy(I-1,J)+dudy(I,J-1)+dudy(I-1,J-1)) )**2)

Expand All @@ -1050,7 +1050,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
! we need to add the FrictWork done by the dissipation operators, since this work
! is done only for numerical stability and is therefore spurious
if (CS%use_GME) then
target_FrictWork_GME(i,j,k) = FWfrac * FrictWorkMax(i,j,k) - FrictWork(i,j,k)
target_FrictWork_GME(i,j,k) = FWfrac * FrictWorkMax(i,j,k) - FrictWork_diss(i,j,k)
endif

endif ; endif
Expand Down
31 changes: 3 additions & 28 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -329,16 +329,6 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo
!$OMP end parallel

! if (associated(MEKE)) then
! do j=js,je ; do i=is,ie
! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + &
! KH_v(i,J,k) + KH_v(i,J-1,k))
!! 0.25*MAX(VarMix%N2_u(I,j,k)+VarMix%N2_u(I-1,j,k)+VarMix%N2_v(i,J,k)+VarMix%N2_v(i,J-1,k),0.0) * &
!! ( (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 )
! enddo; enddo
! endif

if (CS%detangle_interfaces) then
call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, &
CS, int_slope_u, int_slope_v)
Expand All @@ -362,8 +352,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS
! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S
if (use_stored_slopes) then
call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, &
int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, &
VarMix%N2_u, VarMix%N2_v)
int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y)
else
call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, CS, &
int_slope_u, int_slope_v)
Expand Down Expand Up @@ -457,7 +446,7 @@ end subroutine thickness_diffuse
!! Fluxes are limited to give positive definite thicknesses.
!! Called by thickness_diffuse().
subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, MEKE, &
CS, int_slope_u, int_slope_v, slope_x, slope_y, N2_x, N2_y)
CS, int_slope_u, int_slope_v, slope_x, slope_y)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness in H (m or kg/m2)
Expand All @@ -483,10 +472,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
!! density gradients.
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points
real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: N2_x !< Brunt-Vaisala frequency at
!! u points (s-2)
real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: N2_y !< Brunt-Vaisala frequency at
!! v points (s-2)
! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
T, & ! The temperature (or density) in C, with the values in
Expand Down Expand Up @@ -571,7 +556,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics
logical :: present_int_slope_u, present_int_slope_v
logical :: present_slope_x, present_slope_y, calc_derivatives
logical :: present_N2_x, present_N2_y
integer :: is, ie, js, je, nz, IsdB
integer :: i, j, k
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; IsdB = G%IsdB
Expand All @@ -589,8 +573,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
present_int_slope_v = PRESENT(int_slope_v)
present_slope_x = PRESENT(slope_x)
present_slope_y = PRESENT(slope_y)
present_N2_x = PRESENT(N2_x)
present_N2_y = PRESENT(N2_y)

nk_linear = max(GV%nkml, 1)

Expand Down Expand Up @@ -1198,14 +1180,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J)))
if (associated(CS%GMwork)) CS%GMwork(i,j) = Work_h
if (associated(MEKE)) then ; if (associated(MEKE%GM_src)) then
MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
! if (present_N2_x .and. present_N2_y) &
! MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + 0.25*(Kh_u(I,j,k) + KH_u(I-1,j,k) + &
! KH_v(i,J,k) + KH_v(i,J-1,k)) * &
! 0.25*MAX(N2_x(I,j,k)+N2_x(I-1,j,k)+N2_y(i,J,k)+N2_y(i,J-1,k),0.0) * &
! ( (0.5*(slope_x(I,j,k)+slope_x(I-1,j,k)) )**2 + &
! (0.5*(slope_y(i,J,k)+slope_y(i,J-1,k)) )**2 )

MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + MIN(0.0,Work_h)
endif ; endif
enddo ; enddo ; endif

Expand Down

0 comments on commit 172ae02

Please sign in to comment.