diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9dbf7ec3c6..3521131431 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3113,10 +3113,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 @@ -3217,13 +3217,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, & diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 0557ec7cd5..a184cf473f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -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) @@ -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) @@ -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, & @@ -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)), & @@ -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) diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 44a0b0bf5c..fa78938477 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -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)") @@ -837,7 +837,7 @@ 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) @@ -845,7 +845,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, & + 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)") @@ -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, & @@ -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)), & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index f9e4aa0efe..579ddead2d 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -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) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 1c589f509c..65b3bdf50e 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -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) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ab7a6fa5fc..566ef03254 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -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 @@ -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 @@ -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)), & @@ -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)), & @@ -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] @@ -590,6 +597,15 @@ 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, & @@ -597,7 +613,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$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, & !$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, & @@ -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 @@ -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.") diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 4f1dbb89ac..defbd78aa7 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -170,7 +170,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 @@ -474,14 +474,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.") @@ -996,18 +995,47 @@ 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, dz, 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] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer vertical extents [Z ~> m] + 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 @@ -1030,6 +1058,8 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2] real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1] real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1] + real :: Z_to_H ! A local copy of depth to thickness conversion factors or the inverse of the + ! mass-weighted average specific volumes around an interface [H Z-1 ~> nondim or kg m-3] real :: inv_PI3 ! The inverse of pi cubed [nondim] integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -1038,41 +1068,41 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo nz = GV%ke inv_PI3 = 1.0 / ((4.0*atan(1.0))**3) + Z_to_H = GV%Z_to_H ! This will be replaced with a varying value in non-Boussinesq mode. if ((k > 1) .and. (k < nz)) then - ! With USE_QG_LEITH_VISC=True, this might need to change to - ! do j=js-2,je+2 ; do I=is-2,ie+1 - ! but other arrays used here (e.g., h and CS%slope_x) would also need to have wider valid halos. - 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**2 ) + + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**3 ) 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**2 ) - Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) - dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * (GV%Z_to_H * Ih) + + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**3 ) + Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) + if (.not.GV%Boussinesq) & + Z_to_H = ( (h(i,j,k-1) + h(i+1,j,k-1)) + (h(i,j,k) + h(i+1,j,k)) ) / & + ( (dz(i,j,k-1) + dz(i+1,j,k-1)) + (dz(i,j,k) + dz(i+1,j,k)) + GV%dZ_subroundoff) + dslopex_dz(I,j) = 2. * ( slope_x(I,j,k) - slope_x(I,j,k+1) ) * (Z_to_H * Ih) h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - ! With USE_QG_LEITH_VISC=True, this might need to change to - ! do J=js-2,je+1 ; do i=is-2,ie+2 - 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**2 ) + + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**3 ) 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**2 ) - Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) - dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * (GV%Z_to_H * Ih) + + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**3 ) + Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) + if (.not.GV%Boussinesq) & + Z_to_H = ( (h(i,j,k-1) + h(i,j+1,k-1)) + (h(i,j,k) + h(i,j+1,k)) ) / & + ( (dz(i,j,k-1) + dz(i,j+1,k-1)) + (dz(i,j,k) + dz(i,j+1,k)) + GV%dZ_subroundoff) + dslopey_dz(i,J) = 2. * ( slope_y(i,J,k) - slope_y(i,J,k+1) ) * (Z_to_H * Ih) h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - ! With USE_QG_LEITH_VISC=True, this might need to be - ! do J=js-2,je+1 ; do i=is-1,ie+1 - 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 * & ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) & @@ -1080,9 +1110,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo ( ( 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 - ! With USE_QG_LEITH_VISC=True, this might need to be - ! do j=js-1,je+1 ; do I=is-2,ie+1 - 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 * & ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) & @@ -1100,7 +1128,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo + (div_xx_dy(i+1,J) + div_xx_dy(i,J-1))))**2) if (CS%use_beta_in_QG_Leith) then beta_u(I,j) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2)) + (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2)) CS%KH_u_QG(I,j,k) = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), 3.0*beta_u(I,j)) * & CS%Laplac3_const_u(I,j) * inv_PI3 else @@ -1116,7 +1144,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo + (div_xx_dx(I,j+1) + div_xx_dx(I-1,j))))**2) if (CS%use_beta_in_QG_Leith) then beta_v(i,J) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2)) + (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2)) CS%KH_v_QG(i,J,k) = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), 3.0*beta_v(i,J)) * & CS%Laplac3_const_v(i,J) * inv_PI3 else