Skip to content

Commit

Permalink
Merge pull request #134 from gustavo-marques/fix_leith
Browse files Browse the repository at this point in the history 
Fix bugs in Leith add new input parameter
  • Loading branch information
@alperaltuntas
alperaltuntas authored Dec 12, 2019
2 parents 3d17e03 + 78bb4c1 commit 2cccc83
Showing 1 changed file with 26 additions and 22 deletions.
48 changes: 26 additions & 22 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,8 @@ module MOM_hor_visc
!! Default is False to maintain answers with legacy experiments
!! but should be changed to True for new experiments.
logical :: anisotropic !< If true, allow anisotropic component to the viscosity.
logical :: add_LES_viscosity!< If true, adds the viscosity from Smagorinsky and Leith to
!! the background viscosity instead of taking the maximum.
real :: Kh_aniso !< The anisotropic viscosity [L2 T-1 ~> m2 s-1].
logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function
!! of state. This is set depending on ANISOTROPIC_MODE.
Expand Down Expand Up @@ -679,8 +681,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
enddo ; enddo
endif

call pass_var(vort_xy, G%Domain, position=CORNER, complete=.true.)

! Vorticity gradient
do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
Expand All @@ -692,23 +692,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
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

call pass_vector(vort_xy_dy, vort_xy_dx, G%Domain)

if (CS%modified_Leith) then
! 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)
div_xx(i,j) = dudx(i,j) + dvdy(i,j)
enddo ; enddo

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

Expand All @@ -717,7 +711,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + &
(0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2)
enddo ; enddo
do J=js-1,Jeq ; do I=is-1,Ieq
!do J=js-1,Jeq ; do I=is-1,Ieq
do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1
grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + &
(0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2)
enddo ; enddo
Expand All @@ -727,13 +722,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
div_xx_dx(I,j) = 0.0
enddo ; enddo
do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2
div_xx_dy(i,J) = 0.0
enddo ; enddo
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
grad_div_mag_h(i,j) = 0.0
enddo ; enddo
do J=js-1,Jeq ; do I=is-1,Ieq
do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
grad_div_mag_q(I,J) = 0.0
enddo ; enddo

Expand Down Expand Up @@ -802,8 +797,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Determine the Laplacian viscosity at h points, using the
! largest value from several parameterizations.
Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity
if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag )
if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3)
if (CS%add_LES_viscosity) then
if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag
if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3
else
if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag )
if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3)
endif
! All viscosity contributions above are subject to resolution scaling
if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j)
Expand All @@ -827,7 +827,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh
if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)
! if (CS%debug) sh_xx_3d(i,j,k) = sh_xx(i,j)

str_xx(i,j) = -Kh * sh_xx(i,j)
else ! not Laplacian
Expand Down Expand Up @@ -965,8 +964,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! Determine the Laplacian viscosity at q points, using the
! largest value from several parameterizations.
Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity
if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag )
if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3)
if (CS%add_LES_viscosity) then
if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag
if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3
else
if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag )
if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3)
endif
! All viscosity contributions above are subject to resolution scaling
if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j)
Expand All @@ -993,7 +997,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh
if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J)
! if (CS%debug) sh_xy_3d(I,J,k) = sh_xy(I,J)

str_xy(I,J) = -Kh * sh_xy(I,J)
else ! not Laplacian
Expand Down Expand Up @@ -1294,8 +1297,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%Laplacian) then
call hchksum(Kh_h, "Kh_h", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T)
! call Bchksum(sh_xy_3d, "shear_xy", G%HI, haloshift=0, scale=US%s_to_T)
! call hchksum(sh_xx_3d, "shear_xx", G%HI, haloshift=0, scale=US%s_to_T)
endif
if (CS%biharmonic) call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
if (CS%biharmonic) call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
Expand Down Expand Up @@ -1504,6 +1505,9 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE)
call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, &
"If true, allow anistropic viscosity in the Laplacian "//&
"horizontal viscosity.", default=.false.)
call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, &
"If true, adds the viscosity from Smagorinsky and Leith to the "//&
"background viscosity instead of taking the maximum.", default=.false.)
endif
if (CS%anisotropic .or. get_all) then
call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, &
Expand Down

0 comments on commit 2cccc83

Please sign in to comment.