From e27eec5104b272c310f8a4ea5a60f12bec1d5583 Mon Sep 17 00:00:00 2001 From: Elizabeth Yankovsky Date: Wed, 7 Aug 2024 14:28:43 -0400 Subject: [PATCH 1/5] EBT Backscatter Backscatter code using the equivalent barotropic (EBT) mode documented in Yankovsky et al. (2024). --- src/parameterizations/lateral/MOM_MEKE.F90 | 88 ++++- .../lateral/MOM_MEKE_types.F90 | 2 + .../lateral/MOM_hor_visc.F90 | 346 ++++++++++++++++-- .../lateral/MOM_lateral_mixing_coeffs.F90 | 26 +- 4 files changed, 416 insertions(+), 46 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index c3859c3bd1..c2ba54c659 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -55,6 +55,7 @@ module MOM_MEKE logical :: initialized = .false. !< True if this control structure has been initialized. ! Parameters real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim] + real :: MEKE_bhFrCoeff!< Efficiency of conversion of ME into MEKE by the biharmonic dissipation [nondim] real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim] real :: MEKE_GMECoeff !< Efficiency of conversion of MEKE into ME by GME [nondim] real :: MEKE_damping !< Local depth-independent MEKE dissipation rate [T-1 ~> s-1]. @@ -126,8 +127,10 @@ module MOM_MEKE type(diag_ctrl), pointer :: diag => NULL() !< A type that regulates diagnostics output !>@{ Diagnostic handles integer :: id_MEKE = -1, id_Ue = -1, id_Kh = -1, id_src = -1 + integer :: id_src_adv = -1, id_src_mom_K4 = -1, id_src_btm_drag = -1 + integer :: id_src_GM = -1, id_src_mom_lp = -1, id_src_mom_bh = -1 integer :: id_Ub = -1, id_Ut = -1 - integer :: id_GM_src = -1, id_mom_src = -1, id_GME_snk = -1, id_decay = -1 + integer :: id_GM_src = -1, id_mom_src = -1, id_mom_src_bh = -1, id_GME_snk = -1, id_decay = -1 integer :: id_KhMEKE_u = -1, id_KhMEKE_v = -1, id_Ku = -1, id_Au = -1 integer :: id_Le = -1, id_gamma_b = -1, id_gamma_t = -1 integer :: id_Lrhines = -1, id_Leady = -1 @@ -192,6 +195,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h depth_tot, & ! The depth of the water column [H ~> m or kg m-2]. src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3). MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1]. + src_adv, & ! The MEKE source/tendency from the horizontal advection of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_K4, & ! The MEKE source/tendency from the bihamornic of MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_btm_drag, & ! The MEKE source/tendency from the bottom drag acting on MEKE [L2 T-3 ~> W kg-1] (= m2 s-3). + src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). + src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). + ldamping_Strang1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1]. + MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1] drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1]. del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2]. @@ -225,6 +236,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). + real :: sfac ! A factor needed to compute damping due to Strang splitting [nondim]c logical :: use_drag_rate ! Flag to indicate drag_rate is finite integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features @@ -254,6 +266,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%debug) then if (allocated(MEKE%mom_src)) & call hchksum(MEKE%mom_src, 'MEKE mom_src', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + if (allocated(MEKE%mom_src_bh)) & + call hchksum(MEKE%mom_src_bh, 'MEKE mom_src_bh', G%HI, scale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GME_snk)) & call hchksum(MEKE%GME_snk, 'MEKE GME_snk', G%HI, unscale=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (allocated(MEKE%GM_src)) & @@ -387,12 +401,22 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = CS%MEKE_BGsrc + src_adv(i,j) = 0. + src_mom_K4(i,j) = 0. + src_btm_drag(i,j) = 0. + src_GM(i,j) = 0. + src_mom_lp(i,j) = 0. + src_mom_bh(i,j) = 0. enddo ; enddo if (allocated(MEKE%mom_src)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) + !src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) + src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) & + - (CS%MEKE_bhFrCoeff-CS%MEKE_FrCoeff)*I_mass(i,j)*MEKE%mom_src_bh(i,j) + src_mom_lp(i,j) = - CS%MEKE_FrCoeff*I_mass(i,j)*(MEKE%mom_src(i,j)-MEKE%mom_src_bh(i,j)) + src_mom_bh(i,j) = - CS%MEKE_bhFrCoeff*I_mass(i,j)*MEKE%mom_src_bh(i,j) enddo ; enddo endif @@ -414,6 +438,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) + src_GM(i,j) = -CS%MEKE_GMcoeff*I_mass(i,j)*MEKE%GM_src(i,j) enddo ; enddo endif endif @@ -433,6 +458,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie + MEKE_current(i,j) = MEKE%MEKE(i,j) MEKE%MEKE(i,j) = (MEKE%MEKE(i,j) + sdt*src(i,j))*G%mask2dT(i,j) enddo ; enddo @@ -459,6 +485,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! while leaving MEKE unchanged if it is negative MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + ldamping_Strang1(i,j) = ldamping + src_GM(i,j) = src_GM(i,j) / (1.0 + sdt_damp*ldamping) + src_mom_lp(i,j) = src_mom_lp(i,j) / (1.0 + sdt_damp*ldamping) + src_mom_bh(i,j) = src_mom_bh(i,j) / (1.0 + sdt_damp*ldamping) + sfac = ( 1.0 + sdt_damp*ldamping ) + src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) enddo ; enddo if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then @@ -528,6 +560,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h del4MEKE(i,j) = (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) + src_mom_K4(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & + ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & + (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) enddo ; enddo endif ! @@ -595,6 +630,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h MEKE%MEKE(i,j) = MEKE%MEKE(i,j) + (sdt*(G%IareaT(i,j)*I_mass(i,j))) * & ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) + src_adv(i,j) = (G%IareaT(i,j)*I_mass(i,j)) * & + ((MEKE_uflux(I-1,j) - MEKE_uflux(I,j)) + & + (MEKE_vflux(i,J-1) - MEKE_vflux(i,J))) enddo ; enddo endif ! MEKE_KH>0 @@ -625,6 +663,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! while leaving MEKE unchanged if it is negative MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + src_GM(i,j) = src_GM(i,j) / (1.0 + sdt_damp*ldamping) + src_mom_lp(i,j) = src_mom_lp(i,j) / (1.0 + sdt_damp*ldamping) + src_mom_bh(i,j) = src_mom_bh(i,j) / (1.0 + sdt_damp*ldamping) + src_adv(i,j) = src_adv(i,j) / (1.0 + sdt_damp*ldamping) + src_mom_K4(i,j) = src_mom_K4(i,j) / (1.0 + sdt_damp*ldamping) + sfac = ( 1.0 + sdt_damp*ldamping_Strang1(i,j) ) * ( 1.0 + sdt_damp*ldamping ) + src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) enddo ; enddo endif endif ! MEKE_KH>=0 @@ -727,9 +772,16 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (CS%id_KhMEKE_u>0) call post_data(CS%id_KhMEKE_u, Kh_u, CS%diag) if (CS%id_KhMEKE_v>0) call post_data(CS%id_KhMEKE_v, Kh_v, CS%diag) if (CS%id_src>0) call post_data(CS%id_src, src, CS%diag) + if (CS%id_src_adv>0) call post_data(CS%id_src_adv, src_adv, CS%diag) + if (CS%id_src_mom_K4>0) call post_data(CS%id_src_mom_K4, src_mom_K4, CS%diag) + if (CS%id_src_btm_drag>0) call post_data(CS%id_src_btm_drag, src_btm_drag, CS%diag) + if (CS%id_src_GM>0) call post_data(CS%id_src_GM, src_GM, CS%diag) + if (CS%id_src_mom_lp>0) call post_data(CS%id_src_mom_lp, src_mom_lp, CS%diag) + if (CS%id_src_mom_bh>0) call post_data(CS%id_src_mom_bh, src_mom_bh, CS%diag) if (CS%id_decay>0) call post_data(CS%id_decay, MEKE_decay, CS%diag) if (CS%id_GM_src>0) call post_data(CS%id_GM_src, MEKE%GM_src, CS%diag) if (CS%id_mom_src>0) call post_data(CS%id_mom_src, MEKE%mom_src, CS%diag) + if (CS%id_mom_src_bh>0) call post_data(CS%id_mom_src_bh, MEKE%mom_src_bh, CS%diag) if (CS%id_GME_snk>0) call post_data(CS%id_GME_snk, MEKE%GME_snk, CS%diag) if (CS%id_Le>0) call post_data(CS%id_Le, LmixScale, CS%diag) if (CS%id_gamma_b>0) then @@ -1210,6 +1262,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME "The efficiency of the conversion of mean energy into "//& "MEKE. If MEKE_FRCOEFF is negative, this conversion "//& "is not used or calculated.", units="nondim", default=-1.0) + call get_param(param_file, mdl, "MEKE_BHFRCOEFF", CS%MEKE_bhFrCoeff, & + "The efficiency of the conversion of mean energy into "//& + "MEKE by the biharmonic dissipation. If MEKE_bhFRCOEFF is negative, this conversion "//& + "is not used or calculated.", units="nondim", default=-1.0) call get_param(param_file, mdl, "MEKE_GMECOEFF", CS%MEKE_GMECoeff, & "The efficiency of the conversion of MEKE into mean energy "//& "by GME. If MEKE_GMECOEFF is negative, this conversion "//& @@ -1399,6 +1455,20 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (.not. allocated(MEKE%MEKE)) CS%id_Ut = -1 CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, & 'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + !add diagnostics for the terms in the MEKE budget + CS%id_src_adv = register_diag_field('ocean_model', 'MEKE_src_adv', diag%axesT1, Time, & + 'MEKE energy source from the horizontal advection of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_K4 = register_diag_field('ocean_model', 'MEKE_src_mom_K4', diag%axesT1, Time, & + 'MEKE energy source from the biharmonic of MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_btm_drag = register_diag_field('ocean_model', 'MEKE_src_btm_drag', diag%axesT1, Time, & + 'MEKE energy source from the bottom drag acting on MEKE', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_GM = register_diag_field('ocean_model', 'MEKE_src_GM', diag%axesT1, Time, & + 'MEKE energy source from the thickness mixing (GM scheme)', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_lp = register_diag_field('ocean_model', 'MEKE_src_mom_lp', diag%axesT1, Time, & + 'MEKE energy source from the Laplacian of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + CS%id_src_mom_bh = register_diag_field('ocean_model', 'MEKE_src_mom_bh', diag%axesT1, Time, & + 'MEKE energy source from the biharmonic of resolved flows', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) + !end CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, & 'MEKE decay rate', 's-1', conversion=US%s_to_T) CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, & @@ -1409,6 +1479,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME 'MEKE energy available from momentum', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) if (.not. allocated(MEKE%mom_src)) CS%id_mom_src = -1 + CS%id_mom_src_bh = register_diag_field('ocean_model', 'MEKE_mom_src_bh',diag%axesT1, Time, & + 'MEKE energy available from the biharmonic dissipation of momentum', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) + if (.not. allocated(MEKE%mom_src_bh)) CS%id_mom_src_bh = -1 CS%id_GME_snk = register_diag_field('ocean_model', 'MEKE_GME_snk',diag%axesT1, Time, & 'MEKE energy lost to GME backscatter', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) @@ -1742,7 +1816,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct ! Local variables - real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim] + real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_bhFrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim] real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim] logical :: Use_KH_in_MEKE logical :: useMEKE @@ -1754,6 +1828,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) ! Read these parameters to determine what should be in the restarts MEKE_GMcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff) MEKE_FrCoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff) + MEKE_bhFrCoeff = -1. ; call read_param(param_file,"MEKE_bhFRCOEFF",MEKE_bhFrCoeff) MEKE_GMEcoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff) MEKE_KhCoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff) MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) @@ -1770,8 +1845,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2) if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0) - if (MEKE_FrCoeff>=0. .or. MEKE_GMECoeff>=0.) & + if (MEKE_FrCoeff>=0. .or. MEKE_bhFrCoeff>=0. .or. MEKE_GMECoeff>=0.) then allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0) + allocate(MEKE%mom_src_bh(isd:ied,jsd:jed), source=0.0) + endif + if (MEKE_FrCoeff<0.) MEKE_FrCoeff = 0. + if (MEKE_bhFrCoeff<0.) MEKE_bhFrCoeff = 0. if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0) if (MEKE_KhCoeff>=0.) then allocate(MEKE%Kh(isd:ied,jsd:jed), source=0.0) @@ -1817,6 +1896,7 @@ subroutine MEKE_end(MEKE) if (allocated(MEKE%Kh)) deallocate(MEKE%Kh) if (allocated(MEKE%GME_snk)) deallocate(MEKE%GME_snk) if (allocated(MEKE%mom_src)) deallocate(MEKE%mom_src) + if (allocated(MEKE%mom_src_bh)) deallocate(MEKE%mom_src_bh) if (allocated(MEKE%GM_src)) deallocate(MEKE%GM_src) if (allocated(MEKE%MEKE)) deallocate(MEKE%MEKE) end subroutine MEKE_end diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index e51f558ce3..890ee58824 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -11,6 +11,8 @@ module MOM_MEKE_types real, allocatable :: GM_src(:,:) !< MEKE source due to thickness mixing (GM) [R Z L2 T-3 ~> W m-2]. real, allocatable :: mom_src(:,:) !< MEKE source from lateral friction in the !! momentum equations [R Z L2 T-3 ~> W m-2]. + real, allocatable :: mom_src_bh(:,:) !< MEKE source from the biharmonic part of the lateral friction in the + !! momentum equations [R Z L2 T-3 ~> W m-2]. !cyc real, allocatable :: GME_snk(:,:) !< MEKE sink from GME backscatter in the momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: Kh(:,:) !< The MEKE-derived lateral mixing coefficient [L2 T-1 ~> m2 s-1]. real, allocatable :: Kh_diff(:,:) !< Uses the non-MEKE-derived thickness diffusion coefficient to diffuse diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 94afd0c858..0be1674b99 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -50,6 +50,9 @@ module MOM_hor_visc !! limited to guarantee stability. logical :: better_bound_Kh !< If true, use a more careful bounding of the !! Laplacian viscosity to guarantee stability. + logical :: EY24_EBT_BS !! If true, use an equivalent barotropic backscatter + !! with a stabilizing kill switch in MEKE, + !< developed by Yankovsky et al. 2024 logical :: bound_Ah !< If true, the biharmonic coefficient is locally !! limited to guarantee stability. logical :: better_bound_Ah !< If true, use a more careful bounding of the @@ -60,6 +63,9 @@ module MOM_hor_visc !! the viscosity bounds to the theoretical maximum !! for stability without considering other terms [nondim]. !! The default is 0.8. + real :: KS_coef !< A nondimensional coefficient on the biharmonic viscosity that sets the + !! kill switch for backscatter. Default is 1.0. + real :: KS_timescale !< A timescale for computing CFL limit for turning off backscatter (~DT) logical :: backscatter_underbound !< If true, the bounds on the biharmonic viscosity are allowed !! to increase where the Laplacian viscosity is negative (due to !! backscatter parameterizations) beyond the largest timestep-dependent @@ -144,6 +150,7 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & Kh_Max_xx, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1]. Ah_Max_xx, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1]. + Ah_Max_xx_KS, & !< The maximum permitted biharmonic viscosity for kill switch [L4 T-1 ~> m4 s-1]. n1n2_h, & !< Factor n1*n2 in the anisotropic direction tensor at h-points [nondim] n1n1_m_n2n2_h, & !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points [nondim] grid_sp_h2, & !< Harmonic mean of the squares of the grid [L2 ~> m2] @@ -162,6 +169,7 @@ module MOM_hor_visc real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & Kh_Max_xy, & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1]. Ah_Max_xy, & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1]. + Ah_Max_xy_KS, & !< The maximum permitted biharmonic viscosity for kill switch [L4 T-1 ~> m4 s-1]. n1n2_q, & !< Factor n1*n2 in the anisotropic direction tensor at q-points [nondim] n1n1_m_n2n2_q !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points [nondim] @@ -229,8 +237,13 @@ module MOM_hor_visc integer :: id_vort_xy_q = -1, id_div_xx_h = -1 integer :: id_sh_xy_q = -1, id_sh_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 + integer :: id_FrictWork_bh = -1, id_FrictWorkIntz_bh = -1 !cyc integer :: id_FrictWork_GME = -1 integer :: id_normstress = -1, id_shearstress = -1 + integer :: id_visc_limit_h = -1, id_visc_limit_q = -1 + integer :: id_visc_limit_h_flag = -1, id_visc_limit_q_flag = -1 + integer :: id_visc_limit_h_frac = -1, id_visc_limit_q_frac = -1 + integer :: id_BS_coeff_h = -1, id_BS_coeff_q = -1 !>@} end type hor_visc_CS @@ -312,6 +325,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [L2 T-2 ~> m2 s-2] bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] + FrictWorkIntz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R L2 T-3 ~> W m-2] !cyc grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] @@ -320,7 +334,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] m_leithy, & ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2] Ah_sq, & ! The square of the biharmonic viscosity [L8 T-2 ~> m8 s-2] - htot ! The total thickness of all layers [H ~> m or kg m-2] + htot, & ! The total thickness of all layers [H ~> m or kg m-2] + str_xx_BS ! The diagonal term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2] real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] real :: boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] @@ -345,7 +360,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] ! This form guarantees that hq/hu < 4. - GME_effic_q ! The filtered efficiency of the GME terms at q points [nondim] + GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim] + str_xy_BS ! The cross term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2] real :: grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] real :: boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] @@ -355,6 +371,10 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1] sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1] GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1] + visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter + visc_limit_q_flag, & ! determines whether backscatter is shut off + visc_limit_q_frac, & ! determines how close backscatter is to shutting off + BS_coeff_q, & ! A diagnostic array of the backscatter coefficient [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] @@ -367,14 +387,19 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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_bh, & ! work done by the biharmonic MKE dissipation mechanisms [R L2 T-3 ~> W m-2] !cyc FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2] div_xx_h, & ! horizontal divergence [T-1 ~> s-1] sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] - NoSt ! A diagnostic array of normal stress [T-1 ~> s-1]. + NoSt, & ! A diagnostic array of normal stress [T-1 ~> s-1]. + BS_coeff_h ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] - GME_coeff_h ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] + GME_coeff_h, & ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] + visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter + visc_limit_h_flag, & ! determines whether backscatter is shut off + visc_limit_h_frac ! determines how close backscatter is to shutting off real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & @@ -425,6 +450,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim] + real :: tmp ! Fields evaluated on active layers, used for constructing 3D stress fields ! NOTE: The position of these declarations can impact performance, due to the @@ -435,6 +461,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, real, dimension(SZIB_(G),SZJB_(G)) :: & Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1] Kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1] + Kh_BS, & ! Laplacian antiviscosity [L2 T-1 ~> m2 s-1] Shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1] vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1] vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (h or q) [L-1 T-1 ~> m-1 s-1] @@ -451,6 +478,13 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 + visc_limit_h(:,:,:) = 0. + visc_limit_q(:,:,:) = 0. + visc_limit_h_flag(:,:,:) = 0. + visc_limit_q_flag(:,:,:) = 0. + visc_limit_h_frac(:,:,:) = 0. + visc_limit_q_frac(:,:,:) = 0. + m_leithy(:,:) = 0.0 ! Initialize if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then @@ -629,12 +663,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$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, & + !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_bh, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont & !$OMP ) & !$OMP private( & - !$OMP i, j, k, n, & + !$OMP i, j, k, n, tmp, & !$OMP dudx, dudy, dvdx, dvdy, sh_xx, sh_xy, h_u, h_v, & !$OMP Del2u, Del2v, DY_dxBu, DX_dyBu, sh_xx_bt, sh_xy_bt, & !$OMP str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, & @@ -650,7 +684,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, & !$OMP vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, & !$OMP sh_xx_smooth, sh_xy_smooth, & - !$OMP vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy & + !$OMP vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy, & + !$OMP Kh_BS, str_xx_bs, str_xy_bs, bs_coeff_h, bs_coeff_q & + !$OMP ) & + !$OMP firstprivate( & + !$OMP visc_limit_h, visc_limit_h_frac, visc_limit_h_flag, & + !$OMP visc_limit_q, visc_limit_q_frac, visc_limit_q_flag & !$OMP ) do k=1,nz @@ -1139,15 +1178,15 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) enddo ; enddo - if (use_MEKE_Ku) then + if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then ! *Add* the MEKE contribution (which might be negative) if (CS%res_scale_MEKE) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k) enddo ; enddo else do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) enddo ; enddo endif endif @@ -1349,6 +1388,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif + if (CS%EY24_EBT_BS) then + do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh + tmp = CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j) + visc_limit_h(i,j,k) = tmp + visc_limit_h_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j)) + if (Ah(i,j) >= tmp) then + visc_limit_h_flag(i,j,k) = 1. + endif + enddo ; enddo + endif + if ((CS%id_Ah_h>0) .or. CS%debug .or. CS%use_Leithy) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh Ah_h(i,j,k) = Ah(i,j) @@ -1364,7 +1414,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (CS%id_grid_Re_Ah>0) then + if ((CS%id_grid_Re_Ah>0)) then do j=js,je ; do i=is,ie KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) grid_Ah = max(Ah(i,j), CS%min_grid_Ah) @@ -1386,6 +1436,31 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx + ! Backscatter using MEKE + if (CS%EY24_EBT_BS) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_h_flag(i,j,k) > 0) then + Kh_BS(i,j) = 0. + else + Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) + endif + enddo ; enddo + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx_BS(i,j) = -Kh_BS(i,j) * sh_xx(i,j) + enddo ; enddo + + if (CS%id_BS_coeff_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + BS_coeff_h(i,j,k) = Kh_BS(i,j) + enddo ; enddo + endif + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx(i,j) = str_xx(i,j) + str_xx_BS(i,j) + enddo ; enddo + endif ! Backscatter + if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term do J=js-1,Jeq ; do I=is-1,Ieq @@ -1540,10 +1615,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired. - if (use_MEKE_Ku) then + if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then ! *Add* the MEKE contribution (might be negative) - Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & - (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn + Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn endif if (CS%anisotropic) & @@ -1670,6 +1747,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif + if (CS%EY24_EBT_BS) then + do J=js-1,Jeq ; do I=is-1,Ieq + tmp = CS%KS_coef *hrat_min(I,J) * CS%Ah_Max_xy_KS(I,J) + visc_limit_q(I,J,k) = tmp + visc_limit_q_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xy_KS(i,j)) + if (Ah(I,J) >= tmp) then + visc_limit_q_flag(I,J,k) = 1. + endif + enddo ; enddo + endif + ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points if (CS%use_Leithy) then do J=js-1,Jeq ; do I=is-1,Ieq @@ -1694,6 +1782,35 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif ! Get Ah at q points and biharmonic part of str_xy + ! Backscatter using MEKE + if (CS%EY24_EBT_BS) then + !if (grid_Re_Ah(i,j,k) < 1.) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_q_flag(I,J,k) > 0) then + Kh_BS(I,J) = 0. + else + Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) + endif + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy_BS(I,J) = -Kh_BS(I,J) * (sh_xy(I,J)) + enddo ; enddo + + if (CS%id_BS_coeff_q>0) then + do J=js-1,Jeq ; do I=is-1,Ieq + BS_coeff_q(I,J,k) = Kh_BS(I,J) + enddo ; enddo + endif + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = str_xy(I,J) + str_xy_BS(I,J) + enddo ; enddo + endif ! Backscatter + if (CS%use_GME) then ! The wider halo here is to permit one pass of smoothing without a halo update. do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 @@ -1797,31 +1914,57 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) ! This is the old formulation that includes energy diffusion - FrictWork(i,j,k) = GV%H_to_RZ * ( & - (str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - - str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - + 0.25*((str_xy(I,J) * & - ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & - + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & - + str_xy(I-1,J-1) * & - ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & - + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & - + (str_xy(I-1,J) * & - ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & - + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & - + str_xy(I,J-1) * & - ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & - + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) - enddo ; enddo + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork(i,j,k) = 0 + FrictWork_bh(i,j,k) = 0 + else + FrictWork(i,j,k) = GV%H_to_RZ * ( & + (str_xx(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + -str_xx(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + +0.25*((str_xy(I,J)*( & + (u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + +(v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) & + +str_xy(I-1,J-1)*( & + (u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + +(v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1) )) & + +(str_xy(I-1,J)*( & + (u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + +(v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J) ) & + +str_xy(I,J-1)*( & + (u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + +(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) ) + ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion !cyc + FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & + (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + + 0.25*((bhstr_xy(I,J) * & + ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & + + bhstr_xy(I-1,J-1) * & + ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & + + (bhstr_xy(I-1,J) * & + ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & + + bhstr_xy(I,J-1) * & + ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + endif + enddo ; enddo else ; do j=js,je ; do i=is,ie - FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork(i,j,k) = 0 + FrictWork_bh(i,j,k) = 0 + else + FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((str_xx(i,j)*CS%dy2h(i,j) * ( & (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) & - (str_xx(i,j)*CS%dx2h(i,j) * ( & (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) & - + (0.25*(((str_xy(I,J)*( & + + (0.25*(((str_xy(I,J)*( & (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) & - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) & + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) & @@ -1841,10 +1984,42 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) + + ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) + FrictWork_bh(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & + ((bhstr_xx(i,j)*CS%dy2h(i,j) * ( & + (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & + - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) & + - (bhstr_xx(i,j)*CS%dx2h(i,j) * ( & + (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & + - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) & + + (0.25*(((bhstr_xy(I,J)*( & + (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) & + - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) & + + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) & + - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) & + +(bhstr_xy(I-1,J-1)*( & + (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) & + - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) & + + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) & + - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) & + +((bhstr_xy(I-1,J)*( & + (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) & + - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) & + + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) & + - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) & + +(bhstr_xy(I,J-1)*( & + (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) & + - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) & + + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & + - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) + endif enddo ; enddo ; endif endif + + if (CS%use_GME) then if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) @@ -1904,6 +2079,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (k==1) then do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = 0. + MEKE%mom_src_bh(i,j) = 0. !cyc enddo ; enddo if (allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie @@ -1951,12 +2127,29 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, + (str_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + GV%H_to_RZ * ( & !cyc + ((bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + -(bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + + 0.25*(((bhstr_xy(I,J)-RoScl*bhstr_xy(I,J)) * & + ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) & + + (bhstr_xy(I-1,J-1)-RoScl*bhstr_xy(I-1,J-1)) * & + ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & + + ((bhstr_xy(I-1,J)-RoScl*bhstr_xy(I-1,J)) * & + ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & + + (bhstr_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & + ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) enddo ; enddo - endif ! MEKE%backscatter_Ro_c + else !cyc do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) !cyc enddo ; enddo + endif ! MEKE%backscatter_Ro_c if (CS%use_GME .and. allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie @@ -1974,6 +2167,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag) if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag) if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag) + if (CS%id_FrictWork_bh>0) call post_data(CS%id_FrictWork_bh, FrictWork_bh, CS%diag) !cyc if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) @@ -1993,6 +2187,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_dudy_bt > 0) call post_data(CS%id_dudy_bt, dudy_bt, CS%diag) if (CS%id_dvdx_bt > 0) call post_data(CS%id_dvdx_bt, dvdx_bt, CS%diag) endif + if (CS%id_visc_limit_h>0) call post_data(CS%id_visc_limit_h, visc_limit_h, CS%diag) + if (CS%id_visc_limit_q>0) call post_data(CS%id_visc_limit_q, visc_limit_q, CS%diag) + if (CS%id_visc_limit_h_frac>0) call post_data(CS%id_visc_limit_h_frac, visc_limit_h_frac, CS%diag) + if (CS%id_visc_limit_q_frac>0) call post_data(CS%id_visc_limit_q_frac, visc_limit_q_frac, CS%diag) + if (CS%id_visc_limit_h_flag>0) call post_data(CS%id_visc_limit_h_flag, visc_limit_h_flag, CS%diag) + if (CS%id_visc_limit_q_flag>0) call post_data(CS%id_visc_limit_q_flag, visc_limit_q_flag, CS%diag) + + if (CS%EY24_EBT_BS) then + if (CS%id_BS_coeff_h>0) call post_data(CS%id_BS_coeff_h, BS_coeff_h, CS%diag) + if (CS%id_BS_coeff_q>0) call post_data(CS%id_BS_coeff_q, BS_coeff_q, CS%diag) + endif if (CS%debug) then if (CS%Laplacian) then @@ -2013,6 +2218,16 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag) endif + if (CS%id_FrictWorkIntz_bh > 0) then !cyc + do j=js,je + do i=is,ie ; FrictWorkIntz_bh(i,j) = FrictWork_bh(i,j,1) ; enddo + do k=2,nz ; do i=is,ie + FrictWorkIntz_bh(i,j) = FrictWorkIntz_bh(i,j) + FrictWork_bh(i,j,k) + enddo ; enddo + enddo + call post_data(CS%id_FrictWorkIntz_bh, FrictWorkIntz_bh, CS%diag) + endif + if (present(ADp)) then ! Diagnostics of the fractional thicknesses times momentum budget terms ! 3D diagnostics of hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option. @@ -2200,8 +2415,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "If true, the Laplacian coefficient is locally limited "//& "to be stable with a better bounding than just BOUND_KH.", & default=CS%bound_Kh, do_not_log=.not.CS%Laplacian) + call get_param(param_file, mdl, "EY24_EBT_BS", CS%EY24_EBT_BS, & + "If true, use the the backscatter scheme (EBT mode with kill switch)"//& + "developed by Yankovsky et al. (2024). ", & + default=.false., do_not_log=.not.CS%Laplacian) if (.not.CS%Laplacian) CS%bound_Kh = .false. if (.not.CS%Laplacian) CS%better_bound_Kh = .false. + if (.not.(CS%Laplacian.and.use_MEKE)) CS%EY24_EBT_BS = .false. call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& "horizontal viscosity.", default=.false., & @@ -2352,6 +2572,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "viscosity bounds to the theoretical maximum for "//& "stability without considering other terms.", units="nondim", & default=0.8, do_not_log=.not.(CS%better_bound_Ah .or. CS%better_bound_Kh)) + call get_param(param_file, mdl, "KILL_SWITCH_COEF", CS%KS_coef, & + "A nondimensional coefficient on the biharmonic viscosity that "// & + "sets the kill switch for backscatter. Default is 1.0.", units="nondim", & + default=1.0, do_not_log=.not.(CS%EY24_EBT_BS)) call get_param(param_file, mdl, "NOSLIP", CS%no_slip, & "If true, no slip boundary conditions are used; otherwise "//& "free slip boundary conditions are assumed. The "//& @@ -2417,6 +2641,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) fail_if_missing=.true.) Idt = 1.0 / dt endif + call get_param(param_file, mdl, "KILL_SWITCH_TIMESCALE", CS%KS_timescale, & + "A timescale for computing the CFL limit for viscosity "// & + "that determines when backscatter is shut off. Default is DT.", & + default= dt , units="s", scale=US%s_to_T, do_not_log=.not.(CS%EY24_EBT_BS)) + if (CS%no_slip .and. CS%biharmonic) & call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// & "at the same time in MOM.") @@ -2440,11 +2669,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0 ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0 ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0 - if (CS%bound_Kh .or. CS%better_bound_Kh) then + if (CS%bound_Kh .or. CS%better_bound_Kh .or. CS%EY24_EBT_BS) then ALLOC_(CS%Kh_Max_xx(Isd:Ied,Jsd:Jed)) ; CS%Kh_Max_xx(:,:) = 0.0 ALLOC_(CS%Kh_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_Max_xy(:,:) = 0.0 endif - if (CS%Smagorinsky_Kh) then + if (CS%Smagorinsky_Kh .or. CS%EY24_EBT_BS) then ALLOC_(CS%Laplac2_const_xx(isd:ied,jsd:jed)) ; CS%Laplac2_const_xx(:,:) = 0.0 ALLOC_(CS%Laplac2_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac2_const_xy(:,:) = 0.0 endif @@ -2502,6 +2731,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%Ah_Max_xx(isd:ied,jsd:jed)) ; CS%Ah_Max_xx(:,:) = 0.0 ALLOC_(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy(:,:) = 0.0 endif + if (CS%EY24_EBT_BS) then + ALLOC_(CS%Ah_Max_xx_KS(isd:ied,jsd:jed)) ; CS%Ah_Max_xx_KS(:,:) = 0.0 + ALLOC_(CS%Ah_Max_xy_KS(IsdB:IedB,JsdB:JedB)) ; CS%Ah_Max_xy_KS(:,:) = 0.0 + endif if (CS%Smagorinsky_Ah) then ALLOC_(CS%Biharm_const_xx(isd:ied,jsd:jed)) ; CS%Biharm_const_xx(:,:) = 0.0 ALLOC_(CS%Biharm_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Biharm_const_xy(:,:) = 0.0 @@ -2760,8 +2993,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) CS%DX_dyT(i,j)*(G%IdxCv(i,J)*v0v(i,J) + G%IdxCv(i,J-1)*v0v(i,J-1))) * & max(G%IdxCv(i,J)*G%IareaCv(i,J), G%IdxCv(i,J-1)*G%IareaCv(i,J-1)) ) ) CS%Ah_Max_xx(I,J) = 0.0 - if (denom > 0.0) & + if (denom > 0.0) then CS%Ah_Max_xx(I,J) = CS%bound_coef * 0.5 * Idt / denom + if (CS%EY24_EBT_BS) then + CS%Ah_Max_xx_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom) + endif + endif + enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq denom = max( & @@ -2774,8 +3012,13 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) CS%DY_dxBu(I,J)*(v0v(i+1,J)*G%IdyCv(i+1,J) + v0v(i,J)*G%IdyCv(i,J))) * & max(G%IdyCv(i,J)*G%IareaCv(i,J), G%IdyCv(i+1,J)*G%IareaCv(i+1,J)) ) ) CS%Ah_Max_xy(I,J) = 0.0 - if (denom > 0.0) & + if (denom > 0.0) then CS%Ah_Max_xy(I,J) = CS%bound_coef * 0.5 * Idt / denom + if (CS%EY24_EBT_BS) then + CS%Ah_Max_xy_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom) + endif + endif + enddo ; enddo if (CS%debug) then call hchksum(CS%Ah_Max_xx, "Ah_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T) @@ -2874,6 +3117,18 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T) CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, & 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim') + CS%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim') + CS%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, Time, & + 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim') + CS%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') + CS%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, Time, & + 'Value of the biharmonic viscosity limiter at h points', 'nondim') + CS%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, Time, & + 'Value of the biharmonic viscosity limiter at q points', 'nondim') if (CS%id_grid_Re_Ah > 0) & ! Compute the smallest biharmonic viscosity capable of modifying the @@ -2925,6 +3180,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', & 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) endif + + if (CS%EY24_EBT_BS) then + CS%id_BS_coeff_h = register_diag_field('ocean_model', 'BS_coeff_h', diag%axesTL, Time, & + 'Backscatter coefficient at h points', 'm2 s-1') + CS%id_BS_coeff_q = register_diag_field('ocean_model', 'BS_coeff_q', diag%axesBL, Time, & + 'Backscatter coefficient at q points', 'm2 s-1') + endif + CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& 'Integral work done by lateral friction terms. If GME is turned on, this '//& 'includes the GME contribution.', & @@ -2935,6 +3198,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) cmor_field_name='dispkexyfo', & cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',& cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction') + CS%id_FrictWork_bh = register_diag_field('ocean_model','FrictWork_bh',diag%axesTL,Time,& + 'Integral work done by the biharmonic lateral friction terms.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) !cyc + CS%id_FrictWorkIntz_bh = register_diag_field('ocean_model','FrictWorkIntz_bh',diag%axesT1,Time,& + 'Depth integrated work done by the biharmonic lateral friction', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) !cyc end subroutine hor_visc_init @@ -3165,6 +3434,9 @@ subroutine hor_visc_end(CS) if (CS%bound_Ah) then DEALLOC_(CS%Ah_Max_xx) ; DEALLOC_(CS%Ah_Max_xy) endif + if (CS%EY24_EBT_BS) then + DEALLOC_(CS%Ah_Max_xx_KS) ; DEALLOC_(CS%Ah_Max_xy_KS) + endif if (CS%Smagorinsky_Ah) then DEALLOC_(CS%Biharm_const_xx) ; DEALLOC_(CS%Biharm_const_xy) endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index eb1cf7b1cf..2e9620b752 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -113,6 +113,9 @@ module MOM_lateral_mixing_coeffs 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, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim] + real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0. + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & Laplac3_const_u !< Laplacian metric-dependent constants [L3 ~> m3] @@ -228,7 +231,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2] integer :: power_2 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz - integer :: i, j + integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -238,7 +241,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) if (CS%calculate_cg1) then if (.not. allocated(CS%cg1)) call MOM_error(FATAL, & "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.") - if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then + if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) then if (.not. allocated(CS%ebt_struct)) call MOM_error(FATAL, & "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.") if (CS%Resoln_use_ebt) then @@ -258,6 +261,11 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) call do_group_pass(CS%pass_cg1, G%Domain) endif + if (CS%BS_EBT_power>0.) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%BS_struct(i,j,k) = CS%ebt_struct(i,j,k)**CS%BS_EBT_power + enddo ; enddo ; enddo + endif ! Calculate and store the ratio between deformation radius and grid-spacing ! at h-points [nondim]. @@ -1249,6 +1257,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, uses the equivalent barotropic wave speed instead "//& "of first baroclinic wave for calculating the resolution fn.",& default=.false.) + call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, & + "Power to raise EBT vertical structure to when backscatter "// & + "has vertical structure.", units="nondim", default=0.0) call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, & "If true, uses the equivalent barotropic structure "//& "as the vertical structure of thickness diffusivity.",& @@ -1274,7 +1285,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) default=1.0e-17, units="s-1", scale=US%T_to_s) call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, & default=.false., do_not_log=.true.) - CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct + CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct & + .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0. CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE ! Indicate whether to calculate the Eady growth rate CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.) @@ -1298,7 +1310,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.") endif - if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) then + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & + .or. CS%BS_EBT_power>0.) then in_use = .true. call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, & "The depth below which N2 is monotonized to avoid stratification "//& @@ -1307,6 +1320,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) units="m", default=-1.0, scale=GV%m_to_H) allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif + allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0) + CS%BS_struct(:,:,:) = 1.0 if (CS%use_stored_slopes) then if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then @@ -1655,8 +1670,9 @@ end subroutine VarMix_init subroutine VarMix_end(CS) type(VarMix_CS), intent(inout) :: CS - if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct) & + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.) & deallocate(CS%ebt_struct) + if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) if (CS%use_stored_slopes) then deallocate(CS%slope_x) From 90e9ccb180d89243508d8750914e8091efa9c62d Mon Sep 17 00:00:00 2001 From: Elizabeth Yankovsky Date: Fri, 23 Aug 2024 14:07:23 -0400 Subject: [PATCH 2/5] Adding suggestions --- src/parameterizations/lateral/MOM_MEKE.F90 | 22 ++-- .../lateral/MOM_MEKE_types.F90 | 2 +- .../lateral/MOM_hor_visc.F90 | 117 +++++++++--------- 3 files changed, 71 insertions(+), 70 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index c2ba54c659..34334b155a 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -236,7 +236,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). - real :: sfac ! A factor needed to compute damping due to Strang splitting [nondim]c + real :: sfac ! A factor needed to compute damping due to Strang splitting [nondim] + real :: Isfac ! Inverse of sfac [nondim] logical :: use_drag_rate ! Flag to indicate drag_rate is finite integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features @@ -412,7 +413,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (allocated(MEKE%mom_src)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - !src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) src(i,j) = src(i,j) - CS%MEKE_FrCoeff*I_mass(i,j)*MEKE%mom_src(i,j) & - (CS%MEKE_bhFrCoeff-CS%MEKE_FrCoeff)*I_mass(i,j)*MEKE%mom_src_bh(i,j) src_mom_lp(i,j) = - CS%MEKE_FrCoeff*I_mass(i,j)*(MEKE%mom_src(i,j)-MEKE%mom_src_bh(i,j)) @@ -483,12 +483,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (MEKE%MEKE(i,j) < 0.) ldamping = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative + Isfac = 1. / (1. + sdt_damp * ldamping) MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) ldamping_Strang1(i,j) = ldamping - src_GM(i,j) = src_GM(i,j) / (1.0 + sdt_damp*ldamping) - src_mom_lp(i,j) = src_mom_lp(i,j) / (1.0 + sdt_damp*ldamping) - src_mom_bh(i,j) = src_mom_bh(i,j) / (1.0 + sdt_damp*ldamping) + src_GM(i,j) = src_GM(i,j) * Isfac + src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac + src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac sfac = ( 1.0 + sdt_damp*ldamping ) src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) enddo ; enddo @@ -661,13 +662,14 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (MEKE%MEKE(i,j) < 0.) ldamping = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative + Isfac = 1. / (1. + sdt_damp*ldamping) MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) - src_GM(i,j) = src_GM(i,j) / (1.0 + sdt_damp*ldamping) - src_mom_lp(i,j) = src_mom_lp(i,j) / (1.0 + sdt_damp*ldamping) - src_mom_bh(i,j) = src_mom_bh(i,j) / (1.0 + sdt_damp*ldamping) - src_adv(i,j) = src_adv(i,j) / (1.0 + sdt_damp*ldamping) - src_mom_K4(i,j) = src_mom_K4(i,j) / (1.0 + sdt_damp*ldamping) + src_GM(i,j) = src_GM(i,j) * Isfac + src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac + src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac + src_adv(i,j) = src_adv(i,j) * Isfac + src_mom_K4(i,j) = src_mom_K4(i,j) * Isfac sfac = ( 1.0 + sdt_damp*ldamping_Strang1(i,j) ) * ( 1.0 + sdt_damp*ldamping ) src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) enddo ; enddo diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index 890ee58824..a95578848d 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -12,7 +12,7 @@ module MOM_MEKE_types real, allocatable :: mom_src(:,:) !< MEKE source from lateral friction in the !! momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: mom_src_bh(:,:) !< MEKE source from the biharmonic part of the lateral friction in the - !! momentum equations [R Z L2 T-3 ~> W m-2]. !cyc + !! momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: GME_snk(:,:) !< MEKE sink from GME backscatter in the momentum equations [R Z L2 T-3 ~> W m-2]. real, allocatable :: Kh(:,:) !< The MEKE-derived lateral mixing coefficient [L2 T-1 ~> m2 s-1]. real, allocatable :: Kh_diff(:,:) !< Uses the non-MEKE-derived thickness diffusion coefficient to diffuse diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 0be1674b99..d5eaaea468 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -64,8 +64,8 @@ module MOM_hor_visc !! for stability without considering other terms [nondim]. !! The default is 0.8. real :: KS_coef !< A nondimensional coefficient on the biharmonic viscosity that sets the - !! kill switch for backscatter. Default is 1.0. - real :: KS_timescale !< A timescale for computing CFL limit for turning off backscatter (~DT) + !! kill switch for backscatter. Default is 1.0 [nondim]. + real :: KS_timescale !< A timescale for computing CFL limit for turning off backscatter [T ~> s]. logical :: backscatter_underbound !< If true, the bounds on the biharmonic viscosity are allowed !! to increase where the Laplacian viscosity is negative (due to !! backscatter parameterizations) beyond the largest timestep-dependent @@ -237,7 +237,7 @@ module MOM_hor_visc integer :: id_vort_xy_q = -1, id_div_xx_h = -1 integer :: id_sh_xy_q = -1, id_sh_xx_h = -1 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1 - integer :: id_FrictWork_bh = -1, id_FrictWorkIntz_bh = -1 !cyc + integer :: id_FrictWork_bh = -1, id_FrictWorkIntz_bh = -1 integer :: id_FrictWork_GME = -1 integer :: id_normstress = -1, id_shearstress = -1 integer :: id_visc_limit_h = -1, id_visc_limit_q = -1 @@ -325,7 +325,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [L2 T-2 ~> m2 s-2] bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] - FrictWorkIntz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R L2 T-3 ~> W m-2] !cyc + FrictWorkIntz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R L2 T-3 ~> W m-2] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1] @@ -371,9 +371,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1] sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1] GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1] - visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter - visc_limit_q_flag, & ! determines whether backscatter is shut off - visc_limit_q_frac, & ! determines how close backscatter is to shutting off + visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter [nondim] + visc_limit_q_flag, & ! determines whether backscatter is shut off [nondim] + visc_limit_q_frac, & ! determines how close backscatter is to shutting off [nondim] BS_coeff_q, & ! A diagnostic array of the backscatter coefficient [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) :: & @@ -387,7 +387,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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_bh, & ! work done by the biharmonic MKE dissipation mechanisms [R L2 T-3 ~> W m-2] !cyc + FrictWork_bh, & ! work done by the biharmonic 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] sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] @@ -397,9 +397,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] GME_coeff_h, & ! GME coefficient at h-points [L2 T-1 ~> m2 s-1] - visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter - visc_limit_h_flag, & ! determines whether backscatter is shut off - visc_limit_h_frac ! determines how close backscatter is to shutting off + visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter [nondim] + visc_limit_h_flag, & ! determines whether backscatter is shut off [nondim] + visc_limit_h_frac ! determines how close backscatter is to shutting off [nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & @@ -1414,7 +1414,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if ((CS%id_grid_Re_Ah>0)) then + if (CS%id_grid_Re_Ah > 0) then do j=js,je ; do i=is,ie KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) grid_Ah = max(Ah(i,j), CS%min_grid_Ah) @@ -1437,24 +1437,24 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx ! Backscatter using MEKE - if (CS%EY24_EBT_BS) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (visc_limit_h_flag(i,j,k) > 0) then - Kh_BS(i,j) = 0. - else - Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) - endif - enddo ; enddo + if (CS%EY24_EBT_BS) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_h_flag(i,j,k) > 0) then + Kh_BS(i,j) = 0. + else + Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) + endif + enddo ; enddo + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx_BS(i,j) = -Kh_BS(i,j) * sh_xx(i,j) + enddo ; enddo + + if (CS%id_BS_coeff_h>0) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - str_xx_BS(i,j) = -Kh_BS(i,j) * sh_xx(i,j) + BS_coeff_h(i,j,k) = Kh_BS(i,j) enddo ; enddo - - if (CS%id_BS_coeff_h>0) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - BS_coeff_h(i,j,k) = Kh_BS(i,j) - enddo ; enddo - endif + endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = str_xx(i,j) + str_xx_BS(i,j) @@ -1782,29 +1782,28 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif ! Get Ah at q points and biharmonic part of str_xy - ! Backscatter using MEKE - if (CS%EY24_EBT_BS) then - !if (grid_Re_Ah(i,j,k) < 1.) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (visc_limit_q_flag(I,J,k) > 0) then - Kh_BS(I,J) = 0. - else - Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & - (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & - ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & - (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) - endif - enddo ; enddo + ! Backscatter using MEKE + if (CS%EY24_EBT_BS) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (visc_limit_q_flag(I,J,k) > 0) then + Kh_BS(I,J) = 0. + else + Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) + endif + enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - str_xy_BS(I,J) = -Kh_BS(I,J) * (sh_xy(I,J)) - enddo ; enddo + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy_BS(I,J) = -Kh_BS(I,J) * (sh_xy(I,J)) + enddo ; enddo - if (CS%id_BS_coeff_q>0) then - do J=js-1,Jeq ; do I=is-1,Ieq - BS_coeff_q(I,J,k) = Kh_BS(I,J) - enddo ; enddo - endif + if (CS%id_BS_coeff_q>0) then + do J=js-1,Jeq ; do I=is-1,Ieq + BS_coeff_q(I,J,k) = Kh_BS(I,J) + enddo ; enddo + endif do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = str_xy(I,J) + str_xy_BS(I,J) @@ -1915,8 +1914,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) ! This is the old formulation that includes energy diffusion if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork(i,j,k) = 0 - FrictWork_bh(i,j,k) = 0 + FrictWork(i,j,k) = 0. + FrictWork_bh(i,j,k) = 0. else FrictWork(i,j,k) = GV%H_to_RZ * ( & (str_xx(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & @@ -1934,7 +1933,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, (u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & +(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) ) ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion !cyc + ! This is the old formulation that includes energy diffusion FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & @@ -2079,7 +2078,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (k==1) then do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = 0. - MEKE%mom_src_bh(i,j) = 0. !cyc + MEKE%mom_src_bh(i,j) = 0. enddo ; enddo if (allocated(MEKE%GME_snk)) then do j=js,je ; do i=is,ie @@ -2127,7 +2126,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, + (str_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) - MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + GV%H_to_RZ * ( & !cyc + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + GV%H_to_RZ * ( & ((bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & -(bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + 0.25*(((bhstr_xy(I,J)-RoScl*bhstr_xy(I,J)) * & @@ -2143,11 +2142,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) enddo ; enddo - else !cyc + else do j=js,je ; do i=is,ie MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k) - MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) !cyc + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k) enddo ; enddo endif ! MEKE%backscatter_Ro_c @@ -2167,7 +2166,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag) if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag) if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag) - if (CS%id_FrictWork_bh>0) call post_data(CS%id_FrictWork_bh, FrictWork_bh, CS%diag) !cyc + if (CS%id_FrictWork_bh>0) call post_data(CS%id_FrictWork_bh, FrictWork_bh, CS%diag) if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag) if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag) if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag) @@ -2218,7 +2217,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag) endif - if (CS%id_FrictWorkIntz_bh > 0) then !cyc + if (CS%id_FrictWorkIntz_bh > 0) then do j=js,je do i=is,ie ; FrictWorkIntz_bh(i,j) = FrictWork_bh(i,j,1) ; enddo do k=2,nz ; do i=is,ie @@ -3200,10 +3199,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction') CS%id_FrictWork_bh = register_diag_field('ocean_model','FrictWork_bh',diag%axesTL,Time,& 'Integral work done by the biharmonic lateral friction terms.', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) !cyc + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) CS%id_FrictWorkIntz_bh = register_diag_field('ocean_model','FrictWorkIntz_bh',diag%axesT1,Time,& 'Depth integrated work done by the biharmonic lateral friction', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) !cyc + 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2) end subroutine hor_visc_init From c2e69518e5e578a81df263efc6ad4a05c2eeaa07 Mon Sep 17 00:00:00 2001 From: Wenda Zhang Date: Mon, 26 Aug 2024 18:36:08 -0400 Subject: [PATCH 3/5] Modifications to MOM_hor_visc: - Separate FrictWork and FrictWork_bh loops - Simplified the computation for MEKE%mom_src and MEKE%mom_src_bh when MEKE%backscatter_Ro_c /= 0. (cherry picked from commit 8ffc6a8b9b23de54d782831e419eae3872c11ac6) --- .../lateral/MOM_hor_visc.F90 | 90 ++++++++----------- 1 file changed, 39 insertions(+), 51 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index d282ec5116..0b1996d2d1 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1915,8 +1915,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) ! This is the old formulation that includes energy diffusion if (visc_limit_h_flag(i,j,k) > 0) then - FrictWork(i,j,k) = 0. - FrictWork_bh(i,j,k) = 0. + FrictWork(i,j,k) = 0 else FrictWork(i,j,k) = GV%H_to_RZ * ( & (str_xx(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & @@ -1933,29 +1932,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, +str_xy(I,J-1)*( & (u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & +(v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1) )) ) ) - ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) - ! This is the old formulation that includes energy diffusion - FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & - (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - + 0.25*((bhstr_xy(I,J) * & - ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & - + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & - + bhstr_xy(I-1,J-1) * & - ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & - + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & - + (bhstr_xy(I-1,J) * & - ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & - + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & - + bhstr_xy(I,J-1) * & - ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & - + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) endif enddo ; enddo else ; do j=js,je ; do i=is,ie if (visc_limit_h_flag(i,j,k) > 0) then FrictWork(i,j,k) = 0 - FrictWork_bh(i,j,k) = 0 else FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((str_xx(i,j)*CS%dy2h(i,j) * ( & @@ -1985,6 +1966,40 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) & - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) ) + endif + enddo ; enddo ; endif + endif + + if (CS%id_FrictWork_bh>0 .or. CS%id_FrictWorkIntz_bh > 0 .or. allocated(MEKE%mom_src_bh)) then + if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie + ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork_bh(i,j,k) = 0 + else + ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) + ! This is the old formulation that includes energy diffusion !cyc + FrictWork_bh(i,j,k) = GV%H_to_RZ * ( & + (bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & + - bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & + + 0.25*((bhstr_xy(I,J) * & + ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & + + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)) & + + bhstr_xy(I-1,J-1) * & + ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & + + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & + + (bhstr_xy(I-1,J) * & + ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & + + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & + + bhstr_xy(I,J-1) * & + ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & + + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + endif + enddo ; enddo + else ; do j=js,je ; do i=is,ie + if (visc_limit_h_flag(i,j,k) > 0) then + FrictWork_bh(i,j,k) = 0 + else ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v) FrictWork_bh(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( & ((bhstr_xx(i,j)*CS%dy2h(i,j) * ( & @@ -2019,7 +2034,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, - if (CS%use_GME) then if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v) @@ -2112,36 +2126,10 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + GV%H_to_RZ * ( & - ((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - -(str_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - + 0.25*(((str_xy(I,J)-RoScl*bhstr_xy(I,J)) * & - ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & - + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) & - + (str_xy(I-1,J-1)-RoScl*bhstr_xy(I-1,J-1)) * & - ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & - + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & - + ((str_xy(I-1,J)-RoScl*bhstr_xy(I-1,J)) * & - ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & - + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & - + (str_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & - ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & - + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) - MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + GV%H_to_RZ * ( & - ((bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & - -(bhstr_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & - + 0.25*(((bhstr_xy(I,J)-RoScl*bhstr_xy(I,J)) * & - ((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J) & - + (v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J) ) & - + (bhstr_xy(I-1,J-1)-RoScl*bhstr_xy(I-1,J-1)) * & - ((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1) & - + (v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)) ) & - + ((bhstr_xy(I-1,J)-RoScl*bhstr_xy(I-1,J)) * & - ((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J) & - + (v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)) & - + (bhstr_xy(I,J-1)-RoScl*bhstr_xy(I,J-1)) * & - ((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1) & - + (v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)) ) ) ) + MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + (FrictWork(i,j,k) - RoScl*FrictWork_bh(i,j,k)) + MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + & + (FrictWork_bh(i,j,k) - RoScl*FrictWork_bh(i,j,k)) + enddo ; enddo else From 18ef0412f394ca1b319d52c53a095c7292d79d8c Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 6 Sep 2024 17:46:45 -0400 Subject: [PATCH 4/5] MEKE: No-division implementation of src_btm_drag This patch modifies the calculation of `src_btm_drag` to avoid additional division operations. We first note the following variable renames: * `Isfac` is now `damping`, since it describes the net damping or reduction of any damped fields. `sfac` has been removed since it no longer appears in the expressions. * `ldamping` is now `damp_rate`. This is primarily to avoid confusion with `damping`, but also to more accurately describe its role. * `ldamping_Strang1` is replaced to `damp_rate_s1`. It is used for a similar purpose (cacheing of values from the first stage) but now holds a slightly different value. The following modifications were made. 1. The conditional split of `sdt` into `sdt_damping` substeps is quantified by a new term, `damp_step` which is equal to 0.5 or 1. `sdt_damp` is computed as `sdt * damp_step`. The presence of `sdt_damp / sdt` in our expressions is replaced with `damp_step`, which avoids an extra division. 2. The first expression (using new notation) was originally sfac = 1 + sdt_damp * damp_rate D = M * (1 - sfac) / (sdt * sfac) where `D` is `src_btm_drag` and `M` is `MEKE_current(i,j)` This has been transformed to D = - M * (sdt_damp * damp_rate) / (sdt * (1 + sdt_damp * damp_rate)) = - M * (sdt_damp / sft) * damp_rate * (1 / (1 + sdt_damp * damp_rate)) = - M * damp_step * damp_rate * damping This new expression for `D` no longer requires a division. 3. In the second stage expression for `src_btm_drag`, we again use `damp_rate` to replace `ldamping`. We also use `damp_rate1` to denote the original value of `ldaming_Strang1`. sfac = (1 + sdt_damp * damp_rate1) * (1 + sdt_damp * damp_rate) D = M * (1 - sfac) / (sdt * sfac) Using `damping1` to denote the damping from the first stage, `D` can be transformed as follows. D = -M * (sdt_damp * damp_rate + sdt_damp * damp_rate1 + sdt_damp**2 * damp_rate * damp_rate1) / (sdt * (1 + sdt_damp * damp_rate) * (1 + sdt_damp * damp_rate1)) = -M * (sdt_damp / sdt) * (1 / (1 + sdt_damp * damp_rate)) * (damp_rate + damp_rate1 + sdt_damp * damp_rate * damp_rate1) / (1 + sdt_damp * damp_rate1) = -M * damp_step * damping * (damp_rate * (1 + sdt_damp * damp_rate1) + damp_rate1) / (1 + sdt_damp * damp_rate1) = -M * damp_step * damping * (damp_rate + damp_rate1 * damping1) We now store `damp_rate1 * damping1` in `damp_rate_s1(:,:)` rather than just `damping1`, and can reduce this expression to D = -M * damp_step * damping * (damp_rate + damp_rate_s1(i,j)) which requires no division. It also requires no additional read or writes. Assuming there are no mistakes here, this should only cause negligible bitwise differences in the results. Also, comments have been added which explain that we cannot modify the existing expressions for `MEKE%MEKE`, since they would alter the bitwise results of existing runs. --- src/parameterizations/lateral/MOM_MEKE.F90 | 75 ++++++++++++++-------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 34334b155a..b6b0cf00be 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -201,7 +201,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3). src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3). - ldamping_Strang1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1]. + damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1]. MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1] drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1]. @@ -233,11 +233,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: cdrag2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6] real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1] real :: mass_neglect ! A negligible mass [R Z ~> kg m-2]. - real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). - real :: sfac ! A factor needed to compute damping due to Strang splitting [nondim] - real :: Isfac ! Inverse of sfac [nondim] + real :: damp_step ! Size of damping timestep relative to sdt [nondim] + real :: damp_rate ! The MEKE damping rate [T-1 ~> s-1]. + real :: damping ! The net damping of a field after sdt_damp [nondim] logical :: use_drag_rate ! Flag to indicate drag_rate is finite integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features @@ -287,7 +287,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! With a depth-dependent (and possibly strong) damping, it seems ! advisable to use Strang splitting between the damping and diffusion. - sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt + damp_step = 1. + if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5 + sdt_damp = sdt * damp_step ! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg] if (CS%MEKE_advection_factor>0.) then @@ -479,19 +481,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! First stage of Strang splitting !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. + damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative - Isfac = 1. / (1. + sdt_damp * ldamping) - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) - ldamping_Strang1(i,j) = ldamping - src_GM(i,j) = src_GM(i,j) * Isfac - src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac - src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac - sfac = ( 1.0 + sdt_damp*ldamping ) - src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) + + damping = 1. / (1. + sdt_damp * damp_rate) + + ! NOTE: MEKE%MEKE should use `damping` but we must preserve the existing + ! expression for bit reproducibility + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate) + MEKE_decay(i,j) = damp_rate * G%mask2dT(i,j) + + src_GM(i,j) = src_GM(i,j) * damping + src_mom_lp(i,j) = src_mom_lp(i,j) * damping + src_mom_bh(i,j) = src_mom_bh(i,j) * damping + + src_btm_drag(i,j) = - MEKE_current(i,j) * ( & + damp_step * (damp_rate * damping) & + ) + + ! Store the effective damping rate if sdt is split + if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) & + damp_rate_s1(i,j) = damp_rate * damping enddo ; enddo if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then @@ -658,20 +670,27 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. + damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. ! notice that the above line ensures a damping only if MEKE is positive, ! while leaving MEKE unchanged if it is negative - Isfac = 1. / (1. + sdt_damp*ldamping) - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) - src_GM(i,j) = src_GM(i,j) * Isfac - src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac - src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac - src_adv(i,j) = src_adv(i,j) * Isfac - src_mom_K4(i,j) = src_mom_K4(i,j) * Isfac - sfac = ( 1.0 + sdt_damp*ldamping_Strang1(i,j) ) * ( 1.0 + sdt_damp*ldamping ) - src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) ) + + damping = 1. / (1. + sdt_damp * damp_rate) + + ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve + ! the existing expression for bit reproducibility. + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate) + MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j) + + src_GM(i,j) = src_GM(i,j) * damping + src_mom_lp(i,j) = src_mom_lp(i,j) * damping + src_mom_bh(i,j) = src_mom_bh(i,j) * damping + src_adv(i,j) = src_adv(i,j) * damping + src_mom_K4(i,j) = src_mom_K4(i,j) * damping + + src_btm_drag(i,j) = -MEKE_current(i,j) * ( & + damp_step * damping * (damp_rate + damp_rate_s1(i,j)) & + ) enddo ; enddo endif endif ! MEKE_KH>=0 From 174e8a3811a85f5b28913cd5c46509de0d74ace1 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sat, 7 Sep 2024 21:21:13 -0400 Subject: [PATCH 5/5] MEKE: Remove redundant if-block in Strang split In the second stage of Strang splitting, there is a check for `sdt > sdt_damping`, following a check for nonzero `KH` or `K4`. However, `sdt_damp` can only be less than `sdt` when these are nonzero. Otherwise, `sdt_damping` equals `sdt`. Best as I can tell, this check is unnecessary and potentially problematic for optimization. This patch removes the redundant check. --- src/parameterizations/lateral/MOM_MEKE.F90 | 58 +++++++++++----------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index b6b0cf00be..1614acc075 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -659,40 +659,38 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h ! Second stage of Strang splitting if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.0) then - if (sdt>sdt_damp) then - ! Recalculate the drag rate, since MEKE has changed. - if (use_drag_rate) then - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & - cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) - enddo ; enddo - endif + ! Recalculate the drag rate, since MEKE has changed. + if (use_drag_rate) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. - ! notice that the above line ensures a damping only if MEKE is positive, - ! while leaving MEKE unchanged if it is negative - - damping = 1. / (1. + sdt_damp * damp_rate) - - ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve - ! the existing expression for bit reproducibility. - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate) - MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j) - - src_GM(i,j) = src_GM(i,j) * damping - src_mom_lp(i,j) = src_mom_lp(i,j) * damping - src_mom_bh(i,j) = src_mom_bh(i,j) * damping - src_adv(i,j) = src_adv(i,j) * damping - src_mom_K4(i,j) = src_mom_K4(i,j) * damping - - src_btm_drag(i,j) = -MEKE_current(i,j) * ( & - damp_step * damping * (damp_rate + damp_rate_s1(i,j)) & - ) + drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & + cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo endif + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) damp_rate = 0. + ! notice that the above line ensures a damping only if MEKE is positive, + ! while leaving MEKE unchanged if it is negative + + damping = 1. / (1. + sdt_damp * damp_rate) + + ! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve + ! the existing expression for bit reproducibility. + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate) + MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j) + + src_GM(i,j) = src_GM(i,j) * damping + src_mom_lp(i,j) = src_mom_lp(i,j) * damping + src_mom_bh(i,j) = src_mom_bh(i,j) * damping + src_adv(i,j) = src_adv(i,j) * damping + src_mom_K4(i,j) = src_mom_K4(i,j) * damping + + src_btm_drag(i,j) = -MEKE_current(i,j) * ( & + damp_step * damping * (damp_rate + damp_rate_s1(i,j)) & + ) + enddo ; enddo endif ! MEKE_KH>=0 if (CS%debug) then