From 091a9029654e12f3f5f23da92b5cb5c9fb7bf130 Mon Sep 17 00:00:00 2001
From: ignis <ignisization@gmail.com>
Date: Tue, 4 Jun 2019 12:54:22 +0900
Subject: [PATCH] bugfix

---
 code/code_gen/terms.input |   2 +-
 code/m_terms.f90          | 693 +++++++++++++++++++-------------------
 2 files changed, 349 insertions(+), 346 deletions(-)

diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index f784986..1f57d0f 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -4,7 +4,7 @@ c = 1.0 - y
 wrate = rxn_rate(c)
 
 mask = $threshold_min_max(c)
-fsd_orig = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c))) * mask
+fsd_orig = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c)))
 fsd = fsd_orig * mask
 
 sd = ((d2dx(c) + d2dy(c) + d2dz(c)) * $rod + wrate) / fsd_orig
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 6565208..52705d8 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -49,6 +49,7 @@ real*8, allocatable, dimension(:,:,:) :: xyzbuffer8
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer9
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer10
 real*8, allocatable, dimension(:,:,:) :: xyzbuffer11
+real*8, allocatable, dimension(:,:,:) :: xyzbuffer12
 
 contains
 
@@ -97,6 +98,7 @@ 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.
+allocate(xyzbuffer12(nxp,nyp,nzp), stat=ierr) ; xyzbuffer12 = 0.
 
 end subroutine m_terms_init
 
@@ -144,6 +146,7 @@ deallocate(xyzbuffer8)
 deallocate(xyzbuffer9)
 deallocate(xyzbuffer10)
 deallocate(xyzbuffer11)
+deallocate(xyzbuffer12)
 
 end subroutine m_terms_finalize
 
