Skip to content

Horvisc: Reorder if-blocks in rescaled Kh #862

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

Merged
Hallberg-NOAA merged 1 commit into from
Mar 20, 2025
Merged

Horvisc: Reorder if-blocks in rescaled Kh #862

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
78 changes: 51 additions & 27 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 @@ -1513,8 +1513,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
endif ; endif
endif

meke_res_fn = 1.

if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
do J=js-1,Jeq ; do I=is-1,Ieq
sh_xy_sq = sh_xy(I,J)**2
Expand Down Expand Up @@ -1621,40 +1619,55 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,

! All viscosity contributions above are subject to resolution scaling

! NOTE: The following do-block can be decomposed and vectorized after the
! stack size has been reduced.
do J=js-1,Jeq ; do I=is-1,Ieq
if (rescale_Kh) &
if (rescale_Kh) then
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(I,J) = VarMix%Res_fn_q(I,J) * Kh(I,J)
enddo ; enddo
endif

if (CS%res_scale_MEKE) &
meke_res_fn = VarMix%Res_fn_q(I,J)

if (legacy_bound) then
! Older method of bounding for stability
if (legacy_bound) &
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(I,J) = min(Kh(I,J), CS%Kh_Max_xy(I,J))
enddo ; enddo
endif

do J=js-1,Jeq ; do I=is-1,Ieq
Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired.
enddo ; enddo

if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
if (use_kh_struct) then
do J=js-1,Jeq ; do I=is-1,Ieq
meke_res_fn = 1.
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(I,J)

if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
! *Add* the MEKE contribution (might be negative)
if (use_kh_struct) then
Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
(MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
(MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
else
Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + &
MEKE%Ku(i+1,j+1)) + &
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
meke_res_fn = 1.
if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(I,J)

Kh(I,J) = Kh(I,J) + 0.25 * ( &
(MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + &
(MEKE%Ku(i+1,j) + &
MEKE%Ku(i,j+1)) ) * meke_res_fn
endif
enddo ; enddo
endif
endif

if (CS%anisotropic) &
! *Add* the shear component of anisotropic viscosity
Kh(I,J) = Kh(I,J) + CS%Kh_aniso * CS%n1n2_q(I,J)**2
if (CS%anisotropic) then
! *Add* the shear component of anisotropic viscosity
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(I,J) = Kh(I,J) + CS%Kh_aniso * CS%n1n2_q(I,J)**2
enddo ; enddo
endif

do J=js-1,Jeq ; do I=is-1,Ieq
! Newer method of bounding for stability
if ((CS%better_bound_Kh) .and. (CS%better_bound_Ah)) then
visc_bound_rem(I,J) = 1.0
Expand All @@ -1668,23 +1681,34 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
elseif (CS%better_bound_Kh) then
Kh(I,J) = min(Kh(I,J), hrat_min(I,J) * CS%Kh_Max_xy(I,J))
endif
enddo ; enddo

if (CS%use_Leithy) then
! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points
if (CS%use_Leithy) then
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(I,J) = 0.25 * ((Kh_h(i,j,k) + Kh_h(i+1,j+1,k)) + (Kh_h(i,j+1,k) + Kh_h(i+1,j,k)))
end if
enddo ; enddo
end if

if (CS%id_Kh_q>0 .or. CS%debug) &
if (CS%id_Kh_q > 0 .or. CS%debug) then
do J=js-1,Jeq ; do I=is-1,Ieq
Kh_q(I,J,k) = Kh(I,J)
enddo ; enddo
endif

if (CS%id_vort_xy_q>0) &
if (CS%id_vort_xy_q > 0) then
do J=js-1,Jeq ; do I=is-1,Ieq
vort_xy_q(I,J,k) = vort_xy(I,J)
enddo ; enddo
endif

if (CS%id_sh_xy_q>0) &
if (CS%id_sh_xy_q > 0) then
do J=js-1,Jeq ; do I=is-1,Ieq
sh_xy_q(I,J,k) = sh_xy(I,J)
enddo ; enddo
enddo ; enddo
endif

if ( .not. CS%use_Leithy) then
if (.not. CS%use_Leithy) then
do J=js-1,Jeq ; do I=is-1,Ieq
str_xy(I,J) = -Kh(I,J) * sh_xy(I,J)
enddo ; enddo
Expand Down