diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 07538cdec2..00bc24128d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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 @@ -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) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index e8909e24f9..df4e4a32d5 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -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) @@ -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) @@ -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, & @@ -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)), & @@ -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 diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index e6f99cc9d8..3ad3e25d34 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -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) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index fbf416d13d..861c4bed3c 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -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) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 0574297c0c..e2840ac924 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -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 @@ -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 @@ -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)), & @@ -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)) :: & @@ -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] @@ -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, & @@ -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 @@ -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 diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 85049025b6..d663ef0ae2 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -105,8 +105,8 @@ module MOM_lateral_mixing_coeffs !! spacing squared at v [L2 T-2 ~> m2 s-2]. real, allocatable :: Rd_dx_h(:,:) !< Deformation radius over grid spacing [nondim] - real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [nondim] - real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [nondim] + real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim] + real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim] real, allocatable :: ebt_struct(:,:,:) !< Vertical structure function to scale diffusivities with [nondim] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & @@ -142,7 +142,7 @@ module MOM_lateral_mixing_coeffs integer :: Res_fn_power_visc !< The power of dx/Ld in the Kh resolution function. Any !! positive integer power may be used, but even powers !! and especially 2 are coded to be more efficient. - real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [nondim]. + real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [Z L-1 ~> nondim]. ! Leith parameters logical :: use_QG_Leith_GM !< If true, uses the QG Leith viscosity as the GM coefficient @@ -165,7 +165,7 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function -public calc_QG_Leith_viscosity, calc_depth_function +public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function contains @@ -463,14 +463,13 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure ! Local variables - real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & - e ! The interface heights relative to mean sea level [Z ~> m]. - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& "Module must be initialized before it is used.") @@ -514,28 +513,33 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [nondim] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [Z L-1 ~> nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency !! at u-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope [nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope + !! [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure ! Local variables - real :: S2 ! Interface slope squared [nondim] - real :: N2 ! Positive buoyancy frequency or zero [T-2 ~> s-2] + real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] + real :: N2 ! Positive buoyancy frequency or zero [L2 Z-2 T-2 ~> s-2] real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2] real :: H_geom ! The geometric mean of Hup and Hdn [H ~> m or kg m-2]. - real :: S2max ! An upper bound on the squared slopes [nondim] + real :: S2max ! An upper bound on the squared slopes [Z2 L-2 ~> nondim] real :: wNE, wSE, wSW, wNW ! Weights of adjacent points [nondim] real :: H_u(SZIB_(G)), H_v(SZI_(G)) ! Layer thicknesses at u- and v-points [H ~> m or kg m-2] ! Note that at some points in the code S2_u and S2_v hold the running depth ! integrals of the squared slope [H ~> m or kg m-2] before the average is taken. - real :: S2_u(SZIB_(G),SZJ_(G)) ! The thickness-weighted depth average of the squared slope at u points [nondim]. - real :: S2_v(SZI_(G),SZJB_(G)) ! The thickness-weighted depth average of the squared slope at v points [nondim]. + real :: S2_u(SZIB_(G),SZJ_(G)) ! At first the thickness-weighted depth integral of the squared + ! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the + ! squared slope [Z2 L-2 ~> nondim] at u points. + real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared + ! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the + ! squared slope [Z2 L-2 ~> nondim] at v points. integer :: i, j, k, is, ie, js, je, nz, l_seg @@ -545,7 +549,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C if (.not. CS%calculate_Eady_growth_rate) return if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_u is not associated with use_variable_mixing.") - if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:"// & + if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function:R"// & "%SN_v is not associated with use_variable_mixing.") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -563,7 +567,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C ! and midlatitude deformation radii, using calc_resoln_function as a template. !$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW) - do j = js,je + do j=js,je do I=is-1,ie CS%SN_u(I,j) = 0. ; H_u(I) = 0. ; S2_u(I,j) = 0. enddo @@ -599,7 +603,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C enddo !$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW) - do J = js-1,je + do J=js-1,je do i=is,ie CS%SN_v(i,J) = 0. ; H_v(i) = 0. ; S2_v(i,J) = 0. enddo @@ -644,7 +648,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, & scale=US%Z_to_L, haloshift=1) call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, & - scale=US%L_to_Z**2 * US%s_to_T**2, scalar_pair=.true.) + scale=US%L_to_Z**2*US%s_to_T**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & scale=US%s_to_T, scalar_pair=.true.) endif @@ -698,7 +702,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, enddo ; enddo !$OMP parallel do default(shared) private(dnew,dz,weight,l_seg,vint_SN,sum_dz) - do j = G%jsc-1,G%jec+1 + do j=G%jsc-1,G%jec+1 do I=G%isc-1,G%iec vint_SN(I) = 0. sum_dz(I) = dz_neglect @@ -741,7 +745,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, enddo !$OMP parallel do default(shared) private(dnew,dz,weight,l_seg) - do J = G%jsc-1,G%jec + do J=G%jsc-1,G%jec do i=G%isc-1,G%iec+1 vint_SN(i) = 0. sum_dz(i) = dz_neglect @@ -782,14 +786,14 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, enddo enddo - do j = G%jsc,G%jec + do j=G%jsc,G%jec do I=G%isc-1,G%iec CS%SN_u(I,j) = sqrt( SN_cpy(I,j)**2 & + 0.25*( (CS%SN_v(i,J)**2 + CS%SN_v(i+1,J-1)**2) & + (CS%SN_v(i+1,J)**2 + CS%SN_v(i,J-1)**2) ) ) enddo enddo - do J = G%jsc-1,G%jec + do J=G%jsc-1,G%jec do i=G%isc,G%iec CS%SN_v(i,J) = sqrt( CS%SN_v(i,J)**2 & + 0.25*( (SN_cpy(I,j)**2 + SN_cpy(I-1,j+1)**2) & @@ -816,13 +820,13 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop logical, intent(in) :: calculate_slopes !< If true, calculate slopes !! internally otherwise use slopes stored in CS ! Local variables - real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics) - real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics) + real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics) + real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics) real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. - real :: S2 ! Interface slope squared [nondim] - real :: N2 ! Brunt-Vaisala frequency squared [T-2 ~> s-2] + real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim] + real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2] real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2] real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2]. real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times @@ -857,16 +861,16 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop if (calculate_slopes) then ! Calculate the interface slopes E_x and E_y and u- and v- points respectively do j=js-1,je+1 ; do I=is-1,ie - E_x(I,j) = US%Z_to_L*(e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j) + E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j) ! Mask slopes where interface intersects topography if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0. enddo ; enddo do J=js-1,je ; do i=is-1,ie+1 - E_y(i,J) = US%Z_to_L*(e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J) + E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J) ! Mask slopes where interface intersects topography if (min(h(i,J,k),h(i,J+1,k)) < H_cutoff) E_y(I,j) = 0. enddo ; enddo - else + else ! This branch is not used. do j=js-1,je+1 ; do I=is-1,ie E_x(I,j) = CS%slope_x(I,j,k) if (min(h(I,j,k),h(I+1,j,k)) < H_cutoff) E_x(I,j) = 0. @@ -880,22 +884,22 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop ! Calculate N*S*h from this layer and add to the sum do j=js,je ; do I=is-1,ie S2 = ( E_x(I,j)**2 + 0.25*( & - (E_y(I,j)**2+E_y(I+1,j-1)**2)+(E_y(I+1,j)**2+E_y(I,j-1)**2) ) ) + (E_y(I,j)**2+E_y(I+1,j-1)**2) + (E_y(I+1,j)**2+E_y(I,j-1)**2) ) ) Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect) Hup = 2.*h(i+1,j,k)*h(i+1,j,k-1) / (h(i+1,j,k) + h(i+1,j,k-1) + h_neglect) H_geom = sqrt(Hdn*Hup) - N2 = GV%g_prime(k)*US%L_to_Z**2 / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) + N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < H_cutoff) & S2 = 0.0 S2N2_u_local(I,j,k) = (H_geom * GV%H_to_Z) * S2 * N2 enddo ; enddo do J=js-1,je ; do i=is,ie S2 = ( E_y(i,J)**2 + 0.25*( & - (E_x(i,J)**2+E_x(i-1,J+1)**2)+(E_x(i,J+1)**2+E_x(i-1,J)**2) ) ) + (E_x(i,J)**2+E_x(i-1,J+1)**2) + (E_x(i,J+1)**2+E_x(i-1,J)**2) ) ) Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect) Hup = 2.*h(i,j+1,k)*h(i,j+1,k-1) / (h(i,j+1,k) + h(i,j+1,k-1) + h_neglect) H_geom = sqrt(Hdn*Hup) - N2 = GV%g_prime(k)*US%L_to_Z**2 / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) + N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2)) if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < H_cutoff) & S2 = 0.0 S2N2_v_local(i,J,k) = (H_geom * GV%H_to_Z) * S2 * N2 @@ -942,32 +946,60 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop end subroutine calc_slope_functions_using_just_e + +!> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity. +subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + real, intent(in) :: dt !< Time increment [T ~> s] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] + type(VarMix_CS), intent(in) :: CS !< Variable mixing control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] + + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_QG_slopes: "//& + "Module must be initialized before it is used.") + + call find_eta(h, tv, G, GV, US, e, halo_size=3) + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & + slope_x, slope_y, halo=2, OBC=OBC) + +end subroutine calc_QG_slopes + !> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients -subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy) +subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, slope_x, slope_y, & + vort_xy_dx, vort_xy_dy) type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] - integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude - real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence !! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence !! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity !! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] ! Local variables real, dimension(SZI_(G),SZJB_(G)) :: & - dslopey_dz, & ! z-derivative of y-slope at v-points [Z-1 ~> m-1] + dslopey_dz, & ! z-derivative of y-slope at v-points [L-1 ~> m-1] h_at_v, & ! Thickness at v-points [H ~> m or kg m-2] beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1] grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1] grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1] real, dimension(SZIB_(G),SZJ_(G)) :: & - dslopex_dz, & ! z-derivative of x-slope at u-points [Z-1 ~> m-1] + dslopex_dz, & ! z-derivative of x-slope at u-points [L-1 ~> m-1] h_at_u, & ! Thickness at u-points [H ~> m or kg m-2] beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1] grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1] @@ -987,41 +1019,41 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo if ((k > 1) .and. (k < nz)) then - do j=js-1,je+1 ; do I=is-2,Ieq+1 + do j=js-2,je+2 ; do I=is-2,ie+1 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / & ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) & - + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff ) + + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**2 ) h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / & ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) & - + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff ) + + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**2 ) Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z ) - dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * Ih + dslopex_dz(I,j) = 2. * ( slope_x(I,j,k) - slope_x(I,j,k+1) ) * Ih h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - do J=js-2,Jeq+1 ; do i=is-1,ie+1 + do J=js-2,je+1 ; do i=is-2,ie+2 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / & ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) & - + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff ) + + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**2 ) h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / & ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) & - + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff ) + + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**2 ) Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z ) - dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * Ih + dslopey_dz(i,J) = 2. * ( slope_y(i,J,k) - slope_y(i,J,k+1) ) * Ih h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - do J=js-1,je ; do i=is-1,Ieq+1 + do J=js-2,je+1 ; do i=is-1,ie+1 f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) ) - vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * US%L_to_Z * & + vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * & ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) & + ( h_at_u(I-1,j) * dslopex_dz(I-1,j) + h_at_u(I,j+1) * dslopex_dz(I,j+1) ) ) / & ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff) enddo ; enddo - do j=js-1,Jeq+1 ; do I=is-1,ie + do j=js-1,je+1 ; do I=is-2,ie+1 f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) ) - vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * US%L_to_Z * & + vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * & ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) & + ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / & ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) + GV%H_subroundoff) @@ -1236,7 +1268,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If non-zero, is an upper bound on slopes used in the "//& "Visbeck formula for diffusivity. This does not affect the "//& "isopycnal slope calculation used within thickness diffusion.", & - units="nondim", default=0.0) + units="nondim", default=0.0, scale=US%L_to_Z) else CS%Visbeck_S_max = 0. endif