Add a numerical mixing diagnostic to MOM6

src/CMakeLists.txt

@@ -274,6 +274,7 @@ target_sources(mom6shared PRIVATE
   tracer/MOM_tracer_diabatic.F90
   tracer/MOM_tracer_flow_control.F90
   tracer/MOM_tracer_hor_diff.F90
+  tracer/MOM_tracer_numerical_mixing.F90
   tracer/MOM_tracer_registry.F90
   tracer/MOM_tracer_types.F90
   tracer/MOM_tracer_Z_init.F90

src/diagnostics/MOM_diagnostics.F90

@@ -29,7 +29,7 @@ module MOM_diagnostics
 use MOM_interface_heights, only : find_eta, find_col_mass
 use MOM_spatial_means,     only : global_area_mean, global_layer_mean
 use MOM_spatial_means,     only : global_volume_mean, global_area_integral
-use MOM_tracer_registry,   only : tracer_registry_type, post_tracer_transport_diagnostics
+use MOM_tracer_registry,   only : tracer_type, tracer_registry_type, post_tracer_transport_diagnostics
 use MOM_unit_scaling,      only : unit_scale_type
 use MOM_variables,         only : thermo_var_ptrs, ocean_internal_state, p3d
 use MOM_variables,         only : accel_diag_ptrs, cont_diag_ptrs, surface
@@ -164,7 +164,6 @@ module MOM_diagnostics
   !>@}
 end type transport_diag_IDs
 
-
 contains
 !> Diagnostics not more naturally calculated elsewhere are computed here.
 subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
@@ -1654,7 +1653,7 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy
   real :: Idt             ! The inverse of the time interval [T-1 ~> s-1]
   real :: H_to_RZ_dt      ! A conversion factor from accumulated transports to fluxes
                           ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1].
-  integer :: i, j, k, is, ie, js, je, nz
+  integer :: i, j, k, is, ie, js, je, nz, m
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   Idt = 1. / dt_trans
@@ -1705,7 +1704,7 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy
     call post_data(IDs%id_dynamics_h_tendency, h_tend, diag, alt_h=diag_pre_dyn%h_state)
   endif
 
-  call post_tracer_transport_diagnostics(G, GV, Reg, diag_pre_dyn%h_state, diag)
+  call post_tracer_transport_diagnostics(G, GV, Reg, diag_pre_dyn%h_state, diag, uhtr, vhtr, h, dt_trans, Idt)
 
   call diag_restore_grids(diag)
 

src/tracer/MOM_tracer_numerical_mixing.F90

