From b72c633b5f100a0a8075df9b79380d9049cec5a4 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Mon, 3 Jun 2019 14:29:00 +0900
Subject: [PATCH] fsd defined for 0.1 <= c <= 0.9

---
 code/code_gen/terms.input |   4 +-
 code/m_calculate.f90      |  16 ++
 code/m_terms.f90          | 298 ++++++++++++++++++++------------------
 3 files changed, 178 insertions(+), 140 deletions(-)

diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index cbafc1e..3917170 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -2,7 +2,9 @@
 
 c = 1.0 - y
 wrate = rxn_rate(c)
-fsd = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c)))
+
+mask = $threshold_min_max(c)
+fsd = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c))) * mask
 
 sd = ((d2dx(c) + d2dy(c) + d2dz(c)) * $rod + wrate) / fsd
 
diff --git a/code/m_calculate.f90 b/code/m_calculate.f90
index a36bdaa..7f8577a 100644
--- a/code/m_calculate.f90
+++ b/code/m_calculate.f90
@@ -239,4 +239,20 @@ contains
 
      end function rxn_rate
 
+     real function threshold_min_max (c)
+
+      real :: c
+      real :: minc, maxc
+
+      minc = 0.1
+      maxc = 0.9
+
+        if ((c.lt.minc) .or. (c.gt.maxc)) then
+            threshold_min_max = 0.
+        else
+            threshold_min_max = 1.0
+        end if
+
+     end function threshold_min_max
+
 end module m_calculate
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 7a1dcb6..40d3dfa 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -194,49 +194,11 @@ 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))) ) ) )
+! ( threshold_min_max ( c(i,j,k) ) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * xyzbuffer2(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) * xyzbuffer2(i,j,k)
-end do
-end do
-end do
-
-
-! ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer1(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer2(i,j,k) )
-end do
-end do
-end do
-
-
-! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer0(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer2(i,j,k) )
+xyzbuffer2(i,j,k) = ( threshold_min_max ( xyzbuffer11(i,j,k) ) )
 end do
 end do
 end do
@@ -246,7 +208,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( rxn_rate ( xyzbuffer11(i,j,k) ) )
+xyzbuffer1(i,j,k) = ( rxn_rate ( xyzbuffer11(i,j,k) ) )
 end do
 end do
 end do
@@ -287,14 +249,13 @@ end do
 end do
 end do
 
-call d2dx ( xyzbuffer6, xyzbuffer11 )
-call ddx ( xyzbuffer5, xyzbuffer1 )
-call ddz ( xyzbuffer4, xyzbuffer0 )
+call d2dx ( xyzbuffer0, 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))) ) ) ) * mask(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer1(i,j,k) * xyzbuffer2(i,j,k)
+xyzbuffer7(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))) ) ) ) * xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -303,7 +264,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer0(i,j,k) * xyzbuffer2(i,j,k)
+fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * xyzbuffer7(i,j,k)
 end do
 end do
 end do
@@ -312,27 +273,27 @@ 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) * xyzbuffer2(i,j,k)
+fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * xyzbuffer7(i,j,k)
 end do
 end do
 end do
 
 
-! ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
+! ( ( - ddx_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) ) / xyzbuffer2(i,j,k) )
+xyzbuffer6(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer7(i,j,k) )
 end do
 end do
 end do
 
 
-! ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
+! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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) ) )
+xyzbuffer5(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer7(i,j,k) )
 end do
 end do
 end do
@@ -341,7 +302,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_d2dx_c(i) = avg_d2dx_c(i) + xyzbuffer6(i,j,k)
+avg_d2dx_c(i) = avg_d2dx_c(i) + xyzbuffer0(i,j,k)
 end do
 end do
 end do
@@ -355,12 +316,70 @@ end do
 end do
 end do
 
