From a9dd5f2a3df1bae665caf1d70f66725955bb8ca0 Mon Sep 17 00:00:00 2001 From: Wenda Zhang Date: Fri, 26 Jul 2024 23:59:35 -0400 Subject: [PATCH 1/3] Added SQG vertical structure in Varmix to provide vertical profile for diffusivities - added function calc_sqg_struct in MOM_lateral_mixing_coeffs to compute sqg_struct - added sqg_expo to set the exponent of sqg_struct - to use sqg_struct for the backscatter, set BS_use_sqg=true, sqg_expo>0., and BS_EBT_power=0. - if SQG_USE_MEKE=True, use the eddy length scale from MEKE to compute sqg_struct - added eddy length scale Le in MEKE if SQG_USE_MEKE=True - added MEKE%Le into restart file if SQG_USE_MEKE=True - added MEKE in Varmix - registered N2_u and N2_v diagnostics when SQG_EXPO>0 (cherry picked from commit 6d3df0541c33d6f6d1f9fcb695f1a1eb961ec1b3) --- src/core/MOM.F90 | 6 +- src/parameterizations/lateral/MOM_MEKE.F90 | 34 ++++- .../lateral/MOM_MEKE_types.F90 | 1 + .../lateral/MOM_lateral_mixing_coeffs.F90 | 123 +++++++++++++++++- 4 files changed, 153 insertions(+), 11 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e6cd3cf747..7a6d32f1aa 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -750,7 +750,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%VarMix%use_variable_mixing) then call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag) - call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix) + call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) call calc_depth_function(G, CS%VarMix) call disable_averaging(CS%diag) endif @@ -1898,7 +1898,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif @@ -1925,7 +1925,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 9ebe8ae734..4ae2cd4089 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -124,6 +124,7 @@ module MOM_MEKE logical :: debug !< If true, write out checksums of data for debugging integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default), !! read in from a file, or inferred via a neural network + logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure. 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 @@ -400,6 +401,13 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h call hchksum(LmixScale, 'MEKE LmixScale', G%HI, unscale=US%L_to_m) endif + if (allocated(MEKE%Le)) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%Le(i,j) = LmixScale(i,j) + enddo ; enddo + endif + ! Aggregate sources of MEKE (background, frictional and GM) !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -757,7 +765,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif - if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) then + if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & + .or. allocated(MEKE%Le)) then call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_Kh, G%Domain) call cpu_clock_end(CS%id_clock_pass) @@ -1425,6 +1434,9 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME "computing beta in the expression of Rhines scale. Use 1 if full "//& "topographic beta effect is considered; use 0 if it's completely ignored.", & units="nondim", default=0.0) + call get_param(param_file, mdl, "SQG_USE_MEKE", CS%sqg_use_MEKE, & + "If true, the eddy scale of MEKE is used for the SQG vertical structure ",& + default=.false.) ! Nonlocal module parameters call get_param(param_file, mdl, "CDRAG", cdrag, & @@ -1531,6 +1543,7 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME CS%id_clock_pass = cpu_clock_id('(Ocean continuity halo updates)', grain=CLOCK_ROUTINE) + ! Detect whether this instance of MEKE_init() is at the beginning of a run ! or after a restart. If at the beginning, we will initialize MEKE to a local ! equilibrium. @@ -1538,6 +1551,12 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (coldStart) CS%initialize = .false. if (CS%initialize) call MOM_error(WARNING, & "MEKE_init: Initializing MEKE with a local equilibrium balance.") + if (.not.query_initialized(MEKE%Le, "MEKE_Le", restart_CS) .and. allocated(MEKE%Le)) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + MEKE%Le(i,j) = sqrt(G%areaT(i,j)) + enddo ; enddo + endif ! Set up group passes. In the case of a restart, these fields need a halo update now. if (allocated(MEKE%MEKE)) then @@ -1548,8 +1567,10 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (allocated(MEKE%Kh)) call create_group_pass(CS%pass_Kh, MEKE%Kh, G%Domain) if (allocated(MEKE%Ku)) call create_group_pass(CS%pass_Kh, MEKE%Ku, G%Domain) if (allocated(MEKE%Au)) call create_group_pass(CS%pass_Kh, MEKE%Au, G%Domain) + if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain) - if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au)) & + if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & + .or. allocated(MEKE%Le)) & call do_group_pass(CS%pass_Kh, G%Domain) end function MEKE_init @@ -1839,6 +1860,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim] logical :: Use_KH_in_MEKE logical :: useMEKE + logical :: sqg_use_MEKE integer :: isd, ied, jsd, jed ! Determine whether this module will be used @@ -1853,6 +1875,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au) Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE) + sqg_use_MEKE = .false. ; call read_param(param_file,"SQG_USE_MEKE", sqg_use_MEKE) if (.not. useMEKE) return @@ -1884,6 +1907,12 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", & units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T) endif + if (sqg_use_MEKE) then + allocate(MEKE%Le(isd:ied,jsd:jed), source=0.0) + call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, & + longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", & + units="m", conversion=US%L_to_m) + endif if (Use_Kh_in_MEKE) then allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0) call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, & @@ -1918,6 +1947,7 @@ subroutine MEKE_end(MEKE) 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) + if (allocated(MEKE%Le)) deallocate(MEKE%Le) end subroutine MEKE_end !> \namespace mom_meke diff --git a/src/parameterizations/lateral/MOM_MEKE_types.F90 b/src/parameterizations/lateral/MOM_MEKE_types.F90 index a95578848d..38bdd72452 100644 --- a/src/parameterizations/lateral/MOM_MEKE_types.F90 +++ b/src/parameterizations/lateral/MOM_MEKE_types.F90 @@ -24,6 +24,7 @@ module MOM_MEKE_types !! backscatter from unresolved eddies (see Jansen and Held, 2014). real, allocatable :: Au(:,:) !< The MEKE-derived lateral biharmonic viscosity !! coefficient [L4 T-1 ~> m4 s-1]. + real, allocatable :: Le(:,:) !< Eddy length scale [L m] ! Parameters real :: KhTh_fac = 1.0 !< Multiplier to map Kh(MEKE) to KhTh [nondim] diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e3979fe35a..de72cb68d0 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -18,6 +18,8 @@ module MOM_lateral_mixing_coeffs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init use MOM_open_boundary, only : ocean_OBC_type +use MOM_MEKE_types, only : MEKE_type + implicit none ; private @@ -112,9 +114,12 @@ 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 :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim] + real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure 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 :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0 + logical :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & @@ -163,6 +168,7 @@ module MOM_lateral_mixing_coeffs integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1 integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1 integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1 + integer :: id_sqg_struct=-1, id_BS_struct=-1 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. !>@} @@ -173,7 +179,7 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function -public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function +public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function, calc_sqg_struct contains @@ -214,13 +220,14 @@ subroutine calc_depth_function(G, CS) end subroutine calc_depth_function !> Calculates and stores the non-dimensional resolution functions -subroutine calc_resoln_function(h, tv, G, GV, US, CS) +subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure + type(MEKE_type), intent(in) :: MEKE !< MEKE struct ! Local variables ! Depending on the power-function being used, dimensional rescaling may be limited, so some @@ -230,6 +237,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) real :: cg1_v ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1]. real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2] integer :: power_2 + real :: dt !< Time increment [T ~> s] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -261,10 +269,19 @@ 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%sqg_expo>0.0) then + call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) + call pass_var(CS%sqg_struct, 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 + elseif (CS%BS_use_sqg) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%BS_struct(i,j,k) = CS%sqg_struct(i,j,k) + enddo ; enddo ; enddo endif ! Calculate and store the ratio between deformation radius and grid-spacing @@ -460,6 +477,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) if (query_averaging_enabled(CS%diag)) then if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag) + if (CS%id_BS_struct > 0) call post_data(CS%id_BS_struct, CS%BS_struct, CS%diag) endif if (CS%debug) then @@ -470,6 +488,77 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) end subroutine calc_resoln_function +!> Calculates and stores functions of SQG mode +subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv ! s] + type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct + type(MEKE_type), intent(in) :: MEKE !< MEKE struct + ! type(ocean_OBC_type), pointer :: OBC !< Open +! boundaries control structure. + + ! Local variables + real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & + e ! The interface heights relative to mean sea level [Z ~> m]. + real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] + real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] + real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G), SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1] + real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2] + real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m] +! real, dimension(SZK_(GV)) :: zs ! Stretched vertical coordinate [m] + real, dimension(SZI_(G), SZJ_(G)) :: Le ! Eddy length scale [m] + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& + "Module must be initialized before it is used.") + + call find_eta(h, tv, G, GV, US, e, halo_size=2) + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & + CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, & + dzSxN=dzSxN, dzSyN=dzSyN, halo=1) + + do j=js,je ; do i=is,ie + CS%sqg_struct(i,j,1) = 1.0 + enddo ; enddo + if (allocated(MEKE%Le)) then + do j=js,je ; do i=is,ie + Le(i,j) = MEKE%Le(i,j) + f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8) + enddo ; enddo + else + do j=js,je ; do i=is,ie + Le(i,j) = sqrt(G%areaT(i,j)) + f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8) + enddo ; enddo + endif + do k=2,nz ; do j=js,je ; do i=is,ie + N2 = max(0.25*(N2_u(I-1,j,k) + N2_u(I,j,k) + N2_v(i,J-1,k) + N2_v(i,J,k)),0.0) + dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k))*N2**0.5/f(i,j) +! dzs = -N2**0.5/f(i,j)*dzc + CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*dzc/Le(i,j)) + enddo ; enddo ; enddo + + + if (query_averaging_enabled(CS%diag)) then + if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag) + if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag) + if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag) + endif + +end subroutine calc_sqg_struct + !> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al. !! style scaling of diffusivity subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) @@ -1260,6 +1349,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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, "BS_USE_SQG", CS%BS_use_sqg, & + "If true, the SQG vertical structure is used for backscatter "//& + "on the condition that BS_EBT_power=0", & + default=.false.) + if (CS%BS_EBT_power>0.) CS%BS_use_sqg = .false. + call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, & + "Nondimensional exponent coeffecient of the SQG mode "// & + "that is used for the vertical struture of diffusivities.", units="nondim", default=0.0) + if (CS%sqg_expo==0.) CS%BS_use_sqg = .false. 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.",& @@ -1320,6 +1418,10 @@ 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 + if (CS%sqg_expo>0.0) then + allocate(CS%sqg_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 @@ -1335,7 +1437,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) endif endif - if (CS%use_stored_slopes) then + if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then ! CS%calculate_Eady_growth_rate=.true. in_use = .true. allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0) @@ -1418,7 +1520,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 'm2', conversion=US%L_to_m**2) endif - if (CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) then + if (CS%sqg_expo>0.0) then + CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & + 'Vertical structure of SQG mode', 'nondim') + endif + + CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, & + 'Vertical structure of backscatter', 'nondim') + + if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) .or. (CS%sqg_expo>0.0)) then CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', & 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) @@ -1672,9 +1782,10 @@ subroutine VarMix_end(CS) 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 (CS%sqg_expo>0.0) deallocate(CS%sqg_struct) if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) - if (CS%use_stored_slopes) then + if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then deallocate(CS%slope_x) deallocate(CS%slope_y) endif From 83092217654d455e5c864788b9e4d18c527fc29c Mon Sep 17 00:00:00 2001 From: Wenda Zhang Date: Thu, 3 Oct 2024 17:15:44 -0400 Subject: [PATCH 2/3] Compute vertical structures for khth, khtr, backscatter, and kdgl90 all in VarMix - Vertical structures including khth_struct, khtr_struct, BS_struct, and kdgl90_struct are now computed in VarMix - Each diffusivity/viscosity have two vertical structure options, equivalent barotropic (EBT) and surface quasigeostrophic (SQG) mode structures - KHTH_USE_EBT_STRUCT, KHTR_USE_EBT_STRUCT, KDGL90_USE_EBT_STRUCT and BS_EBT_POWER parameters, which already existed, still control whether to use the EBT structure for khth, khtr, kdgl90, and backscatter, respectively - Added KHTH_USE_SQG_STRUCT, KHTR_USE_SQG_STRUCT, KDGL90_USE_SQG_STRUCT and BS_USE_SQG parameters to control whether to use the SQG structure for khth, khtr, kdgl90, and backscatter, respectively - If neither EBT nor SQG is called, no vertical structure will be used for that diffusivity/viscosity - An error will be called if both EBT and SQG structures are called for the same diffusivity/viscosity - Added dt as an input of calc_resoln_function. dt is needed for calc_sqg_struct called in calc_resoln_function --- src/core/MOM.F90 | 6 +- src/parameterizations/lateral/MOM_MEKE.F90 | 4 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 154 +++++++++++++++--- .../lateral/MOM_thickness_diffuse.F90 | 14 +- .../vertical/MOM_vert_friction.F90 | 16 +- src/tracer/MOM_neutral_diffusion.F90 | 32 ++-- src/tracer/MOM_tracer_hor_diff.F90 | 35 ++-- 7 files changed, 190 insertions(+), 71 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7a6d32f1aa..e05d2e1619 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -750,7 +750,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%VarMix%use_variable_mixing) then call enable_averages(cycle_time, Time_start + real_to_time(US%T_to_s*cycle_time), CS%diag) - call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) + call calc_resoln_function(h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt) call calc_depth_function(G, CS%VarMix) call disable_averaging(CS%diag) endif @@ -1898,7 +1898,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif @@ -1925,7 +1925,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS if (.not. skip_diffusion) then if (CS%VarMix%use_variable_mixing) then call pass_var(CS%h, G%Domain) - call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE) + call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix, CS%MEKE, dt_offline) call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 4ae2cd4089..fd5cfbeee3 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -124,7 +124,7 @@ module MOM_MEKE logical :: debug !< If true, write out checksums of data for debugging integer :: eke_src !< Enum specifying whether EKE is stepped forward prognostically (default), !! read in from a file, or inferred via a neural network - logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure. + logical :: sqg_use_MEKE !< If True, use MEKE%Le for the SQG vertical structure. 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 @@ -1912,7 +1912,7 @@ subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) call register_restart_field(MEKE%Le, "MEKE_Le", .false., restart_CS, & longname="Eddy length scale from Mesoscale Eddy Kinetic Energy", & units="m", conversion=US%L_to_m) - endif + endif if (Use_Kh_in_MEKE) then allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0) call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, & diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index de72cb68d0..7d1ccd060e 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -52,6 +52,14 @@ module MOM_lateral_mixing_coeffs !! as the vertical structure of thickness diffusivity. logical :: kdgl90_use_ebt_struct !< If true, uses the equivalent barotropic structure !! as the vertical structure of diffusivity in the GL90 scheme. + logical :: kdgl90_use_sqg_struct !< If true, uses the surface quasigeostrophic structure + !! as the vertical structure of diffusivity in the GL90 scheme. + logical :: khth_use_sqg_struct !< If true, uses the surface quasigeostrophic structure + !! as the vertical structure of thickness diffusivity. + logical :: khtr_use_ebt_struct !< If true, uses the equivalent barotropic structure + !! as the vertical structure of tracer diffusivity. + logical :: khtr_use_sqg_struct !< If true, uses the surface quasigeostrophic structure + !! as the vertical structure of tracer diffusivity. logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first !! baroclinic wave speed and populate CS%cg1. !! This parameter is set depending on other parameters. @@ -117,6 +125,9 @@ module MOM_lateral_mixing_coeffs real, allocatable :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim] real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure to scale diffusivities with [nondim] real, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim] + real, allocatable :: khth_struct(:,:,:) !< Vertical structure function used in thickness diffusivity [nondim] + real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim] + real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [nondim] real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0. real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0 logical :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure. @@ -168,7 +179,8 @@ module MOM_lateral_mixing_coeffs integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1 integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1 integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1 - integer :: id_sqg_struct=-1, id_BS_struct=-1 + integer :: id_sqg_struct=-1, id_BS_struct=-1, id_khth_struct=-1, id_khtr_struct=-1 + integer :: id_kdgl90_struct=-1 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. !>@} @@ -220,7 +232,7 @@ subroutine calc_depth_function(G, CS) end subroutine calc_depth_function !> Calculates and stores the non-dimensional resolution functions -subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) +subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] @@ -228,6 +240,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure type(MEKE_type), intent(in) :: MEKE !< MEKE struct + real, intent(in) :: dt !< Time increment [T ~> s] ! Local variables ! Depending on the power-function being used, dimensional rescaling may be limited, so some @@ -237,7 +250,6 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) real :: cg1_v ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1]. real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2] integer :: power_2 - real :: dt !< Time increment [T ~> s] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -249,7 +261,8 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) 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 .or. CS%BS_EBT_power>0.) then + if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & + .or. CS%khtr_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 @@ -269,12 +282,17 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) call do_group_pass(CS%pass_cg1, G%Domain) endif - if (CS%sqg_expo>0.0) then + if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & + .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) call pass_var(CS%sqg_struct, G%Domain) endif - if (CS%BS_EBT_power>0.) then + if (CS%BS_EBT_power>0. .and. CS%BS_use_sqg) then + call MOM_error(FATAL, & + "calc_resoln_function: BS_EBT_POWER>0. & + and BS_USE_SQG=True cannot be set together") + elseif (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 @@ -284,6 +302,48 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) enddo ; enddo ; enddo endif + if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT & + and KHTH_USE_SQG_STRUCT can be true") + elseif (CS%khth_use_ebt_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%khth_struct(i,j,k) = CS%ebt_struct(i,j,k) + enddo ; enddo ; enddo + elseif (CS%khth_use_sqg_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%khth_struct(i,j,k) = CS%sqg_struct(i,j,k) + enddo ; enddo ; enddo + endif + + if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT & + and KHTR_USE_SQG_STRUCT can be true") + elseif (CS%khtr_use_ebt_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%khtr_struct(i,j,k) = CS%ebt_struct(i,j,k) + enddo ; enddo ; enddo + elseif (CS%khtr_use_sqg_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%khtr_struct(i,j,k) = CS%sqg_struct(i,j,k) + enddo ; enddo ; enddo + endif + + if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT & + and KD_GL90_USE_SQG_STRUCT can be true") + elseif (CS%kdgl90_use_ebt_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%kdgl90_struct(i,j,k) = CS%ebt_struct(i,j,k) + enddo ; enddo ; enddo + elseif (CS%kdgl90_use_sqg_struct) then + do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + CS%kdgl90_struct(i,j,k) = CS%sqg_struct(i,j,k) + enddo ; enddo ; enddo + endif + ! Calculate and store the ratio between deformation radius and grid-spacing ! at h-points [nondim]. if (CS%calculate_rd_dx) then @@ -478,6 +538,9 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE) if (query_averaging_enabled(CS%diag)) then if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag) if (CS%id_BS_struct > 0) call post_data(CS%id_BS_struct, CS%BS_struct, CS%diag) + if (CS%id_khth_struct > 0) call post_data(CS%id_khth_struct, CS%khth_struct, CS%diag) + if (CS%id_khtr_struct > 0) call post_data(CS%id_khtr_struct, CS%khtr_struct, CS%diag) + if (CS%id_kdgl90_struct > 0) call post_data(CS%id_kdgl90_struct, CS%kdgl90_struct, CS%diag) endif if (CS%debug) then @@ -534,18 +597,24 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) do j=js,je ; do i=is,ie Le(i,j) = MEKE%Le(i,j) f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & - G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8) + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8/US%T_to_s) enddo ; enddo else do j=js,je ; do i=is,ie Le(i,j) = sqrt(G%areaT(i,j)) f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & - G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8) + G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8*US%T_to_s) enddo ; enddo endif + if (CS%debug) then + call hchksum(Le, 'SQG length scale', G%HI, unscale=US%L_to_m) + call hchksum(f, 'Coriolis at h point', G%HI, unscale=US%s_to_T) + call uvchksum( 'MEKE LmixScale', dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.) + endif do k=2,nz ; do j=js,je ; do i=is,ie N2 = max(0.25*(N2_u(I-1,j,k) + N2_u(I,j,k) + N2_v(i,J-1,k) + N2_v(i,J,k)),0.0) - dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k))*N2**0.5/f(i,j) + dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k)) * & + N2**0.5/f(i,j)*US%Z_to_L ! dzs = -N2**0.5/f(i,j)*dzc CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*dzc/Le(i,j)) enddo ; enddo ; enddo @@ -1353,19 +1422,33 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) "If true, the SQG vertical structure is used for backscatter "//& "on the condition that BS_EBT_power=0", & default=.false.) - if (CS%BS_EBT_power>0.) CS%BS_use_sqg = .false. call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, & "Nondimensional exponent coeffecient of the SQG mode "// & - "that is used for the vertical struture of diffusivities.", units="nondim", default=0.0) - if (CS%sqg_expo==0.) CS%BS_use_sqg = .false. + "that is used for the vertical struture of diffusivities.", units="nondim", default=1.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.",& default=.false.) + call get_param(param_file, mdl, "KHTH_USE_SQG_STRUCT", CS%khth_use_sqg_struct, & + "If true, uses the surface quasigeostrophic structure "//& + "as the vertical structure of thickness diffusivity.",& + default=.false.) + call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%khtr_use_ebt_struct, & + "If true, uses the equivalent barotropic structure "//& + "as the vertical structure of tracer diffusivity.",& + default=.false.) + call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", CS%khtr_use_sqg_struct, & + "If true, uses the surface quasigeostrophic structure "//& + "as the vertical structure of tracer diffusivity.",& + default=.false.) call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", CS%kdgl90_use_ebt_struct, & "If true, uses the equivalent barotropic structure "//& "as the vertical structure of diffusivity in the GL90 scheme.",& default=.false.) + call get_param(param_file, mdl, "KD_GL90_USE_SQG_STRUCT", CS%kdgl90_use_sqg_struct, & + "If true, uses the equivalent barotropic structure "//& + "as the vertical structure of diffusivity in the GL90 scheme.",& + default=.false.) call get_param(param_file, mdl, "KHTH_SLOPE_CFF", KhTh_Slope_Cff, & "The nondimensional coefficient in the Visbeck formula "//& "for the interface depth diffusivity", units="nondim", default=0.0) @@ -1409,7 +1492,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) endif if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & - .or. CS%BS_EBT_power>0.) then + .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) 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 "//& @@ -1418,13 +1501,23 @@ 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 - if (CS%sqg_expo>0.0) then - allocate(CS%sqg_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%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then + allocate(CS%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0) + endif + + if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then + allocate(CS%khtr_struct(isd:ied, jsd:jed, gv%ke), source=0.0) + endif + + if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then + allocate(CS%kdgl90_struct(isd:ied, jsd:jed, gv%ke), source=0.0) + endif + if (CS%use_stored_slopes) then if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", CS%Visbeck_S_max, & @@ -1520,15 +1613,29 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 'm2', conversion=US%L_to_m**2) endif - if (CS%sqg_expo>0.0) then - CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & + CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & 'Vertical structure of SQG mode', 'nondim') + if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & + .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then + allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, & 'Vertical structure of backscatter', 'nondim') + if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then + CS%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, Time, & + 'Vertical structure of thickness diffusivity', 'nondim') + endif + if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then + CS%id_khtr_struct = register_diag_field('ocean_model', 'khtr_struct', diag%axesTl, Time, & + 'Vertical structure of tracer diffusivity', 'nondim') + endif + if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then + CS%id_kdgl90_struct = register_diag_field('ocean_model', 'kdgl90_struct', diag%axesTl, Time, & + 'Vertical structure of GL90 diffusivity', 'nondim') + endif - if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) .or. (CS%sqg_expo>0.0)) then + if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) ) then CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', & 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) @@ -1780,10 +1887,13 @@ 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 .or. CS%BS_EBT_power>0.) & - deallocate(CS%ebt_struct) - if (CS%sqg_expo>0.0) deallocate(CS%sqg_struct) + if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct & + .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct) + if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct) if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) + if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) deallocate(CS%khth_struct) + if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) deallocate(CS%khtr_struct) + if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) deallocate(CS%kdgl90_struct) if (CS%use_stored_slopes .or. CS%sqg_expo>0.0) then deallocate(CS%slope_x) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index b2f022ad18..daea75e73d 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -179,7 +179,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp ! to layer centers [L2 T-1 ~> m2 s-1] real :: KH_v_lay(SZI_(G),SZJ_(G)) ! Diagnostic of isopycnal height diffusivities at v-points averaged ! to layer centers [L2 T-1 ~> m2 s-1] - logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck + logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_vert_struct, use_Visbeck logical :: use_QG_Leith integer :: i, j, k, is, ie, js, je, nz @@ -198,7 +198,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp endif use_VarMix = .false. ; Resoln_scaled = .false. ; use_stored_slopes = .false. - khth_use_ebt_struct = .false. ; use_Visbeck = .false. ; use_QG_Leith = .false. + khth_use_vert_struct = .false. ; use_Visbeck = .false. ; use_QG_Leith = .false. Depth_scaled = .false. if (VarMix%use_variable_mixing) then @@ -206,7 +206,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp Resoln_scaled = VarMix%Resoln_scaled_KhTh Depth_scaled = VarMix%Depth_scaled_KhTh use_stored_slopes = VarMix%use_stored_slopes - khth_use_ebt_struct = VarMix%khth_use_ebt_struct + khth_use_vert_struct = allocated(VarMix%khth_struct) use_Visbeck = VarMix%use_Visbeck use_QG_Leith = VarMix%use_QG_Leith_GM if (allocated(VarMix%cg1)) cg1 => VarMix%cg1 @@ -298,10 +298,10 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp KH_u(I,j,1) = min(KH_u_CFL(I,j), Khth_loc_u(I,j)) enddo ; enddo - if (khth_use_ebt_struct) then + if (khth_use_vert_struct) then !$OMP do do K=2,nz+1 ; do j=js,je ; do I=is-1,ie - KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i+1,j,k-1) ) enddo ; enddo ; enddo else !$OMP do @@ -394,10 +394,10 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo ; enddo endif - if (khth_use_ebt_struct) then + if (khth_use_vert_struct) then !$OMP do do K=2,nz+1 ; do J=js-1,je ; do i=is,ie - KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i,j+1,k-1) ) enddo ; enddo ; enddo else !$OMP do diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index e241d191b2..f119ccb040 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -616,7 +616,7 @@ subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, do_i, z_i, j, G, GV, CS, Va !! otherwise they are v-points. ! local variables - logical :: kdgl90_use_ebt_struct + logical :: kdgl90_use_vert_struct ! use vertical structure for GL90 coefficient integer :: i, k, is, ie, nz, Isq, Ieq real :: f2 !< Squared Coriolis parameter at a velocity grid point [T-2 ~> s-2]. real :: h_neglect ! A vertical distance that is so small it is usually lost in roundoff error @@ -629,9 +629,9 @@ subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, do_i, z_i, j, G, GV, CS, Va nz = GV%ke h_neglect = GV%dZ_subroundoff - kdgl90_use_ebt_struct = .false. + kdgl90_use_vert_struct = .false. if (VarMix%use_variable_mixing) then - kdgl90_use_ebt_struct = VarMix%kdgl90_use_ebt_struct + kdgl90_use_vert_struct = allocated(VarMix%kdgl90_struct) endif if (work_on_u) then @@ -647,8 +647,9 @@ subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, do_i, z_i, j, G, GV, CS, Va else a_cpl_gl90(I,K) = f2 * CS%kappa_gl90 / GV%g_prime(K) endif - if (kdgl90_use_ebt_struct) then - a_cpl_gl90(I,K) = a_cpl_gl90(I,K) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + if (kdgl90_use_vert_struct) then + a_cpl_gl90(I,K) = a_cpl_gl90(I,K) * 0.5 * & + ( VarMix%kdgl90_struct(i,j,k-1) + VarMix%kdgl90_struct(i+1,j,k-1) ) endif endif ! botfn determines when a point is within the influence of the GL90 bottom boundary layer, @@ -671,8 +672,9 @@ subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, do_i, z_i, j, G, GV, CS, Va else a_cpl_gl90(i,K) = f2 * CS%kappa_gl90 / GV%g_prime(K) endif - if (kdgl90_use_ebt_struct) then - a_cpl_gl90(i,K) = a_cpl_gl90(i,K) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + if (kdgl90_use_vert_struct) then + a_cpl_gl90(i,K) = a_cpl_gl90(i,K) * 0.5 * & + ( VarMix%kdgl90_struct(i,j,k-1) + VarMix%kdgl90_struct(i,j+1,k-1) ) endif endif ! botfn determines when a point is within the influence of the GL90 bottom boundary layer, diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 051ac446db..1aaf7409d2 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -57,8 +57,8 @@ module MOM_neutral_diffusion logical :: tapering = .false. !< If true, neutral diffusion linearly decays towards zero within a !! transition zone defined using boundary layer depths. Only available when !! interior_only=true. - logical :: KhTh_use_ebt_struct !< If true, uses the equivalent barotropic structure - !! as the vertical structure of tracer diffusivity. + logical :: KhTh_use_vert_struct !< If true, uses vertical structure + !! for tracer diffusivity. logical :: use_unmasked_transport_bug !< If true, use an older form for the accumulation of !! neutral-diffusion transports that were unmasked, as used prior to Jan 2018. real, allocatable, dimension(:,:) :: hbl !< Boundary layer depth [H ~> m or kg m-2] @@ -67,7 +67,7 @@ module MOM_neutral_diffusion !! at cell interfaces [nondim] real, allocatable, dimension(:) :: coeff_r !< Non-dimensional coefficient in the right column, !! at cell interfaces [nondim] - ! Array used when KhTh_use_ebt_struct is true + ! Array used when KhTh_use_vert_struct is true real, allocatable, dimension(:,:,:) :: Coef_h !< Coef_x and Coef_y averaged at t-points [L2 ~> m2] ! Positions of neutral surfaces in both the u, v directions real, allocatable, dimension(:,:,:) :: uPoL !< Non-dimensional position with left layer uKoL-1, u-point [nondim] @@ -153,6 +153,7 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, logical :: boundary_extrap ! Indicate whether high-order boundary !! extrapolation should be used within boundary cells. logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm + logical :: KhTh_use_ebt_struct, KhTh_use_sqg_struct if (associated(CS)) then call MOM_error(FATAL, "neutral_diffusion_init called with associated control structure.") @@ -197,10 +198,14 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "a transition zone defined using boundary layer depths. "//& "Only applicable when NDIFF_INTERIOR_ONLY=True", default=.false.) endif - call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, & + call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", KhTh_use_ebt_struct, & "If true, uses the equivalent barotropic structure "//& "as the vertical structure of the tracer diffusivity.",& default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", KhTh_use_sqg_struct, & + "If true, uses the surface geostrophic structure "//& + "as the vertical structure of the tracer diffusivity.",& + default=.false.,do_not_log=.true.) call get_param(param_file, mdl, "NDIFF_USE_UNMASKED_TRANSPORT_BUG", CS%use_unmasked_transport_bug, & "If true, use an older form for the accumulation of neutral-diffusion "//& "transports that were unmasked, as used prior to Jan 2018. This is not "//& @@ -296,7 +301,8 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, endif endif - if (CS%KhTh_use_ebt_struct) & + CS%KhTh_use_vert_struct = KhTh_use_ebt_struct .or. KhTh_use_sqg_struct + if (CS%KhTh_use_vert_struct) & allocate(CS%Coef_h(G%isd:G%ied,G%jsd:G%jed,SZK_(GV)+1), source=0.) ! Store a rescaling factor for use in diagnostic messages. @@ -624,11 +630,11 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx ! Zonal flux of tracer in units that vary between a ! thickness times a concentration ([C H ~> degC m or degC kg m-2] for temperature) or a ! volume or mass times a concentration ([C H L2 ~> degC m3 or degC kg] for temperature), - ! depending on the setting of CS%KhTh_use_ebt_struct. + ! depending on the setting of CS%KhTh_use_vert_struct. real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx ! Meridional flux of tracer in units that vary between a ! thickness times a concentration ([C H ~> degC m or degC kg m-2] for temperature) or a ! volume or mass times a concentration ([C H L2 ~> degC m3 or degC kg] for temperature), - ! depending on the setting of CS%KhTh_use_ebt_struct. + ! depending on the setting of CS%KhTh_use_vert_struct. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tendency ! tendency array for diagnostics ! [H conc T-1 ~> m conc s-1 or kg m-2 conc s-1] ! For temperature these units are @@ -673,7 +679,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) endif endif - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then ! Compute Coef at h points CS%Coef_h(:,:,:) = 0. do j = G%jsc,G%jec ; do i = G%isc,G%iec @@ -706,7 +712,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) vFlx(:,:,:) = 0. ! x-flux - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then if (CS%tapering) then do j = G%jsc,G%jec ; do I = G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then @@ -770,7 +776,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) endif ! y-flux - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then if (CS%tapering) then do J = G%jsc-1,G%jec ; do i = G%isc,G%iec if (G%mask2dCv(i,J)>0.) then @@ -837,7 +843,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) ! Update the tracer concentration from divergence of neutral diffusive flux components, noting ! that uFlx and vFlx use an unexpected sign convention. - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then do j = G%jsc,G%jec ; do i = G%isc,G%iec if (G%mask2dT(i,j)>0.) then if (CS%ndiff_answer_date <= 20240330) then @@ -940,7 +946,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) ! Note sign corresponds to downgradient flux convention. if (tracer%id_dfx_2d > 0) then - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then do j = G%jsc,G%jec ; do I = G%isc-1,G%iec trans_x_2d(I,j) = 0. if (G%mask2dCu(I,j)>0.) then @@ -969,7 +975,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) ! Note sign corresponds to downgradient flux convention. if (tracer%id_dfy_2d > 0) then - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTh_use_vert_struct) then do J = G%jsc-1,G%jec ; do i = G%isc,G%iec trans_y_2d(i,J) = 0. if (G%mask2dCv(i,J)>0.) then diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 855c6bef30..3511b88f39 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -52,7 +52,7 @@ module MOM_tracer_hor_diff real :: max_diff_CFL !< If positive, locally limit the along-isopycnal !! tracer diffusivity to keep the diffusive CFL !! locally at or below this value [nondim]. - logical :: KhTh_use_ebt_struct !< If true, uses the equivalent barotropic structure + logical :: KhTr_use_vert_struct !< If true, uses the equivalent barotropic structure !! as the vertical structure of tracer diffusivity. logical :: Diffuse_ML_interior !< If true, diffuse along isopycnals between !! the mixed layer and the interior. @@ -218,6 +218,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ use_VarMix = VarMix%use_variable_mixing Resoln_scaled = VarMix%Resoln_scaled_KhTr use_Eady = CS%KhTr_Slope_Cff > 0. + CS%KhTr_use_vert_struct = allocated(VarMix%khtr_struct) endif call cpu_clock_begin(id_clock_pass) @@ -422,18 +423,18 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ enddo enddo enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_vert_struct) then do K=2,nz+1 do J=js-1,je do i=is,ie - Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) enddo enddo enddo do k=2,nz+1 do j=js,je do I=is-1,ie - Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) enddo enddo enddo @@ -478,18 +479,18 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ enddo enddo enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_vert_struct) then do K=2,nz+1 do J=js-1,je do i=is,ie - Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) enddo enddo enddo do k=2,nz+1 do j=js,je do I=is-1,ie - Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) enddo enddo enddo @@ -605,11 +606,11 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ do j=js,je ; do I=is-1,ie Kh_u(I,j,:) = G%mask2dCu(I,j)*Kh_u(I,j,1) enddo ; enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_vert_struct) then do K=2,nz+1 do j=js,je do I=is-1,ie - Kh_u(I,j,K) = Kh_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) ) + Kh_u(I,j,K) = Kh_u(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) ) enddo enddo enddo @@ -621,11 +622,11 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ do J=js-1,je ; do i=is,ie Kh_v(i,J,:) = G%mask2dCv(i,J)*Kh_v(i,J,1) enddo ; enddo - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_vert_struct) then do K=2,nz+1 do J=js-1,je do i=is,ie - Kh_v(i,J,K) = Kh_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) ) + Kh_v(i,J,K) = Kh_v(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) ) enddo enddo enddo @@ -647,9 +648,9 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_ (G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37) Kh_h(i,j,:) = normalize*G%mask2dT(i,j)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & (Kh_v(i,J-1,1)+Kh_v(i,J,1))) - if (CS%KhTh_use_ebt_struct) then + if (CS%KhTr_use_vert_struct) then do K=2,nz+1 - Kh_h(i,j,K) = normalize*G%mask2dT(i,j)*VarMix%ebt_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & + Kh_h(i,j,K) = normalize*G%mask2dT(i,j)*VarMix%khtr_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + & (Kh_v(i,J-1,1)+Kh_v(i,J,1))) enddo endif @@ -1630,10 +1631,10 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic call get_param(param_file, mdl, "KHTR", CS%KhTr, & "The background along-isopycnal tracer diffusivity.", & units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s) - call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, & - "If true, uses the equivalent barotropic structure "//& - "as the vertical structure of the tracer diffusivity.",& - default=.false.) +! call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, & +! "If true, uses the equivalent barotropic structure "//& +! "as the vertical structure of the tracer diffusivity.",& +! default=.false.) call get_param(param_file, mdl, "KHTR_SLOPE_CFF", CS%KhTr_Slope_Cff, & "The scaling coefficient for along-isopycnal tracer "//& "diffusivity using a shear-based (Visbeck-like) "//& From 4ed3dca5da817fc3603fc557f8d01592ed1b2134 Mon Sep 17 00:00:00 2001 From: Wenda Zhang Date: Fri, 11 Oct 2024 14:15:40 -0400 Subject: [PATCH 3/3] Bug fixes - Changed ()**0.5 to sqrt in sqg_struct, which solves the dimension error - Added if statement for using BS_struct in MOM_hor_visc - Added if statement for CS%sqg_struct<=0. - Changed the subroundoff of Coriolis parameter to 1e-40 in calc_sqg_struct - Changed parameter BS_USE_SQG to BS_USE_SQG_STRUCT --- src/parameterizations/lateral/MOM_MEKE.F90 | 4 +- .../lateral/MOM_hor_visc.F90 | 66 ++++++--- .../lateral/MOM_lateral_mixing_coeffs.F90 | 130 ++++++++++-------- 3 files changed, 120 insertions(+), 80 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index fd5cfbeee3..4bac09b7cc 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -766,7 +766,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & - .or. allocated(MEKE%Le)) then + .or. allocated(MEKE%Le)) then call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_Kh, G%Domain) call cpu_clock_end(CS%id_clock_pass) @@ -1570,7 +1570,7 @@ logical function MEKE_init(Time, G, GV, US, param_file, diag, dbcomms_CS, CS, ME if (allocated(MEKE%Le)) call create_group_pass(CS%pass_Kh, MEKE%Le, G%Domain) if (allocated(MEKE%Kh) .or. allocated(MEKE%Ku) .or. allocated(MEKE%Au) & - .or. allocated(MEKE%Le)) & + .or. allocated(MEKE%Le)) & call do_group_pass(CS%pass_Kh, G%Domain) end function MEKE_init diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 56a359857f..fa5bb40577 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -446,6 +446,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, logical :: use_MEKE_Ku logical :: use_MEKE_Au logical :: use_cont_huv + logical :: use_kh_struct integer :: is_vort, ie_vort, js_vort, je_vort ! Loop ranges for vorticity terms integer :: is_Kh, ie_Kh, js_Kh, je_Kh ! Loop ranges for thickness point viscosities integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -502,6 +503,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, if (CS%id_FrictWorkIntz > 0) find_FrictWork = .true. if (allocated(MEKE%mom_src)) find_FrictWork = .true. + use_kh_struct = allocated(VarMix%BS_struct) backscat_subround = 0.0 if (find_FrictWork .and. allocated(MEKE%mom_src) .and. (MEKE%backscatter_Ro_c > 0.0) .and. & (MEKE%backscatter_Ro_Pow /= 0.0)) & @@ -663,7 +665,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP is_vort, ie_vort, js_vort, je_vort, & !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, & - !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & + !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, use_kh_struct, & !$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, & @@ -1181,14 +1183,26 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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) * VarMix%BS_struct(i,j,k) - enddo ; enddo + if (use_kh_struct) then + 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) * 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) * VarMix%BS_struct(i,j,k) + enddo ; enddo + endif else - do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) - enddo ; enddo + 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) + 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) + enddo ; enddo + endif endif endif @@ -1443,7 +1457,11 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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) + if (use_kh_struct) then + Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k) + else + Kh_BS(i,j) = MEKE%Ku(i,j) + endif endif enddo ; enddo @@ -1618,10 +1636,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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)*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 + if (use_kh_struct) then + 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 + else + 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 + endif endif if (CS%anisotropic) & @@ -1789,10 +1814,17 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, 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))) ) + if (use_kh_struct) then + 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))) ) + else + Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + & + MEKE%Ku(i+1,j+1)) + & + (MEKE%Ku(i+1,j) + & + MEKE%Ku(i,j+1)) ) + endif endif enddo ; enddo diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 7d1ccd060e..cc8c0c48b8 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -129,8 +129,8 @@ module MOM_lateral_mixing_coeffs real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim] real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [nondim] real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0. - real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 0.0 - logical :: BS_use_sqg !< If true, use sqg_stuct for backscatter vertical structure. + real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 1.0 + logical :: BS_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & @@ -282,31 +282,23 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain) call do_group_pass(CS%pass_cg1, G%Domain) endif - if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & + if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) call pass_var(CS%sqg_struct, G%Domain) endif - if (CS%BS_EBT_power>0. .and. CS%BS_use_sqg) then - call MOM_error(FATAL, & - "calc_resoln_function: BS_EBT_POWER>0. & - and BS_USE_SQG=True cannot be set together") - elseif (CS%BS_EBT_power>0.) then + 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 - elseif (CS%BS_use_sqg) then + elseif (CS%BS_use_sqg_struct) then do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied CS%BS_struct(i,j,k) = CS%sqg_struct(i,j,k) enddo ; enddo ; enddo endif - if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then - call MOM_error(FATAL, & - "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT & - and KHTH_USE_SQG_STRUCT can be true") - elseif (CS%khth_use_ebt_struct) then + if (CS%khth_use_ebt_struct) then do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied CS%khth_struct(i,j,k) = CS%ebt_struct(i,j,k) enddo ; enddo ; enddo @@ -316,11 +308,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) enddo ; enddo ; enddo endif - if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then - call MOM_error(FATAL, & - "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT & - and KHTR_USE_SQG_STRUCT can be true") - elseif (CS%khtr_use_ebt_struct) then + if (CS%khtr_use_ebt_struct) then do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied CS%khtr_struct(i,j,k) = CS%ebt_struct(i,j,k) enddo ; enddo ; enddo @@ -330,11 +318,7 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, dt) enddo ; enddo ; enddo endif - if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then - call MOM_error(FATAL, & - "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT & - and KD_GL90_USE_SQG_STRUCT can be true") - elseif (CS%kdgl90_use_ebt_struct) then + if (CS%kdgl90_use_ebt_struct) then do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied CS%kdgl90_struct(i,j,k) = CS%ebt_struct(i,j,k) enddo ; enddo ; enddo @@ -561,8 +545,6 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) real, intent(in) :: dt !< Time increment [T ~> s] type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct type(MEKE_type), intent(in) :: MEKE !< MEKE struct - ! type(ocean_OBC_type), pointer :: OBC !< Open -! boundaries control structure. ! Local variables real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & @@ -574,13 +556,14 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] real, dimension(SZI_(G), SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1] - real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2] - real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m] -! real, dimension(SZK_(GV)) :: zs ! Stretched vertical coordinate [m] + real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2] + real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [m] + real :: f_subround ! The minimal resolved value of Coriolis parameter to prevent division by zero [T-1 ~> s-1] real, dimension(SZI_(G), SZJ_(G)) :: Le ! Eddy length scale [m] integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + f_subround = 1.0e-40 * US%s_to_T if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& "Module must be initialized before it is used.") @@ -590,34 +573,32 @@ subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE) CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v,dzu=dzu, dzv=dzv, & dzSxN=dzSxN, dzSyN=dzSyN, halo=1) - do j=js,je ; do i=is,ie - CS%sqg_struct(i,j,1) = 1.0 - enddo ; enddo - if (allocated(MEKE%Le)) then - do j=js,je ; do i=is,ie - Le(i,j) = MEKE%Le(i,j) - f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & - G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8/US%T_to_s) - enddo ; enddo + if (CS%sqg_expo<=0.) then + CS%sqg_struct(:,:,:) = 1. else do j=js,je ; do i=is,ie - Le(i,j) = sqrt(G%areaT(i,j)) - f(i,j) = max(0.25*abs(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & - G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)), 1.0e-8*US%T_to_s) + CS%sqg_struct(i,j,1) = 1.0 enddo ; enddo + if (allocated(MEKE%Le)) then + do j=js,je ; do i=is,ie + Le(i,j) = MEKE%Le(i,j) + f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround) + enddo ; enddo + else + do j=js,je ; do i=is,ie + Le(i,j) = sqrt(G%areaT(i,j)) + f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround) + enddo ; enddo + endif + do k=2,nz ; do j=js,je ; do i=is,ie + N2 = max(0.25 * ((N2_u(I-1,j,k) + N2_u(I,j,k)) + (N2_v(i,J-1,k) + N2_v(i,J,k))), 0.0) + dzc = 0.25 * ((dzu(I-1,j,k) + dzu(I,j,k)) + (dzv(i,J-1,k) + dzv(i,J,k))) + CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * & + exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j)))) + enddo ; enddo ; enddo endif - if (CS%debug) then - call hchksum(Le, 'SQG length scale', G%HI, unscale=US%L_to_m) - call hchksum(f, 'Coriolis at h point', G%HI, unscale=US%s_to_T) - call uvchksum( 'MEKE LmixScale', dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.) - endif - do k=2,nz ; do j=js,je ; do i=is,ie - N2 = max(0.25*(N2_u(I-1,j,k) + N2_u(I,j,k) + N2_v(i,J-1,k) + N2_v(i,J,k)),0.0) - dzc = 0.25*(dzu(I-1,j,k) + dzu(I,j,k) + dzv(i,J-1,k) + dzv(i,J,k)) * & - N2**0.5/f(i,j)*US%Z_to_L -! dzs = -N2**0.5/f(i,j)*dzc - CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1)*exp(-CS%sqg_expo*dzc/Le(i,j)) - enddo ; enddo ; enddo if (query_averaging_enabled(CS%diag)) then @@ -1418,7 +1399,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) 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, "BS_USE_SQG", CS%BS_use_sqg, & + call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", CS%BS_use_sqg_struct, & "If true, the SQG vertical structure is used for backscatter "//& "on the condition that BS_EBT_power=0", & default=.false.) @@ -1503,8 +1484,33 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) endif - allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0) - CS%BS_struct(:,:,:) = 1.0 + if (CS%BS_EBT_power>0. .and. CS%BS_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: BS_EBT_POWER>0. & + & and BS_USE_SQG=True cannot be set together") + endif + + if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT & + & and KHTH_USE_SQG_STRUCT can be true") + endif + + if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT & + & and KHTR_USE_SQG_STRUCT can be true") + endif + + if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) then + call MOM_error(FATAL, & + "calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT & + & and KD_GL90_USE_SQG_STRUCT can be true") + endif + + if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then + allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0) + endif if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then allocate(CS%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0) @@ -1615,13 +1621,15 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, & 'Vertical structure of SQG mode', 'nondim') - if (CS%BS_use_sqg .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & - .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then + if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct & + .or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0) endif - CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, & - 'Vertical structure of backscatter', 'nondim') + if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then + CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, & + 'Vertical structure of backscatter', 'nondim') + endif if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then CS%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, Time, & 'Vertical structure of thickness diffusivity', 'nondim') @@ -1888,7 +1896,7 @@ 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 & - .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct) + .or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) deallocate(CS%ebt_struct) if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct) if (allocated(CS%BS_struct)) deallocate(CS%BS_struct) if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) deallocate(CS%khth_struct)