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

Fake PR for QG Leith Visc bug fix #10

Open
wants to merge 2 commits into
base: dev/gfdl
Choose a base branch
from
Open

Fake PR for QG Leith Visc bug fix #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
727b5c2
(*)Corrected the units of VarMix_CS%slope_x
Hallberg-NOAA Jan 25, 2023
567f7bd
+(*)Fix USE_QG_LEITH_VISC = True reproducibility
Hallberg-NOAA Jan 28, 2023
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
17 changes: 8 additions & 9 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2892,19 +2892,19 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
call cpu_clock_begin(id_clock_pass_init)
call create_group_pass(tmp_pass_uv_T_S_h, CS%u, CS%v, G%Domain)
if (use_temperature) then
call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain, halo=1)
call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain, halo=1)
call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain, halo=2)
call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain, halo=2)
endif
call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain, halo=1)
call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain, halo=2)
call do_group_pass(tmp_pass_uv_T_S_h, G%Domain)
call cpu_clock_end(id_clock_pass_init)

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_m)
call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=2, scale=GV%H_to_m)
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 @@ -2981,11 +2981,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

if (CS%split) then
allocate(eta(SZI_(G),SZJ_(G)), source=0.0)
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, &
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)
if (CS%dtbt_reset_period > 0.0) then
CS%dtbt_reset_interval = real_to_time(CS%dtbt_reset_period)
Expand Down
9 changes: 5 additions & 4 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 @@ -776,7 +776,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 hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_m)
call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, scale=GV%H_to_m)
! 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)
call check_redundant("Predictor uh ", uh, vh, G)
Expand All @@ -785,7 +785,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! 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)
call cpu_clock_end(id_clock_horvisc)
Expand Down Expand Up @@ -1192,7 +1192,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 @@ -1206,6 +1206,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 @@ -1424,7 +1425,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)
call set_initialized(CS%diffu, "diffu", restart_CS)
call set_initialized(CS%diffv, "diffv", restart_CS)
else
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_unsplit.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,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
53 changes: 35 additions & 18 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 @@ -13,7 +13,7 @@ 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_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,7 +22,7 @@ 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_variables, only : accel_diag_ptrs, thermo_var_ptrs

implicit none ; private

Expand Down Expand Up @@ -223,11 +223,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)
CS, tv, dt, OBC, BT, TD, ADp)
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 @@ -245,12 +245,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(in) :: 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(in) :: 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), intent(in), optional :: TD !< Thickness diffusion control structure
type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics

! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: &
Expand Down Expand Up @@ -314,9 +317,11 @@ 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]
Expand Down Expand Up @@ -521,12 +526,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,

endif ! use_GME

if (CS%use_QG_Leith_visc .and. ((CS%Leith_Kh) .or. (CS%Leith_Ah))) then
! 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)
call pass_vector(slope_x, slope_y, G%Domain)
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 apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
!$OMP use_MEKE_Ku, use_MEKE_Au, &
!$OMP use_MEKE_Ku, use_MEKE_Au, slope_x, slope_y, &
!$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 @@ -850,9 +861,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.
! This accumulates terms, some of which are in VarMix.
call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, k, div_xx_dx, div_xx_dy, &
vort_xy_dx, vort_xy_dy)
slope_x, slope_y, vort_xy_dx, vort_xy_dy)

endif

Expand Down Expand Up @@ -1936,8 +1947,14 @@ 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 .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then
call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//&
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.")
endif

Expand Down
Loading