diff --git a/ROMS/Nonlinear/set_avg.F b/ROMS/Nonlinear/set_avg.F
index 7886b4fe..3a1c23b9 100755
--- a/ROMS/Nonlinear/set_avg.F
+++ b/ROMS/Nonlinear/set_avg.F
@@ -3865,6 +3865,7 @@ SUBROUTINE set_avg_tile (ng, tile,                                &
      &                          EWperiodic(ng), NSperiodic(ng),         &
      &                          AVERAGE(ng)%avgAKr)
 #    endif
+            END IF
 #   endif
 #  endif
 
diff --git a/ROMS/Nonlinear/step3d_t.F b/ROMS/Nonlinear/step3d_t.F
index e4ed0002..b5c99213 100755
--- a/ROMS/Nonlinear/step3d_t.F
+++ b/ROMS/Nonlinear/step3d_t.F
@@ -1973,23 +1973,34 @@ SUBROUTINE step3d_t_tile (ng, tile,                               &
       END DO
 # ifdef TS_VAR
 !
-!  Compute tracer variance.
+!  Compute numerical and physical destruction of tracer variance.
 !
       itrc=inert(1)
       DO j=JstrR,JendR
         DO i=IstrR,IendR
           DO k=1,N(ng)
-#  if !defined SPLINES_VDIFF || defined TS_MPDATA
+# if !defined SPLINES_VDIFF
+            ! Multiply by oHz to convert to Tunits.
             t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*oHz(i,j,k)
             t(i,j,k,nnew,itrc+1)=t(i,j,k,nnew,itrc+1)*oHz(i,j,k)
-#  endif
+# else
+            ! Multiply by oHz to convert to Tunits. Only necessary for MPDATA. 
+            if ((Hadvection(itrc,ng)%MPDATA)) then
+              t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*oHz(i,j,k)
+              t(i,j,k,nnew,itrc+1)=t(i,j,k,nnew,itrc+1)*oHz(i,j,k)
+            end if
+# endif
+            ! Compute numerical mixing: mnum = ( A{c^2} - (A{c})^2 ) / dt
             t(i,j,k,nnew,itrc+2)=(t(i,j,k,nnew,itrc+1)-                 &
      &                           t(i,j,k,nnew,itrc)*t(i,j,k,nnew,itrc))/&
      &                           dt(ng)
           END DO
           DO k=1,N(ng)-1
+            ! Compute vertical salinity gradient 
             cff=(t(i,j,k+1,nnew,isalt)-t(i,j,k,nnew,isalt))/            &
      &                         (z_r(i,j,k+1)-z_r(i,j,k))
+            ! Compute destruction of vertical salinity variance, 
+            ! mphy = 2*kappa_v*(ds/dz)^2
             AKr(i,j,k)=2.0_r8*cff*cff*(Akt(i,j,k,itemp)-                &
      &                                 Akt_bak(itemp,ng))
           END DO