diff --git a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 index f41b0f7c0f..ae7f9744a9 100644 --- a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 +++ b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 @@ -16,6 +16,7 @@ module MOM_Zanna_Bolton use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE use MOM_ANN, only : ANN_init, ANN_apply, ANN_end, ANN_CS use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_lateral_mixing_coeffs, only : VarMix_CS implicit none ; private @@ -58,13 +59,19 @@ module MOM_Zanna_Bolton !! 1 - memory for eddy forcing !! 2 - memory for large scale flow !! 3 - memory for both eddy forcing and large scale flow + integer :: ZB_memory_dynamical_scale !< Select which dynamical length scale to infer memory time + !! scale from: + !! 0 - no memory + !! 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2) + !! 2 - grid scale + !! 3 - deformation radius real :: ZB_eddy_memory !< The eddy memory time scale [T ~> s] real :: ZB_eddy_memory_const !< The non-dimensional factor for the eddy memory time scale - !! that scales the deformation scale [nondim] + !! that scales the dynamical time scale [nondim] !! Typical range: 1-10 real :: ZB_large_scale_memory !< The memory time scale for the velocity gradients [T ~> s] real :: ZB_large_scale_memory_const !< The non-dimensional factor for the memory time scale - !! for the velocity gradients that scales the deformation scale [nondim] + !! for the velocity gradients that scales the dynamical time scale [nondim] !! Typical range: 1-10 real :: ZB_max_memory !< The maximum eddy memory time scale; default: 1 year [T ~> s] real :: DT !< The baroclinic model time step [T ~> s] @@ -76,8 +83,7 @@ module MOM_Zanna_Bolton !! points including metric terms [T-1 ~> s-1] vort_xy, & !< Vertical vorticity (dv/dx - du/dy) in q (CORNER) !! points including metric terms [T-1 ~> s-1] - hq !< Thickness in CORNER points [H ~> m or kg m-2] - + hq !< Thickness in CORNER points [H ~> m or kg m-2] real, dimension(:,:,:), allocatable :: & Txx, & !< Subgrid stress xx component in h [L2 T-2 ~> m2 s-2] Tyy, & !< Subgrid stress yy component in h [L2 T-2 ~> m2 s-2] @@ -95,7 +101,8 @@ module MOM_Zanna_Bolton real, dimension(:,:), allocatable :: & kappa_h, & !< Scaling coefficient in h points [L2 ~> m2] - kappa_q !< Scaling coefficient in q points [L2 ~> m2] + kappa_q, & !< Scaling coefficient in q points [L2 ~> m2] + Rd_dx_h !< The deformation radius compared with the grid spacing [nondim]. real, allocatable :: & ICoriolis_h(:,:), & !< Inverse Coriolis parameter at h points [T ~> s] @@ -166,6 +173,7 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. logical, intent(out) :: use_ZB2020 !< If true, turns on ZB scheme. + real :: subroundoff_Cor !> A negligible parameter which avoids division by zero !! but small compared to Coriolis parameter [T-1 ~> s-1] @@ -252,6 +260,13 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) "\t 2 - memory for large scale flow\n" //& "\t 3 - memory for both eddy forcing and large scale flow", default=0) + call get_param(param_file, mdl, "ZB_MEMORY_DYNAMICAL_SCALE", CS%ZB_memory_dynamical_scale, & + "Select which dynamical length scale to infer memory time scale from:\n" //& + "\t 0 - no memory\n" //& + "\t 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2)\n" //& + "\t 2 - grid scale\n" //& + "\t 3 - deformation radius", default=0) + call get_param(param_file, mdl, "ZB_EDDY_MEMORY", CS%ZB_eddy_memory, & "The eddy memory time scale. "//& "A value of zero means no memory.", units="s", default=0.0) @@ -304,6 +319,13 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) "MOM_Zanna_Bolton: only one of ZB_EDDY_MEMORY_CONST or ZB_EDDY_MEMORY can be larger than 0") if (CS%ZB_large_scale_memory_const > 0 .AND. CS%ZB_large_scale_memory > 0) call MOM_error(FATAL, & "MOM_Zanna_Bolton: only one of ZB_LARGE_SCALE_MEMORY_CONST or ZB_LARGE_SCALE_MEMORY can be larger than 0") + + if (CS%ZB_eddy_memory_const > 0 .OR. CS%ZB_large_scale_memory_const > 0) then + if (CS%ZB_memory_dynamical_scale == 0) then + call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_MEMORY_DYNAMICAL_SCALE must be greater than 0 if ZB_EDDY_MEMORY_CONST > 0 or ZB_LARGE_SCALE_MEMORY_CONST > 0") + endif + endif ! Register fields for output from this module. CS%diag => diag @@ -408,6 +430,9 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) endif endif + if (CS%ZB_memory_dynamical_scale == 3) then + allocate(CS%Rd_dx_h(SZI_(G),SZJ_(G))); CS%Rd_dx_h(:,:) = 0. + endif ! Precomputing the scaling coefficient @@ -579,7 +604,7 @@ end subroutine ZB_copy_gradient_and_thickness !> Baroclinic Zanna-Bolton-2020 parameterization, see !! eq. 6 in https://laurezanna.github.io/files/Zanna-Bolton-2020.pdf subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & - dx2h, dy2h, dx2q, dy2q) + dx2h, dy2h, dx2q, dy2q, VarMix) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. @@ -605,6 +630,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & real, dimension(SZIB_(G),SZJB_(G)), & intent(in) :: dx2q, & !< dx^2 at q points [L2 ~> m2] dy2q !< dy^2 at q points [L2 ~> m2] + type(VarMix_CS), target, intent(in) :: VarMix !< Variable mixing coefficients real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & ZB2020u !< Zonal acceleration due to convergence of @@ -631,7 +657,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & ! Introduce memory into velocity gradients if specified if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then - if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.false.) + if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy=.false.) endif ! Compute the stress tensor given the @@ -646,7 +672,7 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & call filter_stress(G, GV, CS) if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then - if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.true.) + if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy=.true.) endif ! Compute acceleration from a given stress tensor @@ -698,7 +724,7 @@ end subroutine Zanna_Bolton_2020 -subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) +subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, VarMix, eddy) type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -708,15 +734,17 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + type(VarMix_CS), target, intent(in) :: VarMix !< Variable mixing coefficients logical, intent(in) :: eddy !< If true, step forward memory on eddy forcing; ! if false, step forward memory on velocity gradients ! Local variables real :: sh_xy_h !< Shearing strain interpolated to h point [T-1 ~> s-1] + real :: vort_xy_h !< Relative vorticity interpolated to h point [T-1 ~> s-1] real :: sh_xx_q !< Horizontal tension interpolated to q point [T-1 ~> s-1] real :: D !< deformation rate [T-1 ~> s-1] - real :: tau !< eddy memory time scale tau [T ~> s] + real :: tau, tau0 !< eddy memory time scale tau [T ~> s] integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq integer :: i, j, k character(len=160) :: mesg ! The text of an error message @@ -759,19 +787,44 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) call advection(u, v, G, GV, CS, q=CS%vort_xy_with_memory) endif endif + + if (CS%ZB_memory_dynamical_scale == 3) then + if (allocated(VarMix%Rd_dx_h) .AND. allocated(CS%Rd_dx_h)) then + do j=js-1,je+1 ; do i=is-1,ie+1 + CS%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j) + enddo ; enddo + endif + endif do k=1,nz - ! relaxation: this is the Eulerian part of the averaging do j=js-1,je+1 ; do i=is-1,ie+1 - ! compute eddy memory time scale at cell center - sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) & - + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) ) - D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h) + ! compute dynamical time scale tau0 + if (CS%ZB_memory_dynamical_scale > 0) then + ! compute eddy memory time scale at cell center + sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) & + + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) ) + vort_xy_h = 0.25 * ( (CS%vort_xy(I-1,J-1,k) + CS%vort_xy(I,J,k)) & + + (CS%vort_xy(I-1,J,k) + CS%vort_xy(I,J-1,k)) ) + if (CS%ZB_memory_dynamical_scale == 1) then + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h) + tau0 = 1 / (D + 1.0/CS%ZB_max_memory) + elseif (CS%ZB_memory_dynamical_scale == 2) then + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h) + tau0 = 1 / (D + 1.0/CS%ZB_max_memory) + elseif (CS%ZB_memory_dynamical_scale == 3) then + + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h) + tau0 = CS%Rd_dx_h(i,j) / (D + CS%Rd_dx_h(i,j)/CS%ZB_max_memory) + endif + + endif + + ! relaxation: this is the Eulerian part of the averaging if (eddy) then if (CS%ZB_eddy_memory_const > 0.0) then - tau = CS%ZB_eddy_memory_const / (D + 1.0/CS%ZB_max_memory) + tau = CS%ZB_eddy_memory_const * tau0 else tau = CS%ZB_eddy_memory endif @@ -788,7 +841,7 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) else if (CS%ZB_large_scale_memory_const > 0.0) then - tau = CS%ZB_large_scale_memory_const / (D + 1.0/CS%ZB_max_memory) + tau = CS%ZB_large_scale_memory_const * tau0 else tau = CS%ZB_large_scale_memory endif diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index d8b229041e..3c913467c4 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1692,7 +1692,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, if (CS%use_ZB2020) then call Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS%ZB2020, & - CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q) + CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q, VarMix) endif end subroutine horizontal_viscosity diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 1bf416b00a..97b393baa9 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1109,6 +1109,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The ! default value is roughly (pi / (the age of the universe)). logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use + integer :: ZB_memory_dynamical_scale integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use ! for remapping. Values below 20190101 recover the remapping @@ -1180,6 +1181,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, the viscosity contribution from MEKE is scaled by "//& "the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere. if (.not.use_MEKE) Resoln_scaled_MEKE_visc = .false. + call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, & "If true, uses the equivalent barotropic wave speed instead "//& "of first baroclinic wave for calculating the resolution fn.",& @@ -1211,6 +1213,13 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE + call get_param(param_file, mdl, "ZB_MEMORY_DYNAMICAL_SCALE", ZB_memory_dynamical_scale, & + "Select which dynamical length scale to infer memory time scale from:\n" //& + "\t 0 - no memory\n" //& + "\t 1 - deformation scale (old), which probably differs from option 1 by factor of sqrt(2)\n" //& + "\t 2 - grid scale\n" //& + "\t 3 - deformation radius", default=0) + CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (ZB_memory_dynamical_scale == 3) ! 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.) call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, &