@@ -0,0 +1,160 @@
+!> Functions and routines involved in calculating numerical mixing of tracers due to advection schemes.
+module MOM_tracer_numerical_mixing
+
+use MOM_grid,          only : ocean_grid_type
+use MOM_tracer_types,  only : tracer_type
+use MOM_verticalGrid,  only : verticalGrid_type
+
+implicit none ; private
+
+#include <MOM_memory.h>
+
+public advection_scheme_variance_production
+
+contains
+
+!< Calculate the spurious variance production of tracer `Tr` due to the advection schemes.
+subroutine advection_scheme_variance_production(G, GV, Tr, h_diag, h, dt_trans, Idt, uhtr, vhtr, asvp)
+
+  type(ocean_grid_type),                        intent(in) :: G         !< Ocean grid structure
+  type(verticalGrid_type),                      intent(in) :: GV        !< Ocean vertical grid structure
+  type(tracer_type),                            intent(in) :: Tr        !< Pointer to the tracer regsitry
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in) :: h_diag    !< Previous thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in) :: h         !< Updated layer thicknesses [H ~> m or kg m-2]
+  real,                                         intent(in) :: dt_trans  !< The transport time interval [T ~> s]
+  real,                                         intent(in) :: Idt       !< Inverse time interval [T-1 ~> s-1]
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),   intent(in) :: uhtr      !< Accumulated zonal thickness fluxes
+                                                                        !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),   intent(in) :: vhtr      !< Accumulated meridional thickness fluxes
+                                                                        !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: asvp      !< Advection scheme varianve production
+                                                                        !! diagnostic [CU2 H T-1 ~> conc2 m s-1]
+
+  call thickness_weighted_variance_advection(G, GV, Tr, h_diag, h, dt_trans, Idt, asvp)
+  call thickness_weighted_variance_flux_divergence(G, Gv, Tr, uhtr, vhtr, Idt, asvp)
+
+end subroutine advection_scheme_variance_production
+
+!< Subroutine to calculate the thickness weighted variance advection over the transport timestep.
+subroutine thickness_weighted_variance_advection(G, GV, Tr, h_diag, h, dt, Idt, asvp)
+
+  type(ocean_grid_type),                        intent(in) :: G       !< Ocean grid structure
+  type(verticalGrid_type),                      intent(in) :: GV      !< Ocean vertical grid structure
+  type(tracer_type),                            intent(in) :: Tr      !< Pointer to the tracer registry
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in) :: h_diag  !< Previous thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),    intent(in) :: h       !< Updated thicknesses [H ~> m or kg m-2]
+  real,                                         intent(in) :: dt      !< Transport time interval [T ~> s]
+  real,                                         intent(in) :: Idt     !< Inverse time interval [T-1 ~> s-1]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: asvp    !< Array to store asvp in
+                                                                      !! [CU2 H T-1 ~> conc2 m s-1]
+
+  !< Local variables
+  integer :: is, ie, js, je, nz  !< Grid cell centre and layer indexes
+  integer :: i, j, k             !< Counters
+  real :: Ih                     !< Inverse updated thickness [H-1 ~> m-1]
+  real :: ht_prev                !< Thickness weighted tracer prior to dynamics [CU H ~> conc m]
+  real :: ht_adv                 !< Thickness weighted tracer after lateral advection [CU H ~> conc m]
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+  do k=1,nz ; do j=js,je ; do i=is,ie
+    Ih = 1 / h(i,j,k)
+    ht_prev = h_diag(i,j,k) * Tr%t_prev(i,j,k)
+    ht_adv = ht_prev + dt * Tr%advection_xy(i,j,k)
+    asvp(i,j,k) = ( (Ih * (ht_adv*ht_adv)) - (ht_prev * Tr%t_prev(i,j,k)) ) * Idt
+  enddo ; enddo ; enddo
+
+end subroutine thickness_weighted_variance_advection
+
+!< Subroutine to calculate the divergence of the variance flux due to the advection scheme
+subroutine thickness_weighted_variance_flux_divergence(G, Gv, Tr, uhtr, vhtr, Idt, asvp)
+
+  type(ocean_grid_type),                         intent(in) :: G     !< Ocean grid structure
+  type(verticalGrid_type),                       intent(in) :: GV    !< Ocean vertical grid structure
+  type(tracer_type),                             intent(in) :: Tr    !< Pointer to the tracer registry
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),    intent(in) :: uhtr  !< Accumulated zonal thickness fluxes
+                                                                     !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),    intent(in) :: vhtr  !< Accumulated meridional thickness fluxes
+                                                                     !! used to advect tracers [H L2 ~> m3 or kg]
+  real,                                          intent(in) :: Idt   !< Inverse time interval [T-1 ~> s-1]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: asvp  !< Array to store asvp in
+                                                                     !! [CU2 H T-1 ~> conc2 m s-1]
+
+  !< Local variables
+  integer :: is, ie, js, je, nz                           !< Grid cell centre and layer indexes
+  integer :: i, j, k                                      !< Counters
+  real :: east, west, north, south                        !< east, west, north and south faces of h point
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: x_upwind  !< Zonal upwind values for tracer [CU ~> conc]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: y_upwind  !< Meridional upwind values for tracer [CU ~> conc]
+
+  x_upwind(:,:,:) = 0.
+  y_upwind(:,:,:) = 0.
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+  call zonal_upwind_values(G, GV, Tr, uhtr, x_upwind)
+  call meridional_upwind_values(G, GV, Tr, vhtr, y_upwind)
+
+  do k=1,nz ;  do j=js,je ; do i=is,ie
+     east = (2 * (Tr%ad_x(I,j,k)  *x_upwind(I,j,k)))   - ((Idt*uhtr(I,j,k))   * (x_upwind(I,j,k)  *x_upwind(I,j,k)))
+     west = (2 * (Tr%ad_x(I-1,j,k)*x_upwind(I-1,j,k))) - ((Idt*uhtr(I-1,j,k)) * (x_upwind(I-1,j,k)*x_upwind(I-1,j,k)))
+    north = (2 * (Tr%ad_y(i,J,k)  *y_upwind(i,J,k)))   - ((Idt*vhtr(i,J,k))   * (y_upwind(i,J,k)  *y_upwind(i,J,k)))
+    south = (2 * (Tr%ad_y(i,J-1,k)*y_upwind(i,J-1,k))) - ((Idt*vhtr(i,J-1,k)) * (y_upwind(i,J-1,k)*y_upwind(i,J-1,k)))
+    asvp(i,j,k) = asvp(i,j,k) + (((east - west) + (north - south)) * G%IareaT(i,j))
+  enddo ; enddo; enddo
+
+end subroutine thickness_weighted_variance_flux_divergence
+
+!< Subroutine to calculate upwind values in zonal direction
+subroutine zonal_upwind_values(G, GV, Tr, uhtr, x_upwind)
+
+  type(ocean_grid_type),                         intent(in) :: G         !< Ocean grid structure
+  type(verticalGrid_type),                       intent(in) :: GV        !< Ocean vertical grid structure
+  type(tracer_type),                             intent(in) :: Tr        !< Pointer to the tracer registry
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),    intent(in) :: uhtr      !< Accumulated zonal thickness fluxes
+                                                                         !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: x_upwind  !< zonal upwind tracer [CU ~> conc]
+
+  !< Local variables
+  integer :: is, ie, js, je, nz  !< Grid cell centre indexes
+  integer :: i, j, k             !< Counters
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+  do k=1,nz ;  do j=js,je ; do I=is-1,ie
+    if (uhtr(I,j,k) >= 0) then
+      x_upwind(I,j,k) = Tr%t_prev(i,j,k)
+    elseif (uhtr(I,j,k) < 0) then
+      x_upwind(I,j,k) = Tr%t_prev(i+1,j,k)
+    endif
+  enddo ; enddo ; enddo
+
+end subroutine zonal_upwind_values
+
+!< Subroutine to calculate upwind values in the meridional direction
+subroutine meridional_upwind_values(G, GV, Tr, vhtr, y_upwind)
+
+  type(ocean_grid_type),                        intent(in) :: G         !< Ocean grid structure
+  type(verticalGrid_type),                      intent(in) :: GV        !< Ocean vertical grid structure
+  type(tracer_type),                            intent(in) :: Tr        !< Pointer to tracer registry
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),   intent(in) :: vhtr      !< Accumulated meridional thickness fluxes used
+                                                                        !! to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),intent(inout) :: y_upwind  !< Meridional upwind tracer values [CU ~> conc]
+
+  !< Local variables
+  integer :: is, ie, js, je, nz  !< Grid cell centre indexes
+  integer :: i, j, k             !< Counters
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+  do k=1,nz ; do J=js-1,je ; do i=is,ie
+    if (vhtr(i,J,k) >= 0) then
+      y_upwind(i,J,k) = Tr%t_prev(i,j,k)
+    elseif (vhtr(i,J,k) < 0) then
+      y_upwind(i,J,k) = Tr%t_prev(i,j+1,k)
+    endif
+  enddo ; enddo ; enddo
+
+end subroutine meridional_upwind_values
+
+end module MOM_tracer_numerical_mixing