-call d2dz ( xyzbuffer0, xyzbuffer11 )
+call d2dz ( xyzbuffer4, xyzbuffer11 )
+call ddz ( xyzbuffer11, xyzbuffer5 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_ny(i) = fsd_avg_ny(i) + xyzbuffer9(i,j,k) * xyzbuffer2(i,j,k)
+fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer6(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_nz(i) = fsd_avg_nz(i) + xyzbuffer5(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_v(i) = fsd_avg_v(i) + v(i,j,k) * xyzbuffer7(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) )
+end do
+end do
+end do
+
+
+! ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
+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) ) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_d2dz_c(i) = avg_d2dz_c(i) + xyzbuffer4(i,j,k)
+end do
+end do
+end do
+
+call ddx ( xyzbuffer3, 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)
 end do
 end do
 end do
@@ -370,7 +389,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer11(i,j,k) = ( ( xyzbuffer6(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer0(i,j,k) )
+xyzbuffer2(i,j,k) = ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) )
 end do
 end do
 end do
@@ -379,26 +398,18 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_d2dz_c(i) = avg_d2dz_c(i) + xyzbuffer0(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) + xyzbuffer11(i,j,k)
+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
-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) )
+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
@@ -407,7 +418,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_fsd(i) = avg_fsd(i) + xyzbuffer2(i,j,k)
+avg_lapc(i) = avg_lapc(i) + xyzbuffer2(i,j,k)
 end do
 end do
 end do
@@ -416,7 +427,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer8(i,j,k) * xyzbuffer2(i,j,k)
+fsd_avg_sd(i) = fsd_avg_sd(i) + xyzbuffer9(i,j,k) * xyzbuffer7(i,j,k)
 end do
 end do
 end do
@@ -426,17 +437,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer0(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer1(i,j,k) )
-end do
-end do
-end do
-
-call ddy ( xyzbuffer7, xyzbuffer9 )
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer0(i,j,k) * xyzbuffer2(i,j,k)
+xyzbuffer2(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer9(i,j,k) ) * xyzbuffer6(i,j,k) )
 end do
 end do
 end do
@@ -446,7 +447,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer7(i,j,k) ) + xyzbuffer4(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer5(i,j,k) ) + xyzbuffer11(i,j,k) )
 end do
 end do
 end do
@@ -455,7 +456,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer6(i,j,k) * xyzbuffer2(i,j,k)
+fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer1(i,j,k) * xyzbuffer7(i,j,k)
 end do
 end do
 end do
@@ -465,7 +466,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( dabs ( xyzbuffer6(i,j,k) ) )
+xyzbuffer0(i,j,k) = ( dabs ( xyzbuffer1(i,j,k) ) )
 end do
 end do
 end do
@@ -474,7 +475,16 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer5(i,j,k) * xyzbuffer2(i,j,k)
+fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer2(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_absk(i) = fsd_avg_absk(i) + xyzbuffer0(i,j,k) * xyzbuffer7(i,j,k)
 end do
 end do
 end do
@@ -623,41 +633,11 @@ call ddx ( xyzbuffer8, xyzbuffer10 )
 call ddy ( xyzbuffer3, xyzbuffer10 )
 call ddz ( xyzbuffer2, xyzbuffer10 )
 
-! ( 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))) ) ) )
+! ( threshold_min_max ( 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))) ) ) )
-end do
-end do
-end do
-
-
-! ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer0(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer1(i,j,k) )
-end do
-end do
-end do
-
-
-! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( - xyzbuffer2(i,j,k) ) / xyzbuffer1(i,j,k) )
-end do
-end do
-end do
-
-
-! (((nx(i,j,k) - fsd_avg_nx(i)))*((nx(i,j,k) - fsd_avg_nx(i))))
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer6(i,j,k) = (((xyzbuffer0(i,j,k) - fsd_avg_nx(i)))*((xyzbuffer0(i,j,k) - fsd_avg_nx(i))))
+xyzbuffer1(i,j,k) = ( threshold_min_max ( xyzbuffer10(i,j,k) ) )
 end do
 end do
 end do
@@ -667,19 +647,59 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( rxn_rate ( xyzbuffer10(i,j,k) ) )
+xyzbuffer0(i,j,k) = ( rxn_rate ( xyzbuffer10(i,j,k) ) )
 end do
 end do
 end do
 
