diff --git a/code/code_gen/post.py b/code/code_gen/post.py
index 487d3bc..def389a 100644
--- a/code/code_gen/post.py
+++ b/code/code_gen/post.py
@@ -582,7 +582,7 @@ class Stage4():
         img = np.zeros((len(l1), len(l1)))
         for i, v in enumerate(l1):
             for j in range(i, len(l1)):
-                img[i,i:j] = img[i,i:j] + (1 if v in g[l1[j]] else 0)
+                img[i,i:j+1] = img[i,i:j+1] + (1 if v in g[l1[j]] else 0)
         return img > 0
 
 
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 07c7cfa..7a1dcb6 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -48,6 +48,7 @@ real*8, allocatable, dimension(:,:,:) :: xyzbuffer7
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer8
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer9
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer10
+real*8, allocatable, dimension(:,:,:) :: xyzbuffer11
 
 contains
 
@@ -95,6 +96,7 @@ allocate(xyzbuffer7(nxp,nyp,nzp), stat=ierr) ; xyzbuffer7 = 0.
 allocate(xyzbuffer8(nxp,nyp,nzp), stat=ierr) ; xyzbuffer8 = 0.
 allocate(xyzbuffer9(nxp,nyp,nzp), stat=ierr) ; xyzbuffer9 = 0.
 allocate(xyzbuffer10(nxp,nyp,nzp), stat=ierr) ; xyzbuffer10 = 0.
+allocate(xyzbuffer11(nxp,nyp,nzp), stat=ierr) ; xyzbuffer11 = 0.
 
 end subroutine m_terms_init
 
@@ -141,6 +143,7 @@ deallocate(xyzbuffer7)
 deallocate(xyzbuffer8)
 deallocate(xyzbuffer9)
 deallocate(xyzbuffer10)
+deallocate(xyzbuffer11)
 
 end subroutine m_terms_finalize
 
@@ -181,21 +184,21 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = ( 1.0 - y(i,j,k) )
+xyzbuffer11(i,j,k) = ( 1.0 - y(i,j,k) )
 end do
 end do
 end do
 
