diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index 3917170..f784986 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -4,13 +4,14 @@ c = 1.0 - y
 wrate = rxn_rate(c)
 
 mask = $threshold_min_max(c)
-fsd = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c))) * mask
+fsd_orig = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c))) * mask
+fsd = fsd_orig * mask
 
-sd = ((d2dx(c) + d2dy(c) + d2dz(c)) * $rod + wrate) / fsd
+sd = ((d2dx(c) + d2dy(c) + d2dz(c)) * $rod + wrate) / fsd_orig
 
-nx = - ddx(c) / fsd
-ny = - ddy(c) / fsd
-nz = - ddz(c) / fsd
+nx = - ddx(c) / fsd_orig
+ny = - ddy(c) / fsd_orig
+nz = - ddz(c) / fsd_orig
 
 vn = u * nx + v * ny + w * nz
 
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 40d3dfa..6565208 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -261,25 +261,7 @@ end do
 end do
 
 
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * xyzbuffer7(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * xyzbuffer7(i,j,k)
-end do
-end do
-end do
-
-
-! ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddx_c(i,j,k) ) / fsd_orig(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
@@ -289,7 +271,7 @@ end do
 end do
 
 
-! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddz_c(i,j,k) ) / fsd_orig(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
@@ -319,10 +301,11 @@ end do
 call d2dz ( xyzbuffer4, xyzbuffer11 )
 call ddz ( xyzbuffer11, xyzbuffer5 )
 
+! ( fsd_orig(i,j,k) * mask(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer6(i,j,k) * xyzbuffer7(i,j,k)
+xyzbuffer9(i,j,k) = ( xyzbuffer7(i,j,k) * xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -331,7 +314,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer5(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer6(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -340,17 +323,35 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_v(i) = fsd_avg_v(i) + v(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer5(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
 
 
-! ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer9(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer7(i,j,k) )
+fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+
+! ( ( - ddy_c(i,j,k) ) / fsd_orig(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer3(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer7(i,j,k) )
 end do
 end do
 end do
@@ -360,7 +361,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer6(i,j,k) ) + ( v(i,j,k) * xyzbuffer9(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer5(i,j,k) ) )
+xyzbuffer2(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer6(i,j,k) ) + ( v(i,j,k) * xyzbuffer3(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer5(i,j,k) ) )
 end do
 end do
 end do
@@ -374,12 +375,50 @@ end do
 end do
 end do
 
-call ddx ( xyzbuffer3, xyzbuffer6 )
+call ddx ( xyzbuffer5, xyzbuffer6 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_ny(i) = fsd_avg_ny(i) + xyzbuffer9(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_ny(i) = fsd_avg_ny(i) + xyzbuffer3(i,j,k) * xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_v(i) = fsd_avg_v(i) + v(i,j,k) * xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+
+! ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd_orig(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer8(i,j,k) = ( ( ( ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) ) * rod ) + xyzbuffer1(i,j,k) ) / xyzbuffer7(i,j,k) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_fsd(i) = avg_fsd(i) + xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+call ddy ( xyzbuffer1, xyzbuffer3 )
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer8(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -389,7 +428,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) )
+xyzbuffer7(i,j,k) = ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) )
 end do
 end do
 end do
@@ -398,36 +437,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_fsd(i) = avg_fsd(i) + xyzbuffer7(i,j,k)
-end do
-end do
-end do
-
-call ddy ( xyzbuffer5, xyzbuffer9 )
-
-! ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer9(i,j,k) = ( ( ( ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) ) * rod ) + xyzbuffer1(i,j,k) ) / xyzbuffer7(i,j,k) )
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_lapc(i) = avg_lapc(i) + xyzbuffer2(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer9(i,j,k) * xyzbuffer7(i,j,k)
+avg_lapc(i) = avg_lapc(i) + xyzbuffer7(i,j,k)
 end do
 end do
 end do
@@ -437,7 +447,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer9(i,j,k) ) * xyzbuffer6(i,j,k) )
+xyzbuffer4(i,j,k) = ( ( xyzbuffer2(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer6(i,j,k) )
 end do
 end do
 end do
@@ -447,7 +457,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer1(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer5(i,j,k) ) + xyzbuffer11(i,j,k) )
+xyzbuffer10(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer1(i,j,k) ) + xyzbuffer11(i,j,k) )
 end do
 end do
 end do
@@ -456,7 +466,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer1(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer10(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -466,7 +476,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer0(i,j,k) = ( dabs ( xyzbuffer1(i,j,k) ) )
+xyzbuffer8(i,j,k) = ( dabs ( xyzbuffer10(i,j,k) ) )
 end do
 end do
 end do
@@ -475,7 +485,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer2(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer4(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -484,7 +494,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer0(i,j,k) * xyzbuffer7(i,j,k)
+fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer8(i,j,k) * xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -664,21 +674,43 @@ end do
 end do
 
 
-! ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddy_c(i,j,k) ) / fsd_orig(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer6(i,j,k) )
+xyzbuffer5(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer6(i,j,k) )
+end do
+end do
+end do
+
+call d2dz ( xyzbuffer4, xyzbuffer10 )
+call ddy ( xyzbuffer10, xyzbuffer5 )
+
+! ( fsd_orig(i,j,k) * mask(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer3(i,j,k) = ( xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k) )
 end do
 end do
 end do
 
 
-! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddx_c(i,j,k) ) / fsd_orig(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer4(i,j,k) = ( ( - xyzbuffer2(i,j,k) ) / xyzbuffer6(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer6(i,j,k) )
+end do
+end do
+end do
+
+
+! ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd_orig(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer8(i,j,k) = ( ( ( ( ( xyzbuffer7(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer4(i,j,k) ) * rod ) + xyzbuffer0(i,j,k) ) / xyzbuffer6(i,j,k) )
 end do
 end do
 end do
@@ -688,28 +720,27 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = (((xyzbuffer5(i,j,k) - fsd_avg_nx(i)))*((xyzbuffer5(i,j,k) - fsd_avg_nx(i))))
+xyzbuffer0(i,j,k) = (((xyzbuffer1(i,j,k) - fsd_avg_nx(i)))*((xyzbuffer1(i,j,k) - fsd_avg_nx(i))))
 end do
 end do
 end do
 
-call d2dz ( xyzbuffer2, xyzbuffer10 )
-call ddz ( xyzbuffer1, xyzbuffer4 )
+call ddx ( xyzbuffer7, xyzbuffer1 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer8(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer0(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
 
 
-! ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddz_c(i,j,k) ) / fsd_orig(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer6(i,j,k) )
+xyzbuffer4(i,j,k) = ( ( - xyzbuffer2(i,j,k) ) / xyzbuffer6(i,j,k) )
 end do
 end do
 end do
@@ -719,29 +750,18 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer5(i,j,k) ) + ( v(i,j,k) * xyzbuffer10(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer4(i,j,k) ) )
+xyzbuffer9(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer1(i,j,k) ) + ( v(i,j,k) * xyzbuffer5(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer4(i,j,k) ) )
 end do
 end do
 end do
 
-call ddx ( xyzbuffer3, xyzbuffer5 )
-
-! ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer4(i,j,k) = ( ( ( ( ( xyzbuffer7(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer2(i,j,k) ) * rod ) + xyzbuffer0(i,j,k) ) / xyzbuffer6(i,j,k) )
-end do
-end do
-end do
-
-call ddy ( xyzbuffer9, xyzbuffer10 )
+call ddz ( xyzbuffer2, xyzbuffer4 )
 
 ! ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer4(i,j,k) ) * xyzbuffer5(i,j,k) )
+xyzbuffer0(i,j,k) = ( ( xyzbuffer9(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer1(i,j,k) )
 end do
 end do
 end do
@@ -751,7 +771,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer0(i,j,k) = (((xyzbuffer2(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer2(i,j,k) - fsd_avg_t1(i))))
+xyzbuffer6(i,j,k) = (((xyzbuffer0(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer0(i,j,k) - fsd_avg_t1(i))))
 end do
 end do
 end do
@@ -761,7 +781,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer1(i,j,k) )
+xyzbuffer5(i,j,k) = ( ( xyzbuffer7(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -770,7 +790,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer0(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer6(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
@@ -780,7 +800,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer4(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_nx(i)) )
+xyzbuffer4(i,j,k) = ( (xyzbuffer0(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer1(i,j,k) - fsd_avg_nx(i)) )
 end do
 end do
 end do
@@ -790,7 +810,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = (((xyzbuffer7(i,j,k) - fsd_avg_divn(i)))*((xyzbuffer7(i,j,k) - fsd_avg_divn(i))))
+xyzbuffer10(i,j,k) = (((xyzbuffer5(i,j,k) - fsd_avg_divn(i)))*((xyzbuffer5(i,j,k) - fsd_avg_divn(i))))
 end do
 end do
 end do
@@ -800,7 +820,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer7(i,j,k) ) )
+xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer5(i,j,k) ) )
 end do
 end do
 end do
@@ -809,7 +829,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer10(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer10(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
@@ -819,7 +839,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer7(i,j,k) - fsd_avg_divn(i)) )
+xyzbuffer8(i,j,k) = ( (xyzbuffer0(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_divn(i)) )
 end do
 end do
 end do
@@ -829,7 +849,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer3(i,j,k) = (((xyzbuffer9(i,j,k) - fsd_avg_absk(i)))*((xyzbuffer9(i,j,k) - fsd_avg_absk(i))))
+xyzbuffer2(i,j,k) = (((xyzbuffer9(i,j,k) - fsd_avg_absk(i)))*((xyzbuffer9(i,j,k) - fsd_avg_absk(i))))
 end do
 end do
 end do
@@ -838,7 +858,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer4(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer4(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
@@ -847,7 +867,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer3(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer2(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
@@ -856,7 +876,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer8(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer8(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do
@@ -866,7 +886,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer1(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
+xyzbuffer1(i,j,k) = ( (xyzbuffer0(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
 end do
 end do
 end do
@@ -875,7 +895,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer1(i,j,k) * xyzbuffer6(i,j,k)
+fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer1(i,j,k) * xyzbuffer3(i,j,k)
 end do
 end do
 end do