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

Author: Wendazhang33

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

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, &

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,14 +232,15 @@ 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]
   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
+  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)