-call d2dx ( xyzbuffer4, xyzbuffer10 )
-call ddx ( xyzbuffer8, xyzbuffer0 )
-call ddz ( xyzbuffer2, xyzbuffer7 )
+call d2dx ( xyzbuffer7, xyzbuffer10 )
+
+! ( ( 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))) ) ) ) * mask(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer6(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))) ) ) ) * xyzbuffer1(i,j,k) )
+end do
+end do
+end do
+
+
+! ( ( - ddx_c(i,j,k) ) / fsd(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) )
+end do
+end do
+end do
+
+
+! ( ( - ddz_c(i,j,k) ) / fsd(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) )
+end do
+end do
+end do
+
+
+! (((nx(i,j,k) - fsd_avg_nx(i)))*((nx(i,j,k) - fsd_avg_nx(i))))
+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))))
+end do
+end do
+end do
+
+call d2dz ( xyzbuffer2, xyzbuffer10 )
+call ddz ( xyzbuffer1, xyzbuffer4 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer8(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -689,7 +709,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer1(i,j,k) )
+xyzbuffer10(i,j,k) = ( ( - xyzbuffer3(i,j,k) ) / xyzbuffer6(i,j,k) )
 end do
 end do
 end do
@@ -699,29 +719,29 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer3(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer0(i,j,k) ) + ( v(i,j,k) * xyzbuffer6(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer7(i,j,k) ) )
+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) ) )
 end do
 end do
 end do
 
-call d2dz ( xyzbuffer7, xyzbuffer10 )
+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
-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) )
+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, xyzbuffer6 )
+call ddy ( xyzbuffer9, xyzbuffer10 )
 
 ! ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer10(i,j,k) ) * xyzbuffer0(i,j,k) )
+xyzbuffer2(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer4(i,j,k) ) * xyzbuffer5(i,j,k) )
 end do
 end do
 end do
@@ -731,7 +751,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer6(i,j,k) = (((xyzbuffer7(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer7(i,j,k) - fsd_avg_t1(i))))
+xyzbuffer0(i,j,k) = (((xyzbuffer2(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer2(i,j,k) - fsd_avg_t1(i))))
 end do
 end do
 end do
@@ -741,7 +761,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer2(i,j,k) )
+xyzbuffer7(i,j,k) = ( ( xyzbuffer3(i,j,k) + xyzbuffer9(i,j,k) ) + xyzbuffer1(i,j,k) )
 end do
 end do
 end do
@@ -750,7 +770,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer6(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer0(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -760,7 +780,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer4(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer0(i,j,k) - fsd_avg_nx(i)) )
+xyzbuffer4(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_nx(i)) )
 end do
 end do
 end do
@@ -770,7 +790,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = (((xyzbuffer5(i,j,k) - fsd_avg_divn(i)))*((xyzbuffer5(i,j,k) - fsd_avg_divn(i))))
+xyzbuffer10(i,j,k) = (((xyzbuffer7(i,j,k) - fsd_avg_divn(i)))*((xyzbuffer7(i,j,k) - fsd_avg_divn(i))))
 end do
 end do
 end do
@@ -780,7 +800,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer5(i,j,k) ) )
+xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer7(i,j,k) ) )
 end do
 end do
 end do
@@ -789,7 +809,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) * xyzbuffer1(i,j,k)
+fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer10(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -799,7 +819,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_divn(i)) )
+xyzbuffer8(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer7(i,j,k) - fsd_avg_divn(i)) )
 end do
 end do
 end do
@@ -818,7 +838,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) * xyzbuffer1(i,j,k)
+fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer4(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -827,7 +847,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) * xyzbuffer1(i,j,k)
+fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer3(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -836,7 +856,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) * xyzbuffer1(i,j,k)
+fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer8(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do
@@ -846,7 +866,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( (xyzbuffer7(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
+xyzbuffer1(i,j,k) = ( (xyzbuffer2(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
 end do
 end do
 end do
@@ -855,7 +875,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer2(i,j,k) * xyzbuffer1(i,j,k)
+fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer1(i,j,k) * xyzbuffer6(i,j,k)
 end do
 end do
 end do