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.

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 USE_QG_LEITH_VISC = True reproducibility #604

Merged
merged 2 commits into from
May 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
10 changes: 5 additions & 5 deletions src/core/MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -3114,10 +3114,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &

if (CS%debug) then
call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1)
call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_MKS)
call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=2, scale=GV%H_to_MKS)
if (use_temperature) then
call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=1, scale=US%C_to_degC)
call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=1, scale=US%S_to_ppt)
call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=2, scale=US%C_to_degC)
call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=2, scale=US%S_to_ppt)
endif
endif
endif
Expand Down Expand Up @@ -3221,13 +3221,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
if (CS%split) then
allocate(eta(SZI_(G),SZJ_(G)), source=0.0)
if (CS%use_alt_split) then
call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, &
call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, &
G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, restart_CSp, &
CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, &
CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, &
CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil)
else
call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, &
call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, &
G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, &
CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, &
CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, &
Expand Down
9 changes: 5 additions & 4 deletions src/core/MOM_dynamics_split_RK2.F90
Original file line number Diff line number Diff line change
Expand Up @@ -840,7 +840,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
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 hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS)
call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, 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)
call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
Expand All @@ -849,7 +849,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
! diffu = horizontal viscosity terms (u_av)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
MEKE, Varmix, G, GV, US, CS%hor_visc, &
MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, &
OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
call cpu_clock_end(id_clock_horvisc)
Expand Down Expand Up @@ -1296,7 +1296,7 @@ end subroutine remap_dyn_split_RK2_aux_vars

!> This subroutine initializes all of the variables that are used by this
!! dynamic core, including diagnostics and the cpu clocks.
subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &
diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
VarMix, MEKE, thickness_diffuse_CSp, &
OBC, update_OBC_CSp, ALE_CSp, set_visc, &
Expand All @@ -1310,6 +1310,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
Expand Down Expand Up @@ -1518,7 +1519,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
.not. query_initialized(CS%diffv, "diffv", restart_CS)) then
call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, &
OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
call set_initialized(CS%diffu, "diffu", restart_CS)
call set_initialized(CS%diffv, "diffv", restart_CS)
Expand Down
9 changes: 5 additions & 4 deletions src/core/MOM_dynamics_split_RK2b.F90
Original file line number Diff line number Diff line change
Expand Up @@ -558,7 +558,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
! diffu = horizontal viscosity terms (u_av)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, &
OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%AD_pred)
tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%AD_pred)
call cpu_clock_end(id_clock_horvisc)
if (showCallTree) call callTree_wayPoint("done with predictor horizontal_viscosity (step_MOM_dyn_split_RK2b)")

Expand Down Expand Up @@ -837,15 +837,15 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
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 hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS)
call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, 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)
call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
endif

! diffu = horizontal viscosity terms (u_av)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, &
call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, &
OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%ADp)
call cpu_clock_end(id_clock_horvisc)
if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2b)")
Expand Down Expand Up @@ -1222,7 +1222,7 @@ end subroutine remap_dyn_split_RK2b_aux_vars

!> This subroutine initializes all of the variables that are used by this
!! dynamic core, including diagnostics and the cpu clocks.
subroutine initialize_dyn_split_RK2b(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &
diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
VarMix, MEKE, thickness_diffuse_CSp, &
OBC, update_OBC_CSp, ALE_CSp, set_visc, &
Expand All @@ -1236,6 +1236,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, uh, vh, eta, Time, G, GV, US, para
intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_unsplit.F90
Original file line number Diff line number Diff line change
Expand Up @@ -263,7 +263,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
! diffu = horizontal viscosity terms (u,h)
call enable_averages(dt, Time_local, CS%diag)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc)
call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt)
call cpu_clock_end(id_clock_horvisc)
call disable_averaging(CS%diag)

Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_unsplit_RK2.F90
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
call enable_averages(dt,Time_local, CS%diag)
call cpu_clock_begin(id_clock_horvisc)
call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, &
G, GV, US, CS%hor_visc)
G, GV, US, CS%hor_visc, tv, dt)
call cpu_clock_end(id_clock_horvisc)
call disable_averaging(CS%diag)
call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass)
Expand Down
68 changes: 42 additions & 26 deletions src/parameterizations/lateral/MOM_hor_visc.F90
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@ module MOM_hor_visc
use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_Leith_viscosity
use MOM_interface_heights, only : thickness_to_dz
use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_slopes, calc_QG_Leith_viscosity
use MOM_barotropic, only : barotropic_CS, barotropic_get_tav
use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH
use MOM_io, only : MOM_read_data, slasher
Expand All @@ -22,9 +23,9 @@ module MOM_hor_visc
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE
use MOM_unit_scaling, only : unit_scale_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_variables, only : accel_diag_ptrs
use MOM_Zanna_Bolton, only : ZB2020_lateral_stress, ZB2020_init, ZB2020_end, &
ZB2020_CS, ZB2020_copy_gradient_and_thickness
use MOM_variables, only : accel_diag_ptrs, thermo_var_ptrs
use MOM_Zanna_Bolton, only : ZB2020_lateral_stress, ZB2020_init, ZB2020_end
use MOM_Zanna_Bolton, only : ZB2020_CS, ZB2020_copy_gradient_and_thickness