@@ -184,164 +187,31 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer11(i,j,k) = ( 1.0 - y(i,j,k) )
+xyzbuffer12(i,j,k) = ( 1.0 - y(i,j,k) )
+end do
+end do
+end do
+
+call d2dy ( xyzbuffer11, xyzbuffer12 )
+call ddx ( xyzbuffer10, xyzbuffer12 )
+call ddy ( xyzbuffer9, xyzbuffer12 )
+call ddz ( xyzbuffer8, xyzbuffer12 )
+
+! ( 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
+xyzbuffer3(i,j,k) = ( sqrt ( ( ( ((xyzbuffer10(i,j,k))*(xyzbuffer10(i,j,k))) + ((xyzbuffer9(i,j,k))*(xyzbuffer9(i,j,k))) ) + ((xyzbuffer8(i,j,k))*(xyzbuffer8(i,j,k))) ) ) )
 end do
 end do
 end do
 
-call d2dy ( xyzbuffer10, xyzbuffer11 )
-call ddx ( xyzbuffer9, xyzbuffer11 )
-call ddy ( xyzbuffer8, xyzbuffer11 )
-call ddz ( xyzbuffer3, xyzbuffer11 )
 
 ! ( threshold_min_max ( c(i,j,k) ) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = ( threshold_min_max ( xyzbuffer11(i,j,k) ) )
-end do
-end do
-end do
-
-
-! ( rxn_rate ( c(i,j,k) ) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer1(i,j,k) = ( rxn_rate ( xyzbuffer11(i,j,k) ) )
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_c(i) = avg_c(i) + xyzbuffer11(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_d2dy_c(i) = avg_d2dy_c(i) + xyzbuffer10(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_ddx_c(i) = avg_ddx_c(i) + xyzbuffer9(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_ddz_c(i) = avg_ddz_c(i) + xyzbuffer3(i,j,k)
-end do
-end do
-end do
-
-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
-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
-
-
-! ( ( - ddx_c(i,j,k) ) / fsd_orig(i,j,k) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer6(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer7(i,j,k) )
-end do
-end do
-end do
-
-
-! ( ( - ddz_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) = ( ( - xyzbuffer3(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_d2dx_c(i) = avg_d2dx_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_ddy_c(i) = avg_ddy_c(i) + xyzbuffer8(i,j,k)
-end do
-end do
-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
-xyzbuffer9(i,j,k) = ( xyzbuffer7(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_nx(i) = fsd_avg_nx(i) + xyzbuffer6(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_nz(i) = fsd_avg_nz(i) + xyzbuffer5(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_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)
+xyzbuffer2(i,j,k) = ( threshold_min_max ( xyzbuffer12(i,j,k) ) )
 end do
 end do
 end do
@@ -351,7 +221,149 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer3(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer7(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer3(i,j,k) )
+end do
+end do
+end do
+
+
+! ( rxn_rate ( c(i,j,k) ) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer0(i,j,k) = ( rxn_rate ( xyzbuffer12(i,j,k) ) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_c(i) = avg_c(i) + xyzbuffer12(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_d2dy_c(i) = avg_d2dy_c(i) + xyzbuffer11(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_ddx_c(i) = avg_ddx_c(i) + xyzbuffer10(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_ddz_c(i) = avg_ddz_c(i) + xyzbuffer8(i,j,k)
+end do
+end do
+end do
+
+call d2dx ( xyzbuffer7, xyzbuffer12 )
+call ddy ( xyzbuffer6, xyzbuffer1 )
+
+! ( fsd_orig(i,j,k) * mask(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer5(i,j,k) = ( xyzbuffer3(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_ny(i) = fsd_avg_ny(i) + xyzbuffer1(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
+fsd_avg_u(i) = fsd_avg_u(i) + u(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
+fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * xyzbuffer5(i,j,k)
+end do
+end do
+end do
+
+
+! ( ( - 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) = ( ( - xyzbuffer10(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
+avg_d2dx_c(i) = avg_d2dx_c(i) + xyzbuffer7(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_ddy_c(i) = avg_ddy_c(i) + xyzbuffer9(i,j,k)
+end do
+end do
+end do
+
+call d2dz ( xyzbuffer10, xyzbuffer12 )
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_nx(i) = fsd_avg_nx(i) + xyzbuffer4(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
+fsd_avg_v(i) = fsd_avg_v(i) + v(i,j,k) * xyzbuffer5(i,j,k)
+end do
+end do
+end do
+
+
+! ( ( - ddz_c(i,j,k) ) / fsd_orig(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer9(i,j,k) = ( ( - xyzbuffer8(i,j,k) ) / xyzbuffer3(i,j,k) )
 end do
 end do
 end do
@@ -361,7 +373,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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) ) )
+xyzbuffer8(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer4(i,j,k) ) + ( v(i,j,k) * xyzbuffer1(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer9(i,j,k) ) )
 end do
 end do
 end do
@@ -370,55 +382,17 @@ 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)
+avg_d2dz_c(i) = avg_d2dz_c(i) + xyzbuffer10(i,j,k)
 end do
 end do
 end do
 
-call ddx ( xyzbuffer5, xyzbuffer6 )
+call ddx ( xyzbuffer2, xyzbuffer4 )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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)
+fsd_avg_nz(i) = fsd_avg_nz(i) + xyzbuffer9(i,j,k) * xyzbuffer5(i,j,k)
 end do
 end do
 end do
@@ -428,7 +402,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer7(i,j,k) = ( ( xyzbuffer0(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer4(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( xyzbuffer7(i,j,k) + xyzbuffer11(i,j,k) ) + xyzbuffer10(i,j,k) )
 end do
 end do
 end do
@@ -437,7 +411,36 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_lapc(i) = avg_lapc(i) + xyzbuffer7(i,j,k)
+avg_fsd(i) = avg_fsd(i) + xyzbuffer5(i,j,k)
+end do
+end do
+end do
+
+call ddz ( xyzbuffer12, xyzbuffer9 )
+
+! ( ( ( ( ( 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
+xyzbuffer9(i,j,k) = ( ( ( ( ( xyzbuffer7(i,j,k) + xyzbuffer11(i,j,k) ) + xyzbuffer10(i,j,k) ) * rod ) + xyzbuffer0(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
+avg_lapc(i) = avg_lapc(i) + xyzbuffer1(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) * xyzbuffer5(i,j,k)
 end do
 end do
 end do
@@ -447,7 +450,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer4(i,j,k) = ( ( xyzbuffer2(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer6(i,j,k) )
+xyzbuffer3(i,j,k) = ( ( xyzbuffer8(i,j,k) + xyzbuffer9(i,j,k) ) * xyzbuffer4(i,j,k) )
 end do
 end do
 end do
@@ -457,7 +460,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer1(i,j,k) ) + xyzbuffer11(i,j,k) )
+xyzbuffer1(i,j,k) = ( ( xyzbuffer2(i,j,k) + xyzbuffer6(i,j,k) ) + xyzbuffer12(i,j,k) )
 end do
 end do
 end do
@@ -466,7 +469,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer10(i,j,k) * xyzbuffer9(i,j,k)
+fsd_avg_divn(i) = fsd_avg_divn(i) + xyzbuffer1(i,j,k) * xyzbuffer5(i,j,k)
 end do
 end do
 end do
@@ -476,7 +479,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( dabs ( xyzbuffer10(i,j,k) ) )
+xyzbuffer0(i,j,k) = ( dabs ( xyzbuffer1(i,j,k) ) )
 end do
 end do
 end do
@@ -485,7 +488,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer4(i,j,k) * xyzbuffer9(i,j,k)
+fsd_avg_t1(i) = fsd_avg_t1(i) + xyzbuffer3(i,j,k) * xyzbuffer5(i,j,k)
 end do
 end do
 end do
@@ -494,7 +497,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) * xyzbuffer9(i,j,k)
+fsd_avg_absk(i) = fsd_avg_absk(i) + xyzbuffer0(i,j,k) * xyzbuffer5(i,j,k)
 end do
 end do
 end do
@@ -633,64 +636,21 @@ integer :: i, j, k
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer10(i,j,k) = ( 1.0 - y(i,j,k) )
+xyzbuffer12(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 ( xyzbuffer11, xyzbuffer12 )
+call ddx ( xyzbuffer10, xyzbuffer12 )
+call ddy ( xyzbuffer9, xyzbuffer12 )
+call ddz ( xyzbuffer8, xyzbuffer12 )
 
-! ( threshold_min_max ( c(i,j,k) ) )
+! ( 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) = ( threshold_min_max ( xyzbuffer10(i,j,k) ) )
-end do
-end do
-end do
-
-
-! ( rxn_rate ( c(i,j,k) ) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer0(i,j,k) = ( rxn_rate ( xyzbuffer10(i,j,k) ) )
-end do
-end do
-end do
-
-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
-
-
-! ( ( - 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) = ( ( - 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) )
+xyzbuffer3(i,j,k) = ( sqrt ( ( ( ((xyzbuffer10(i,j,k))*(xyzbuffer10(i,j,k))) + ((xyzbuffer9(i,j,k))*(xyzbuffer9(i,j,k))) ) + ((xyzbuffer8(i,j,k))*(xyzbuffer8(i,j,k))) ) ) )
 end do
 end do
 end do
@@ -700,37 +660,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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
-
-
-! (((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
-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 ddx ( xyzbuffer7, xyzbuffer1 )
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer0(i,j,k) * xyzbuffer3(i,j,k)
+xyzbuffer2(i,j,k) = ( ( - xyzbuffer10(i,j,k) ) / xyzbuffer3(i,j,k) )
 end do
 end do
 end do
@@ -740,7 +670,70 @@ end do
 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) ) / xyzbuffer3(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
+xyzbuffer0(i,j,k) = (((xyzbuffer2(i,j,k) - fsd_avg_nx(i)))*((xyzbuffer2(i,j,k) - fsd_avg_nx(i))))
+end do
+end do
+end do
+
+
+! ( rxn_rate ( c(i,j,k) ) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer7(i,j,k) = ( rxn_rate ( xyzbuffer12(i,j,k) ) )
+end do
+end do
+end do
+
+call d2dx ( xyzbuffer6, xyzbuffer12 )
+call ddx ( xyzbuffer5, xyzbuffer2 )
+call ddz ( xyzbuffer4, xyzbuffer1 )
+
+! ( threshold_min_max ( c(i,j,k) ) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer10(i,j,k) = ( threshold_min_max ( xyzbuffer12(i,j,k) ) )
+end do
+end do
+end do
+
+call d2dz ( xyzbuffer8, xyzbuffer12 )
+
+! ( fsd_orig(i,j,k) * mask(i,j,k) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer12(i,j,k) = ( xyzbuffer3(i,j,k) * xyzbuffer10(i,j,k) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_t6(i) = fsd_avg_t6(i) + xyzbuffer0(i,j,k) * xyzbuffer12(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
+xyzbuffer10(i,j,k) = ( ( - xyzbuffer9(i,j,k) ) / xyzbuffer3(i,j,k) )
 end do
 end do
 end do
@@ -750,28 +743,18 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-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) ) )
+xyzbuffer9(i,j,k) = ( ( ( u(i,j,k) * xyzbuffer2(i,j,k) ) + ( v(i,j,k) * xyzbuffer10(i,j,k) ) ) + ( w(i,j,k) * xyzbuffer1(i,j,k) ) )
 end do
 end do
 end do
 
-call ddz ( xyzbuffer2, xyzbuffer4 )
+call ddy ( xyzbuffer1, xyzbuffer10 )
 
-! ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
+! ( ( ( ( ( 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
-xyzbuffer0(i,j,k) = ( ( xyzbuffer9(i,j,k) + xyzbuffer8(i,j,k) ) * xyzbuffer1(i,j,k) )
-end do
-end do
-end do
-
-
-! (((t1(i,j,k) - fsd_avg_t1(i)))*((t1(i,j,k) - fsd_avg_t1(i))))
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer6(i,j,k) = (((xyzbuffer0(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer0(i,j,k) - fsd_avg_t1(i))))
+xyzbuffer0(i,j,k) = ( ( ( ( ( xyzbuffer6(i,j,k) + xyzbuffer11(i,j,k) ) + xyzbuffer8(i,j,k) ) * rod ) + xyzbuffer7(i,j,k) ) / xyzbuffer3(i,j,k) )
 end do
 end do
 end do
@@ -781,55 +764,17 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer5(i,j,k) = ( ( xyzbuffer7(i,j,k) + xyzbuffer10(i,j,k) ) + xyzbuffer2(i,j,k) )
+xyzbuffer7(i,j,k) = ( ( xyzbuffer5(i,j,k) + xyzbuffer1(i,j,k) ) + xyzbuffer4(i,j,k) )
 end do
 end do
 end do
 
 
+! ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer6(i,j,k) * xyzbuffer3(i,j,k)
-end do
-end do
-end do
-
-
-! ( (t1(i,j,k) - fsd_avg_t1(i)) * (nx(i,j,k) - fsd_avg_nx(i)) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-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
-
-
-! (((divn(i,j,k) - fsd_avg_divn(i)))*((divn(i,j,k) - fsd_avg_divn(i))))
-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))))
-end do
-end do
-end do
-
-
-! ( dabs ( divn(i,j,k) ) )
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-xyzbuffer9(i,j,k) = ( dabs ( xyzbuffer5(i,j,k) ) )
-end do
-end do
-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) * xyzbuffer3(i,j,k)
+xyzbuffer6(i,j,k) = ( ( xyzbuffer9(i,j,k) + xyzbuffer0(i,j,k) ) * xyzbuffer2(i,j,k) )
 end do
 end do
 end do
@@ -839,17 +784,27 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer8(i,j,k) = ( (xyzbuffer0(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer5(i,j,k) - fsd_avg_divn(i)) )
+xyzbuffer5(i,j,k) = ( (xyzbuffer6(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer7(i,j,k) - fsd_avg_divn(i)) )
 end do
 end do
 end do
 
 
-! (((absk(i,j,k) - fsd_avg_absk(i)))*((absk(i,j,k) - fsd_avg_absk(i))))
+! (((t1(i,j,k) - fsd_avg_t1(i)))*((t1(i,j,k) - fsd_avg_t1(i))))
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer2(i,j,k) = (((xyzbuffer9(i,j,k) - fsd_avg_absk(i)))*((xyzbuffer9(i,j,k) - fsd_avg_absk(i))))
+xyzbuffer4(i,j,k) = (((xyzbuffer6(i,j,k) - fsd_avg_t1(i)))*((xyzbuffer6(i,j,k) - fsd_avg_t1(i))))
+end do
+end do
+end do
+
+
+! ( dabs ( divn(i,j,k) ) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer11(i,j,k) = ( dabs ( xyzbuffer7(i,j,k) ) )
 end do
 end do
 end do
@@ -858,7 +813,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) * xyzbuffer3(i,j,k)
+fsd_avg_t3(i) = fsd_avg_t3(i) + xyzbuffer5(i,j,k) * xyzbuffer12(i,j,k)
 end do
 end do
 end do
@@ -867,7 +822,27 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_avg_t8(i) = fsd_avg_t8(i) + xyzbuffer2(i,j,k) * xyzbuffer3(i,j,k)
+fsd_avg_t5(i) = fsd_avg_t5(i) + xyzbuffer4(i,j,k) * xyzbuffer12(i,j,k)
+end do
+end do
+end do
+
+
+! ( (t1(i,j,k) - fsd_avg_t1(i)) * (nx(i,j,k) - fsd_avg_nx(i)) )
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer10(i,j,k) = ( (xyzbuffer6(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer2(i,j,k) - fsd_avg_nx(i)) )
+end do
+end do
+end do
+
+
+! (((divn(i,j,k) - fsd_avg_divn(i)))*((divn(i,j,k) - fsd_avg_divn(i))))
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer9(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
@@ -876,7 +851,16 @@ 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) * xyzbuffer3(i,j,k)
+fsd_avg_t2(i) = fsd_avg_t2(i) + xyzbuffer10(i,j,k) * xyzbuffer12(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+fsd_avg_t7(i) = fsd_avg_t7(i) + xyzbuffer9(i,j,k) * xyzbuffer12(i,j,k)
 end do
 end do
 end do
@@ -886,7 +870,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-xyzbuffer1(i,j,k) = ( (xyzbuffer0(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer9(i,j,k) - fsd_avg_absk(i)) )
+xyzbuffer8(i,j,k) = ( (xyzbuffer6(i,j,k) - fsd_avg_t1(i)) * (xyzbuffer11(i,j,k) - fsd_avg_absk(i)) )
 end do
 end do
 end do
@@ -895,7 +879,26 @@ 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) * xyzbuffer3(i,j,k)
+fsd_avg_t4(i) = fsd_avg_t4(i) + xyzbuffer8(i,j,k) * xyzbuffer12(i,j,k)
+end do
+end do
+end do
+
+
+! (((absk(i,j,k) - fsd_avg_absk(i)))*((absk(i,j,k) - fsd_avg_absk(i))))
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+xyzbuffer3(i,j,k) = (((xyzbuffer11(i,j,k) - fsd_avg_absk(i)))*((xyzbuffer11(i,j,k) - fsd_avg_absk(i))))
+end do
+end do
+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) * xyzbuffer12(i,j,k)
 end do
 end do
 end do