src/tracer/MOM_tracer_registry.F90

@@ -14,6 +14,7 @@ module MOM_tracer_registry
 use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data, safe_alloc_ptr
 use MOM_diag_mediator, only : diag_grid_storage
 use MOM_diag_mediator, only : diag_copy_storage_to_diag, diag_save_grids, diag_restore_grids
+use MOM_domains,       only : pass_var
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe
 use MOM_file_parser,   only : get_param, log_version, param_file_type
 use MOM_hor_index,     only : hor_index_type
@@ -25,6 +26,7 @@ module MOM_tracer_registry
 use MOM_unit_scaling,  only : unit_scale_type
 use MOM_verticalGrid,  only : verticalGrid_type
 use MOM_tracer_types,  only : tracer_type, tracer_registry_type
+use MOM_tracer_numerical_mixing,  only : advection_scheme_variance_production
 
 implicit none ; private
 
@@ -382,6 +384,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
           "flux from the horizontal boundary diffusion scheme", trim(flux_units), &
           v_extensive=.true., &
           x_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T)
+      Tr%id_advection_scheme_variance_production = register_diag_field("ocean_model", &
+          trim(shortnm)//"_advection_scheme_variance_production", diag%axesTL, Time, &
+          "Spurious variance production of "//trim(shortnm)//" variance due to advection", &
+          trim(Tr%units)//"2 m s-1", conversion=(TR%conc_scale**2)*GV%H_to_MKS*US%s_to_T)
     else
       Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", &
           diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), &
