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

Jump to bottom

*+Make Leith viscosity runs layout-invariant #286

Merged
merged 1 commit into from
Jul 16, 2024
Merged

*+Make Leith viscosity runs layout-invariant #286

Show file tree
Hide file tree
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
15 changes: 8 additions & 7 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ module MOM_dynamics_split_RK2
use MOM_debugging, only : check_redundant
use MOM_grid, only : ocean_grid_type
use MOM_hor_index, only : hor_index_type
use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS
use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil
use MOM_hor_visc, only : hor_visc_init, hor_visc_end
use MOM_interface_heights, only : thickness_to_dz, find_col_avg_SpV
use MOM_lateral_mixing_coeffs, only : VarMix_CS
Expand Down Expand Up @@ -401,7 +401,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
logical :: showCallTree, sym

integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: cont_stencil, obc_stencil
integer :: cont_stencil, obc_stencil, vel_stencil

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand Down Expand Up @@ -468,19 +468,20 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
if (associated(CS%OBC)) then
if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3
endif
vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc))
call cpu_clock_begin(id_clock_pass)
call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1)
call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, &
To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil))
call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil))
call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=2)
call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=vel_stencil)
call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil)

call create_group_pass(CS%pass_uv, u, v, G%Domain, halo=max(2,cont_stencil))
call create_group_pass(CS%pass_h, h, G%Domain, halo=max(2,cont_stencil))
call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=vel_stencil)
call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=vel_stencil)
call cpu_clock_end(id_clock_pass)
!--- end set up for group halo pass

Expand Down Expand Up @@ -841,7 +842,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s

if (CS%debug) then
call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=vel_stencil, symmetric=sym, scale=US%L_T_to_m_s)
call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS)
! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s)
Expand Down
75 changes: 51 additions & 24 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 @@ -31,7 +31,7 @@ module MOM_hor_visc

#include <MOM_memory.h>

public horizontal_viscosity, hor_visc_init, hor_visc_end
public horizontal_viscosity, hor_visc_init, hor_visc_end, hor_visc_vel_stencil

!> Control structure for horizontal viscosity
type, public :: hor_visc_CS ; private
Expand Down Expand Up @@ -1198,10 +1198,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
if (CS%Smagorinsky_Ah) then
if (CS%bound_Coriolis) then
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) &
+ CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) &
)
+ CS%Biharm_const2_xx(i,j) * Shear_mag(i,j))
Ah(i,j) = max(Ah(i,j), AhSm)
enddo ; enddo
else
Expand Down Expand Up @@ -1432,10 +1431,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

! Pass the velocity gradients and thickness to ZB2020
if (CS%use_ZB2020) then
call ZB2020_copy_gradient_and_thickness( &
sh_xx, sh_xy, vort_xy, &
hq, &
G, GV, CS%ZB2020, k)
call ZB2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, G, GV, CS%ZB2020, k)
endif

if (CS%Laplacian) then
Expand Down Expand Up @@ -1575,8 +1571,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%bound_Coriolis) then
do J=js-1,Jeq ; do I=is-1,Ieq
AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) &
+ CS%Biharm_const2_xy(I,J) * Shear_mag(I,J) &
)
+ CS%Biharm_const2_xy(I,J) * Shear_mag(I,J))
Ah(I,J) = max(Ah(I,J), AhSm)
enddo ; enddo
else
Expand Down Expand Up @@ -1605,8 +1600,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! *Add* the MEKE contribution
do J=js-1,Jeq ; do I=is-1,Ieq
Ah(I,J) = Ah(I,J) + 0.25 * ( &
(MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) &
)
(MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) )
enddo ; enddo
endif

Expand Down Expand Up @@ -1897,11 +1891,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

if (CS%debug) then
if (CS%Laplacian) then
! In symmetric memory mode, Kh_h should also be valid with a haloshift of 1.
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(Kh_q, "Kh_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**2*US%s_to_T)
endif
if (CS%biharmonic) then
! In symmetric memory mode, Ah_h should also be valid with a haloshift of 1.
call hchksum(Ah_h, "Ah_h", G%HI, haloshift=0, scale=US%L_to_m**4*US%s_to_T)
call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, symmetric=.true., scale=US%L_to_m**4*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)
endif

if (CS%id_FrictWorkIntz > 0) then
Expand Down Expand Up @@ -2403,14 +2401,31 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
ALLOC_(CS%m_leithy_max(isd:ied,jsd:jed)) ; CS%m_leithy_max(:,:) = 0.0
endif
if (CS%Re_Ah > 0.0) then
ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)); CS%Re_Ah_const_xx(:,:) = 0.0
ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)); CS%Re_Ah_const_xy(:,:) = 0.0
ALLOC_(CS%Re_Ah_const_xx(isd:ied,jsd:jed)) ; CS%Re_Ah_const_xx(:,:) = 0.0
ALLOC_(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Re_Ah_const_xy(:,:) = 0.0
endif
endif
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J)
CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
enddo ; enddo

if (((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) .and. &
((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3))) call MOM_error(FATAL, &
"The minimum halo size is 3 when a Leith viscosity is being used.")
if (CS%use_Leithy) then
do J=js-3,Jeq+2 ; do I=is-3,Ieq+2
CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
enddo ; enddo
elseif ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
enddo ; enddo
else
do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
enddo ; enddo
endif

do j=js-2,Jeq+2 ; do i=is-2,Ieq+2
CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j)
CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j)
Expand Down Expand Up @@ -2541,12 +2556,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
endif
endif
if (CS%Leith_Ah) then
CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3)
CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3)
endif
if (CS%use_Leithy) then
CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6
CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j))
CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2
CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6
CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j))
CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2
endif
CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah
Expand All @@ -2571,12 +2586,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
endif
endif
if ((CS%Leith_Ah) .or. (CS%use_Leithy))then
CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3)
CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3)
endif
CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah
if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = &
MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale)
MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale)
if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then
CS%Ah_Max_xy(I,J) = Ah_Limit * (grid_sp_q2 * grid_sp_q2)
CS%Ah_bg_xy(I,J) = MIN(CS%Ah_bg_xy(I,J), CS%Ah_Max_xy(I,J))
Expand Down Expand Up @@ -2822,6 +2837,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)

end subroutine hor_visc_init

!> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size
function hor_visc_vel_stencil(CS) result(stencil)
type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity
integer :: stencil !< The horizontal viscosity velocity stencil size with the current settings.

stencil = 2

if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
stencil = 3
endif
end function hor_visc_vel_stencil

!> Calculates factors in the anisotropic orientation tensor to be align with the grid.
!! With n1=1 and n2=0, this recovers the approach of Large et al, 2001.
subroutine align_aniso_tensor_to_grid(CS, n1, n2)
Expand Down
Loading