Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix bugs in Leith add new input parameter #134

Merged
merged 2 commits into from
Dec 12, 2019
Merged

Fix bugs in Leith add new input parameter #134

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8c95cc0
Fix bugs in Leith add new input parameter
gustavo-marques Dec 10, 2019
78bb4c1
Delete quotes from doxygen comment
gustavo-marques Dec 11, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 26 additions & 22 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Open in desktop
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