@@ -407,6 +413,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
           "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), &
           flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, &
           x_cell_method='sum')
+      Tr%id_advection_scheme_variance_production = register_diag_field("ocean_model", &
+          trim(shortnm)//"_advection_scheme_variance_production", diag%axesTL, Time, &
+          "Spurious variance production of "//trim(shortnm)//" variance due to advection", &
+          trim(Tr%units)//"2 m s-1", conversion=(TR%conc_scale**2)*GV%H_to_MKS*US%s_to_T)
     endif
     Tr%id_zint = register_diag_field("ocean_model", trim(shortnm)//"_zint", &
         diag%axesT1, Time, &
@@ -418,8 +428,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
         trim(units) // " m")
     Tr%id_surf = register_diag_field("ocean_model", trim(shortnm)//"_SURF", &
         diag%axesT1, Time, "Surface values of "// trim(longname), trim(units))
-    if (Tr%id_adx > 0) call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz)
-    if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz)
+    if ((Tr%id_adx > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) &
+      call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz)
+    if ((Tr%id_ady > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) &
+      call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz)
     if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz)
     if (Tr%id_dfy > 0) call safe_alloc_ptr(Tr%df_y,isd,ied,JsdB,JedB,nz)
     if (Tr%id_hbd_dfx > 0) call safe_alloc_ptr(Tr%hbd_dfx,IsdB,IedB,jsd,jed,nz)
@@ -468,17 +480,17 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
         diag%axesT1, Time, &
         'Vertical sum of horizontal convergence of residual mean advective fluxes of '//&
         trim(lowercase(flux_longname)), conv_units, conversion=Tr%conv_scale*US%s_to_T)
-    if ((Tr%id_adv_xy > 0) .or. (Tr%id_adv_xy_2d > 0)) &
+    if ((Tr%id_adv_xy > 0) .or. (Tr%id_adv_xy_2d > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) &
       call safe_alloc_ptr(Tr%advection_xy,isd,ied,jsd,jed,nz)
 
     Tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', &
         diag%axesTL, Time, &
         'Net time tendency for '//trim(lowercase(longname)), &
         trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T)
 