implicit none ; private

Expand Down Expand Up @@ -237,11 +238,11 @@ module MOM_hor_visc
!!
!! To work, the following fields must be set outside of the usual
!! is:ie range before this subroutine is called:
!! u[is-2:ie+2,js-2:je+2]
!! v[is-2:ie+2,js-2:je+2]
!! h[is-1:ie+1,js-1:je+1]
!! u(is-2:ie+2,js-2:je+2)
!! v(is-2:ie+2,js-2:je+2)
!! h(is-1:ie+1,js-1:je+1) or up to h(is-2:ie+2,js-2:je+2) with some Leith options.
subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
CS, OBC, BT, TD, ADp, hu_cont, hv_cont)
CS, tv, dt, OBC, BT, TD, ADp, hu_cont, hv_cont)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
Expand All @@ -259,12 +260,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
!! related to Mesoscale Eddy Kinetic Energy.
type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure
type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type
type(barotropic_CS), intent(in), optional :: BT !< Barotropic control structure
type(thickness_diffuse_CS), intent(in), optional :: TD !< Thickness diffusion control structure
type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
!! thermodynamic variables
real, intent(in) :: dt !< Time increment [T ~> s]
type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type
type(barotropic_CS), optional, intent(in) :: BT !< Barotropic control structure
type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure
type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
Expand Down Expand Up @@ -341,12 +345,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]
ShSt ! A diagnostic array of shear stress [T-1 ~> s-1].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
KH_u_GME !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
KH_u_GME, & !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
slope_x !< Isopycnal slope in i-direction [Z L-1 ~> nondim]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
KH_v_GME !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
KH_v_GME, & !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
slope_y !< Isopycnal slope in j-direction [Z L-1 ~> nondim]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
Ah_h, & ! biharmonic viscosity at thickness points [L4 T-1 ~> m4 s-1]
Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1]
dz, & ! Height change across layers [Z ~> m]
FrictWork, & ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2]
FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2]
div_xx_h, & ! horizontal divergence [T-1 ~> s-1]
Expand Down Expand Up @@ -590,14 +597,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
call pass_vector(u_smooth, v_smooth, G%Domain)
endif

if (CS%use_QG_Leith_visc .and. ((CS%Leith_Kh) .or. (CS%Leith_Ah))) then
call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=2)
! Calculate isopycnal slopes that will be used for some forms of viscosity.
call calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, VarMix, OBC)
! If the following halo update is added, the calculations in calc_QG_slopes could work on just
! the computational domains, and some halo updates outside of this routine could be smaller.
! call pass_vector(slope_x, slope_y, G%Domain, halo=2)
endif

!$OMP parallel do default(none) &
!$OMP shared( &
!$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, &
!$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
!$OMP is_vort, ie_vort, js_vort, je_vort, &
!$OMP is_Kh, ie_Kh, js_Kh, je_Kh, &
!$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
!$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, &
!$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, &
!$OMP backscat_subround, GME_effic_h, GME_effic_q, &
!$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, &
!$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, &
Expand Down Expand Up @@ -1008,9 +1024,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
(0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2 )
enddo ; enddo

! This accumulates terms, some of which are in VarMix, so rescaling can not be done here.
call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, k, div_xx_dx, div_xx_dy, &
vort_xy_dx, vort_xy_dy)
! This accumulates terms, some of which are in VarMix.
call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, &
slope_x, slope_y, vort_xy_dx, vort_xy_dy)

endif

Expand Down Expand Up @@ -2207,12 +2223,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, &
"If true, use QG Leith nonlinear eddy viscosity.", &
default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) )
if (CS%use_QG_Leith_visc) then
call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//&
"should not be used until a number of bugs are fixed. Specifically it does not "//&
"reproduce across PE count or layout, and may use arrays that have not been properly "//&
"set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.")
endif
! if (CS%use_QG_Leith_visc) then
! call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//&
! "should not be used until a number of bugs are fixed. Specifically it does not "//&
! "reproduce across PE count or layout, and may use arrays that have not been properly "//&
! "set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.")
! endif
if (CS%use_QG_Leith_visc .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then
call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//&
"LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.")
Expand Down
Loading
Loading