-call d2dy ( xyzbuffer9, xyzbuffer10 )
-call ddx ( xyzbuffer8, xyzbuffer10 )
-call ddy ( xyzbuffer3, xyzbuffer10 )
-call ddz ( xyzbuffer2, xyzbuffer10 )
+call d2dy ( xyzbuffer10, xyzbuffer11 )
+call ddx ( xyzbuffer9, xyzbuffer11 )
+call ddy ( xyzbuffer8, xyzbuffer11 )
+call ddz ( xyzbuffer3, xyzbuffer11 )
 
 ! ( sqrt ( ( ( ((ddx_c(i,j,k))*(ddx_c(i,j,k))) + ((ddy_c(i,j,k))*(ddy_c(i,j,k))) ) + ((ddz_c(i,j,k))*(ddz_c(i,j,k))) ) ) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer1(i,j,k) = ( sqrt ( ( ( ((xyzbuffer8(i,j,k))*(xyzbuffer8(i,j,k))) + ((xyzbuffer3(i,j,k))*(xyzbuffer3(i,j,k))) ) + ((xyzbuffer2(i,j,k))*(xyzbuffer2(i,j,k))) ) ) )
+xyzbuffer2(i,j,k) = ( sqrt ( ( ( ((xyzbuffer9(i,j,k))*(xyzbuffer9(i,j,k))) + ((xyzbuffer8(i,j,k))*(xyzbuffer8(i,j,k))) ) + ((xyzbuffer3(i,j,k))*(xyzbuffer3(i,j,k))) ) ) )
 end do
 end do
 end do
@@ -204,7 +207,7 @@ 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) * xyzbuffer1(i,j,k)
+fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -213,7 +216,7 @@ 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) * xyzbuffer1(i,j,k)
+fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -223,7 +226,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer0(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -233,7 +236,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( - xyzbuffer2(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer0(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -243,7 +246,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = ( rxn_rate ( xyzbuffer10(i,j,k) ) )
+xyzbuffer7(i,j,k) = ( rxn_rate ( xyzbuffer11(i,j,k) ) )
 end do
 end do
 end do
@@ -252,7 +255,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_c(i) = avg_c(i) + xyzbuffer10(i,j,k)
+avg_c(i) = avg_c(i) + xyzbuffer11(i,j,k)
 end do
 end do
 end do
@@ -261,7 +264,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_d2dy_c(i) = avg_d2dy_c(i) + xyzbuffer9(i,j,k)
+avg_d2dy_c(i) = avg_d2dy_c(i) + xyzbuffer10(i,j,k)
 end do
 end do
 end do
@@ -270,7 +273,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_ddx_c(i) = avg_ddx_c(i) + xyzbuffer8(i,j,k)
+avg_ddx_c(i) = avg_ddx_c(i) + xyzbuffer9(i,j,k)
 end do
 end do
 end do
@@ -279,19 +282,19 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_ddz_c(i) = avg_ddz_c(i) + xyzbuffer2(i,j,k)
+avg_ddz_c(i) = avg_ddz_c(i) + xyzbuffer3(i,j,k)
 end do
 end do
 end do
 
-call d2dx ( xyzbuffer5, xyzbuffer10 )
-call ddx ( xyzbuffer4, xyzbuffer0 )
-call ddz ( xyzbuffer8, xyzbuffer7 )
+call d2dx ( xyzbuffer6, xyzbuffer11 )
+call ddx ( xyzbuffer5, xyzbuffer1 )
+call ddz ( xyzbuffer4, xyzbuffer0 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer0(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer1(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -300,7 +303,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer7(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer0(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -309,7 +312,7 @@ 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) * xyzbuffer1(i,j,k)
+fsd_avg_v(i) = fsd_avg_v(i) + v(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -319,7 +322,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer9(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -329,7 +332,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer0(i,j,k) ) + ( v(i,j,k) * xyzbuffer2(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer7(i,j,k) ) )
+xyzbuffer3(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer1(i,j,k) ) + ( v(i,j,k) * xyzbuffer9(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer0(i,j,k) ) )
 end do
 end do
 end do
@@ -338,7 +341,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_d2dx_c(i) = avg_d2dx_c(i) + xyzbuffer5(i,j,k)
+avg_d2dx_c(i) = avg_d2dx_c(i) + xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -347,17 +350,17 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_ddy_c(i) = avg_ddy_c(i) + xyzbuffer3(i,j,k)
+avg_ddy_c(i) = avg_ddy_c(i) + xyzbuffer8(i,j,k)
 end do
 end do
 end do
 
-call d2dz ( xyzbuffer10, xyzbuffer10 )
+call d2dz ( xyzbuffer0, xyzbuffer11 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_ny(i) = fsd_avg_ny(i) + xyzbuffer2(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_ny(i) = fsd_avg_ny(i) + xyzbuffer9(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -367,7 +370,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer3(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer10(i,j,k) )
+xyzbuffer11(i,j,k) = ( ( xyzbuffer6(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer0(i,j,k) )
 end do
 end do
 end do
@@ -376,7 +379,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_d2dz_c(i) = avg_d2dz_c(i) + xyzbuffer10(i,j,k)
+avg_d2dz_c(i) = avg_d2dz_c(i) + xyzbuffer0(i,j,k)
 end do
 end do
 end do
@@ -385,7 +388,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_lapc(i) = avg_lapc(i) + xyzbuffer3(i,j,k)
+avg_lapc(i) = avg_lapc(i) + xyzbuffer11(i,j,k)
 end do
 end do
 end do
@@ -395,7 +398,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = ( ( ( ( ( xyzbuffer5(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer10(i,j,k) ) * rod ) + xyzbuffer6(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer8(i,j,k) = ( ( ( ( ( xyzbuffer6(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer0(i,j,k) ) * rod ) + xyzbuffer7(i,j,k) ) / xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -404,7 +407,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_fsd(i) = avg_fsd(i) + xyzbuffer1(i,j,k)
+avg_fsd(i) = avg_fsd(i) + xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -413,7 +416,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer8(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -423,17 +426,17 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( ( xyzbuffer7(i,j,k) + xyzbuffer6(i,j,k) ) * xyzbuffer0(i,j,k) )
+xyzbuffer0(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer1(i,j,k) )
 end do
 end do
 end do
 
-call ddy ( xyzbuffer10, xyzbuffer2 )
+call ddy ( xyzbuffer7, xyzbuffer9 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer5(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer0(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -443,7 +446,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer9(i,j,k) = ( ( xyzbuffer4(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer8(i,j,k) )
+xyzbuffer6(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer7(i,j,k) ) + xyzbuffer4(i,j,k) )
 end do
 end do
 end do
@@ -452,7 +455,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer9(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer6(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -462,7 +465,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( dabs ( xyzbuffer9(i,j,k) ) )
+xyzbuffer5(i,j,k) = ( dabs ( xyzbuffer6(i,j,k) ) )
 end do
 end do
 end do
@@ -471,7 +474,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer8(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer5(i,j,k) * xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -707,18 +710,18 @@ call d2dz ( xyzbuffer7, xyzbuffer10 )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( ( ( ( ( xyzbuffer4(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer7(i,j,k) ) * rod ) + xyzbuffer5(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer10(i,j,k) = ( ( ( ( ( xyzbuffer4(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer7(i,j,k) ) * rod ) + xyzbuffer5(i,j,k) ) / xyzbuffer1(i,j,k) )
 end do
 end do
 end do
 
-call ddy ( xyzbuffer4, xyzbuffer6 )
+call ddy ( xyzbuffer9, xyzbuffer6 )
 
 ! ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer5(i,j,k) ) * xyzbuffer0(i,j,k) )
+xyzbuffer7(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer10(i,j,k) ) * xyzbuffer0(i,j,k) )
 end do
 end do
 end do
@@ -728,7 +731,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer9(i,j,k) = (((xyzbuffer10(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer10(i,j,k) - fsd_avg_t1(i))))
+xyzbuffer6(i,j,k) = (((xyzbuffer7(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer7(i,j,k) - fsd_avg_t1(i))))
 end do
 end do
 end do
@@ -738,7 +741,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer4(i,j,k) ) + xyzbuffer2(i,j,k) )
+xyzbuffer5(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -747,7 +750,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer9(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do
@@ -757,7 +760,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer3(i,j,k) = ( (xyzbuffer10(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer0(i,j,k) - fsd_avg_nx(i)) )
+xyzbuffer4(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer0(i,j,k) - fsd_avg_nx(i)) )
 end do
 end do
 end do
@@ -767,7 +770,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = (((xyzbuffer8(i,j,k) - fsd_avg_divn(i)))*((xyzbuffer8(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
@@ -777,7 +780,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer0(i,j,k) = ( dabs ( xyzbuffer8(i,j,k) ) )
+xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer5(i,j,k) ) )
 end do
 end do
 end do
@@ -786,7 +789,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer2(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer10(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do
@@ -796,7 +799,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( (xyzbuffer10(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer8(i,j,k) - fsd_avg_divn(i)) )
+xyzbuffer8(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_divn(i)) )
 end do
 end do
 end do
@@ -806,7 +809,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = (((xyzbuffer0(i,j,k) - fsd_avg_absk(i)))*((xyzbuffer0(i,j,k) - fsd_avg_absk(i))))
+xyzbuffer3(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
@@ -815,7 +818,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer3(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer4(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do
@@ -824,7 +827,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer3(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do
@@ -833,7 +836,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer7(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer8(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do
@@ -843,7 +846,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( (xyzbuffer10(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer0(i,j,k) - fsd_avg_absk(i)) )
+xyzbuffer2(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
 end do
 end do
 end do
@@ -852,7 +855,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer5(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer2(i,j,k) * xyzbuffer1(i,j,k)
 end do
 end do
 end do