-    if (Tr%id_tendency > 0) then
+    if ((Tr%id_tendency > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) then
       call safe_alloc_ptr(Tr%t_prev,isd,ied,jsd,jed,nz)
-      do k=1,nz ; do j=js,je ; do i=is,ie
+      do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
         Tr%t_prev(i,j,k) = Tr%t(i,j,k)
       enddo ; enddo ; enddo
     endif
@@ -733,6 +745,12 @@ subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt)
       enddo ; enddo ; enddo
       call post_data(Tr%id_tendency, work3d, diag, alt_h=diag_prev%h_state)
     endif
+    if (Tr%id_advection_scheme_variance_production > 0) then
+      call pass_var(Tr%t, G%Domain, halo=2)
+      do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
+        tr%t_prev(i,j,k) =  Tr%t(i,j,k)
+      enddo ; enddo ; enddo
+    endif
     if ((Tr%id_trxh_tendency > 0) .or. (Tr%id_trxh_tendency_2d > 0)) then
       do k=1,nz ; do j=js,je ; do i=is,ie
         work3d(i,j,k)     = (Tr%t(i,j,k)*h(i,j,k) - Tr%Trxh_prev(i,j,k)) * Idt
@@ -754,23 +772,38 @@ subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt)
 end subroutine post_tracer_diagnostics_at_sync
 
 !> Post the advective and diffusive tendencies
-subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag)
+subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag, uhtr, vhtr, h, dt_trans, Idt)
   type(ocean_grid_type),      intent(in) :: G    !< The ocean's grid structure
   type(verticalGrid_type),    intent(in) :: GV   !< The ocean's vertical grid structure
   type(tracer_registry_type), pointer    :: Reg  !< pointer to the tracer registry
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                               intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2]
   type(diag_ctrl),            intent(in) :: diag !< structure to regulate diagnostic output
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
+                              intent(in) :: uhtr !< Accumulated zonal thickness fluxes
+                                                 !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
+                              intent(in) :: vhtr !< Accumulated meridional thickness fluxes
+                                                 !! used to advect tracers [H L2 ~> m3 or kg]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                              intent(in) :: h   !< The updated layer thicknesses [H ~> m or kg m-2]
+  real,                       intent(in) :: dt_trans  !< The transport time interval [T ~> s]
+  real,                       intent(in) :: Idt       !< The inverse of the time interval [T-1 ~> s-1]
 
   integer :: i, j, k, is, ie, js, je, nz, m, khi
   real    :: work2d(SZI_(G),SZJ_(G))      ! The vertically integrated convergence of lateral advective
                                           ! tracer fluxes [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
   real    :: frac_under_100m(SZI_(G),SZJ_(G),SZK_(GV)) ! weights used to compute 100m vertical integrals [nondim]
   real    :: ztop(SZI_(G),SZJ_(G)) ! position of the top interface [H ~> m or kg m-2]
   real    :: zbot(SZI_(G),SZJ_(G)) ! position of the bottom interface [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV))   :: asvp ! Advection scheme varianve production of a
+                                                      ! tracer [CU2 H T-1 ~> conc2 m s-1]
+  real :: H_to_RZ_dt      ! A conversion factor from accumulated transports to fluxes
+                          ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1].
   type(tracer_type), pointer :: Tr=>NULL()
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+  H_to_RZ_dt = GV%H_to_RZ * Idt
 
   ! If any tracers are posting 100m vertical integrals, compute weights
   frac_under_100m(:,:,:) = 0.0
@@ -821,6 +854,12 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag)
       call post_data(Tr%id_adv_xy_2d, work2d, diag)
     endif
 
+    if (Tr%id_advection_scheme_variance_production > 0) then
+      asvp(:,:,:) = 0.
+      call advection_scheme_variance_production(G, GV, Tr, h_diag, h, dt_trans, Idt, uhtr, vhtr, asvp)
+      call post_data(Tr%id_advection_scheme_variance_production, asvp, diag, alt_h=h_diag)
+    endif
+
     ! A few diagnostics introduce with MARBL driver
     ! Compute full-depth vertical integral
     if (Tr%id_zint > 0) then