diff --git a/code/code_gen/code_gen.py b/code/code_gen/code_gen.py
index ddcbf88..5fdfa30 100644
--- a/code/code_gen/code_gen.py
+++ b/code/code_gen/code_gen.py
@@ -95,6 +95,24 @@ avg_array_divide  = "{0} = {0} / denum {1}"
 
 real_array_diff  = "call {0[0]} ( {0[0]}_{0[1]}, {0[1]} )"
 
+class FortranCode:
+    def __init__ (self, exp):
+        self.exp = exp
+
+    def __repr__ (self):
+        return self.exp
+
+    '''
+    def __add__ (self, other):
+        self.exp = "( {} + {} )".format(self.exp, other.exp)
+        return self
+
+    def __sub__ (self, other):
+        self.exp = "( {} - {} )".format(self.exp, other.exp)
+        return self
+    '''
+
+
 @v_args(inline=True)    # Affects the signatures of the methods
 class ToFortran(Transformer):
 
@@ -102,29 +120,38 @@ class ToFortran(Transformer):
         self.primary = primary_set
         self.derivatives = {}
         self.dependency = {}
+        self.codes = {}
 
     def number(self, numeral):
-        return (str(float(numeral)), [])
+        return (FortranCode(str(float(numeral))), [])
 
     def env(self, name):
-        return (name.value, [])
+        return (FortranCode(name.value), [])
 
     def var(self, name):
-        return (name + "(i,j,k)",
+        return (FortranCode(name + "(i,j,k)"),
                 [name.value] if name.value not in self.primary else [])
 
     def fluc(self, name):
-        fmt = "({0}(i,j,k) - {{}}avg_{0}(i))"
-        return (fmt.format(name),
+        fmt = "({0}(i,j,k) - {{0}}avg_{0}(i))"
+        return (FortranCode(fmt.format(name)),
                 [name.value] if name.value not in self.primary else [])
 
     def dnx (self, partial, b):
         signature = "{}_{}".format(partial.data, b)
+        fcode = FortranCode(signature + "(i,j,k)")
 
         self.derivatives[signature] = (partial.data, b.value)
         self.dependency[signature] = [b.value]
+        self.codes[signature] = {
+            "decl" : real_array_decl.format(signature),
+            "alloc" : real_array_alloc.format(signature),
+            "free" : real_array_free.format(signature),
+            "calc" : real_array_diff.format((partial.data, b)),
+            "avg" : "",
+        }
 
-        return (signature + "(i,j,k)", [signature])
+        return (fcode, [signature])
 
     def icall (self, a, (b, dep)):
         fcode = "({0})".format(b)
@@ -190,17 +217,20 @@ class CalculateTree(Transformer):
         self.primary = []
         self.derived = {}
         self.averaged = {}
+        self.averages = []
         self.derivatives = {}
         self.dependency = {}
+        self.fluctuation = {}
         self.definitions = {}
         self.exp_parser = ToFortran([])
+        self.codes = {}
 
     def varlist(self, *args):
 
         for arg in args:
             self.primary.append(arg.value)
             self.dependency[arg.value] = []
-        # return self.primary
+            self.fluctuation[arg.value] = False
 
         return ""
 
@@ -211,6 +241,15 @@ class CalculateTree(Transformer):
 
         code, dep = self.exp_parser.transform(vdef)
         self.dependency[vname.value] = dep
+        self.fluctuation[vname.value] = "{0}" in code
+
+        self.codes[vname.value] = {
+            "decl" : real_array_decl.format(vname.value),
+            "alloc" : real_array_alloc.format(vname.value),
+            "free" : real_array_free.format(vname.value),
+            "calc" : real_array_loop.format((vname.value, code)),
+            "avg" : "",
+        }
 
         return ""
 
@@ -223,24 +262,55 @@ class CalculateTree(Transformer):
         w = str(weight)
         for v in vlist:
             avg_var = ( "" if w == str(None) else w + "_" ) + "avg_" + v
-            self.dependency[avg_var] = [str(v)]
-            # average_names.append(avg_var)
-            # self.depsDict[avg_var] = [v]
-            # self.flucDict[avg_var] = False
+            self.averages.append(avg_var)
+            self.fluctuation[avg_var] = False
 
-        '''
-        fmt = "avg_{}"
-        if weight is not None:
-            wvalue, wdep = self.var(weight)
-            self.averaged[fmt.format(weight)] = wvalue, None
-            self.dependency[fmt.format(weight)] = wdep
-            fmt = weight + "_" + fmt
+            var = str(v)
+
+            if self.fluctuation[str(v)]:
+                if w == str(None):
+                    var = var + "_"
+                    self.dependency[avg_var] =  [str(v)+"_"]
+
+                    self.dependency[str(v)+"_"] = self.dependency[str(v)]
+                    self.fluctuation[str(v)+"_"] = True
+
+                    self.codes[var] = {}
+                    for k,val in self.codes[str(v)].items():
+                        self.codes[var][k] = val.replace(str(v), var)
+
+                    self.codes[var]["calc"] = self.codes[var]["calc"].format("")
+
+                else:
+                    var = var + "_" + w
+                    self.dependency[avg_var] =  [str(v)+"_"+w, w]
+
+                    self.dependency[str(v)+"_"+w] = self.dependency[str(v)]
+                    self.fluctuation[str(v)+"_"+w] = True
+
+                    self.codes[var] = {}
+                    for k,val in self.codes[str(v)].items():
+                        self.codes[var][k] = val.replace(str(v), var)
+
+                    self.codes[var]["calc"] = self.codes[var]["calc"].format(w+"_")
+            else:
+                self.dependency[avg_var] = [str(v)] + ( [] if w == str(None) else [w] )
+
+            wfmt = "* {}(i,j,k)"
+            pWeight = (wfmt.format(w) if w != str(None) else "")
+
+            meanw = "/ avg_{}"
+            dWeight = (meanw.format(w) if w != str(None) else "")
+
+            self.codes[avg_var] = {
+            "decl" : avg_array_decl.format(avg_var),
+            "alloc" : avg_array_alloc.format(avg_var),
+            "free" : real_array_free.format(avg_var),
+            "calc" : avg_array_sum.format(avg_var, var+"(i,j,k)", pWeight),
+            "avg" : avg_array_divide.format(avg_var, dWeight)
+            }
 
-        for i, (value, dep) in enumerate(args):
-            self.averaged[fmt.format(i)] = value, weight
-            self.dependency[fmt.format(i)] = dep + (wdep if weight is not None else [])
 
-        '''
         return ""
 
 
@@ -256,6 +326,12 @@ class CalculateTree(Transformer):
                 + self.exp_parser.dependency.items()
                 )
 
+    def has_fluc (self):
+        return dict(
+                self.fluctuation.items()
+                + [(k, False) for k, v in self.exp_parser.dependency.items()]
+                )
+
     def array_decl (self):
 
         f_code = ""
@@ -368,8 +444,6 @@ class CalculateTree(Transformer):
 
 tf=CalculateTree()
 
-ft=ToFortran(['u','v','w','y'])
-
 
 calc_parser = Lark(calc_grammar, parser='lalr' ) # , transformer=tf)
 calc = calc_parser.parse
@@ -395,18 +469,15 @@ class FortranProgram:
 
 
         dg = self.parser.dep_graph()
-        fd = {}
+        fd = self.parser.has_fluc()
+
+        '''
         for v in dg.iterkeys():
             fd[v] = False
 
         for n, d in (self.parser.definitions.iteritems()):
             fd[n] = CheckPass.check(d)
-
-        average_names = []
-        for w, vlist in self.parser.averaged.iteritems():
-            for v in vlist:
-                avg_var = ( "" if w == "None" else w + "_" ) + "avg_" + v
-                average_names.append(avg_var)
+        '''
 
         def isFluc (a):
             for x in dg[a]:
@@ -414,6 +485,12 @@ class FortranProgram:
 
             return fd[a]
 
+        average_names = self.parser.averages
+
+        self.pass1avg = filter(lambda x: not isFluc(x), average_names)
+
+        self.pass2avg = filter(isFluc, average_names)
+
         pass1var = filter(lambda x: not isFluc(x), average_names)
 
         pass2var = filter(isFluc, average_names)
@@ -435,14 +512,85 @@ class FortranProgram:
 
             return c
 
-        self.pass1set = dep_closer(set(pass1var))
+        self.pass1set = dep_closer(set(pass1var)) - set(self.parser.primary)
 
-        self.pass2set = dep_closer(set(pass2var))
+        self.pass2set = dep_closer(set(pass2var)) - set(self.parser.primary)
 
-        print self.pass1set
+        self.pass1list = self.sort_vars(dg, self.pass1set)
 
-        print self.pass2set
+        self.pass2list = self.sort_vars(dg, self.pass2set)
 
+        self.codes = dict(self.parser.codes.items() + self.parser.exp_parser.codes.items())
+
+
+    def sort_vars (self, dependency, group):
+        order = []
+        remain = set(group)
+
+        while len(remain) > 0:
+            for v in remain:
+                if len(set(dependency[v]) & remain) == 0:
+                    order.append(v)
+                    remain.remove(v)
+                    break
+
+        return order
+
+    def print_program (self):
+
+        import StringIO
+
+        output = StringIO.StringIO()
+        output.write('First line.\n')
+        print >>output, 'Second line.'
+
+        decl = StringIO.StringIO()
+        alloc = StringIO.StringIO()
+        free = StringIO.StringIO()
+        calc1 = StringIO.StringIO()
+        avg1 = StringIO.StringIO()
+        calc2 = StringIO.StringIO()
+        avg2 = StringIO.StringIO()
+
+        for var in self.pass1set | self.pass2set:
+            print >>decl, self.codes[var]["decl"]
+
+        for var in self.pass1set | self.pass2set:
+            print >>alloc, self.codes[var]["alloc"]
+
+        for var in self.pass1set | self.pass2set:
+            print >>free, self.codes[var]["free"]
+
+        for var in self.pass1list :
+            print >>calc1, self.codes[var]["calc"]
+
+        for var in self.pass1avg :
+            print >>avg1, self.codes[var]["avg"]
+
+        for var in self.pass2list:
+            print >>calc2, self.codes[var]["calc"]
+
+        for var in self.pass2avg:
+            print >>avg2, self.codes[var]["avg"]
+
+        md = {}
+        md["module_name"] = "terms"
+        md["module_data"] = decl.getvalue()
+        md["module_init"] = alloc.getvalue()
+        md["module_finalize"] = free.getvalue()
+        md["module_pass1"] = calc1.getvalue()
+        md["module_pass1_avg"] = avg1.getvalue()
+        md["module_pass2"] = calc2.getvalue()
+        md["module_pass2_avg"] = avg2.getvalue()
+
+        return md
+
+
+    def calculate_pass1 (self):
+        return
+
+    def calculate_pass2 (self):
+        return
 
 
 
@@ -464,39 +612,11 @@ def test():
         terms_raw = ((inputfile.read()))
         parsed_tree = (calc(terms_raw))
 
-    '''
-    tf.transform(parsed_tree)
 
-    namelist, deflist = zip(*list(tf.definitions.iteritems()))
-
-    hasFluc = dict([ (n, hf) for n, hf in zip(namelist, map(CheckPass.check, deflist)) ])
-
-    codes, deps = zip(*map ( ft.transform, deflist ))
-
-    depsDict = dict(zip(namelist, deps))
-    
-    visited = {n: False for n in namelist}
-    
-    def isFluc (a, graph, visit, hf):
-        try:
-            if visit[a]:
-                pass
-            elif len(graph[a]) < 1:
-                visit[a] = True
-            else:
-                hf[a] = any([hf[a]] + [isFluc(x, graph, visit, hf) for x in graph[a]])
-                visit[a] = True
-
-            return hf[a]
-        except KeyError:
-            return False
-
-    flucDict = {n: isFluc(n, depsDict, visited, hasFluc) for n in namelist}
-    '''
+    fp = FortranProgram(terms_raw)
 
 
-    FortranProgram(terms_raw)
-
+    print mod_form.format( fp.print_program() )
 
     '''
     for f, ts in zip(namelist, deps):
diff --git a/code/code_gen/resources/m_template.f90 b/code/code_gen/resources/m_template.f90
index fb49c64..8d6e143 100644
--- a/code/code_gen/resources/m_template.f90
+++ b/code/code_gen/resources/m_template.f90
@@ -47,4 +47,25 @@ denum=real(nfiles*nyp*nzp)
 end subroutine m_{0[module_name]}_average_pass1
 
 
+subroutine m_{0[module_name]}_calculate_pass2
+
+integer :: i, j, k
+
+{0[module_pass2]}
+
+end subroutine m_{0[module_name]}_calculate_pass2
+
+
+subroutine m_{0[module_name]}_average_pass2 (nfiles)
+
+integer :: nfiles
+real*8 :: denum
+
+denum=real(nfiles*nyp*nzp)
+
+{0[module_pass2_avg]}
+
+end subroutine m_{0[module_name]}_average_pass2
+
+
 end module m_{0[module_name]}
diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index 42e4e8b..1eafc9b 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -8,12 +8,10 @@ k = (sqr(u')+sqr(v')+sqr(w'))/2.0
 
 tflux = u' * c_auto'
 
-avg { u, v, w, c_auto, tflux }
+avg { u, v, w, c_auto, tflux, y, fsd_auto}
 
 avg c_auto { u, v, w }
 
-avg y { u, v, w }
+avg y { u, v, w, tflux, c_auto }
 
 avg fsd_auto { u }
-
-avg temp {  }
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 1237c1e..97a2702 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -6,25 +6,29 @@ use m_calculate
 
 implicit none
 
-real*8, allocatable, dimension(:,:,:) :: c_auto
+real*8, allocatable, dimension(:) :: fsd_auto_avg_u
 real*8, allocatable, dimension(:,:,:) :: fsd_auto
-real*8, allocatable, dimension(:,:,:) :: wrate_auto
-real*8, allocatable, dimension(:,:,:) :: ddz_c_auto
-real*8, allocatable, dimension(:,:,:) :: ddy_c_auto
-real*8, allocatable, dimension(:,:,:) :: ddx_c_auto
-real*8, allocatable, dimension(:) :: avg_fsd_auto
-real*8, allocatable, dimension(:) :: fsd_auto_avg_0
 real*8, allocatable, dimension(:) :: avg_c_auto
-real*8, allocatable, dimension(:) :: c_auto_avg_2
-real*8, allocatable, dimension(:) :: y_avg_2
-real*8, allocatable, dimension(:) :: c_auto_avg_1
-real*8, allocatable, dimension(:) :: y_avg_0
-real*8, allocatable, dimension(:) :: y_avg_1
-real*8, allocatable, dimension(:) :: avg_0
-real*8, allocatable, dimension(:) :: avg_1
-real*8, allocatable, dimension(:) :: avg_2
-real*8, allocatable, dimension(:) :: c_auto_avg_0
+real*8, allocatable, dimension(:,:,:) :: ddz_c_auto
+real*8, allocatable, dimension(:,:,:) :: tflux_y
+real*8, allocatable, dimension(:) :: c_auto_avg_w
+real*8, allocatable, dimension(:) :: c_auto_avg_v
+real*8, allocatable, dimension(:) :: c_auto_avg_u
+real*8, allocatable, dimension(:) :: y_avg_w
+real*8, allocatable, dimension(:) :: avg_fsd_auto
+real*8, allocatable, dimension(:) :: y_avg_v
+real*8, allocatable, dimension(:,:,:) :: c_auto
+real*8, allocatable, dimension(:) :: y_avg_u
+real*8, allocatable, dimension(:,:,:) :: ddy_c_auto
+real*8, allocatable, dimension(:) :: y_avg_c_auto
+real*8, allocatable, dimension(:) :: avg_u
+real*8, allocatable, dimension(:) :: avg_v
+real*8, allocatable, dimension(:) :: avg_w
 real*8, allocatable, dimension(:) :: avg_y
+real*8, allocatable, dimension(:,:,:) :: tflux_
+real*8, allocatable, dimension(:) :: y_avg_tflux
+real*8, allocatable, dimension(:,:,:) :: ddx_c_auto
+real*8, allocatable, dimension(:) :: avg_tflux
 
 
 contains
@@ -33,25 +37,29 @@ subroutine m_terms_init
 
 integer :: ierr
 
-allocate(c_auto(nxp,nyp,nzp), stat=ierr) ; c_auto = 0.
+allocate(fsd_auto_avg_u(nxp), stat=ierr) ; fsd_auto_avg_u = 0.
 allocate(fsd_auto(nxp,nyp,nzp), stat=ierr) ; fsd_auto = 0.
-allocate(wrate_auto(nxp,nyp,nzp), stat=ierr) ; wrate_auto = 0.
-allocate(ddz_c_auto(nxp,nyp,nzp), stat=ierr) ; ddz_c_auto = 0.
-allocate(ddy_c_auto(nxp,nyp,nzp), stat=ierr) ; ddy_c_auto = 0.
-allocate(ddx_c_auto(nxp,nyp,nzp), stat=ierr) ; ddx_c_auto = 0.
-allocate(avg_fsd_auto(nxp), stat=ierr) ; avg_fsd_auto = 0.
-allocate(fsd_auto_avg_0(nxp), stat=ierr) ; fsd_auto_avg_0 = 0.
 allocate(avg_c_auto(nxp), stat=ierr) ; avg_c_auto = 0.
-allocate(c_auto_avg_2(nxp), stat=ierr) ; c_auto_avg_2 = 0.
-allocate(y_avg_2(nxp), stat=ierr) ; y_avg_2 = 0.
-allocate(c_auto_avg_1(nxp), stat=ierr) ; c_auto_avg_1 = 0.
-allocate(y_avg_0(nxp), stat=ierr) ; y_avg_0 = 0.
-allocate(y_avg_1(nxp), stat=ierr) ; y_avg_1 = 0.
-allocate(avg_0(nxp), stat=ierr) ; avg_0 = 0.
-allocate(avg_1(nxp), stat=ierr) ; avg_1 = 0.
-allocate(avg_2(nxp), stat=ierr) ; avg_2 = 0.
-allocate(c_auto_avg_0(nxp), stat=ierr) ; c_auto_avg_0 = 0.
+allocate(ddz_c_auto(nxp,nyp,nzp), stat=ierr) ; ddz_c_auto = 0.
+allocate(tflux_y(nxp,nyp,nzp), stat=ierr) ; tflux_y = 0.
+allocate(c_auto_avg_w(nxp), stat=ierr) ; c_auto_avg_w = 0.
+allocate(c_auto_avg_v(nxp), stat=ierr) ; c_auto_avg_v = 0.
+allocate(c_auto_avg_u(nxp), stat=ierr) ; c_auto_avg_u = 0.
+allocate(y_avg_w(nxp), stat=ierr) ; y_avg_w = 0.
+allocate(avg_fsd_auto(nxp), stat=ierr) ; avg_fsd_auto = 0.
+allocate(y_avg_v(nxp), stat=ierr) ; y_avg_v = 0.
+allocate(c_auto(nxp,nyp,nzp), stat=ierr) ; c_auto = 0.
+allocate(y_avg_u(nxp), stat=ierr) ; y_avg_u = 0.
+allocate(ddy_c_auto(nxp,nyp,nzp), stat=ierr) ; ddy_c_auto = 0.
+allocate(y_avg_c_auto(nxp), stat=ierr) ; y_avg_c_auto = 0.
+allocate(avg_u(nxp), stat=ierr) ; avg_u = 0.
+allocate(avg_v(nxp), stat=ierr) ; avg_v = 0.
+allocate(avg_w(nxp), stat=ierr) ; avg_w = 0.
 allocate(avg_y(nxp), stat=ierr) ; avg_y = 0.
+allocate(tflux_(nxp,nyp,nzp), stat=ierr) ; tflux_ = 0.
+allocate(y_avg_tflux(nxp), stat=ierr) ; y_avg_tflux = 0.
+allocate(ddx_c_auto(nxp,nyp,nzp), stat=ierr) ; ddx_c_auto = 0.
+allocate(avg_tflux(nxp), stat=ierr) ; avg_tflux = 0.
 
 
 end subroutine m_terms_init
@@ -59,25 +67,29 @@ end subroutine m_terms_init
 
 subroutine m_terms_finalize
 
-deallocate(c_auto)
+deallocate(fsd_auto_avg_u)
 deallocate(fsd_auto)
-deallocate(wrate_auto)
-deallocate(ddz_c_auto)
-deallocate(ddy_c_auto)
-deallocate(ddx_c_auto)
-deallocate(avg_fsd_auto)
-deallocate(fsd_auto_avg_0)
 deallocate(avg_c_auto)
-deallocate(c_auto_avg_2)
-deallocate(y_avg_2)
-deallocate(c_auto_avg_1)
-deallocate(y_avg_0)
-deallocate(y_avg_1)
-deallocate(avg_0)
-deallocate(avg_1)
-deallocate(avg_2)
-deallocate(c_auto_avg_0)
+deallocate(ddz_c_auto)
+deallocate(tflux_y)
+deallocate(c_auto_avg_w)
+deallocate(c_auto_avg_v)
+deallocate(c_auto_avg_u)
+deallocate(y_avg_w)
+deallocate(avg_fsd_auto)
+deallocate(y_avg_v)
+deallocate(c_auto)
+deallocate(y_avg_u)
+deallocate(ddy_c_auto)
+deallocate(y_avg_c_auto)
+deallocate(avg_u)
+deallocate(avg_v)
+deallocate(avg_w)
 deallocate(avg_y)
+deallocate(tflux_)
+deallocate(y_avg_tflux)
+deallocate(ddx_c_auto)
+deallocate(avg_tflux)
 
 
 end subroutine m_terms_finalize
@@ -91,7 +103,16 @@ integer :: i, j, k
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_1(i) = avg_1(i) + v(i,j,k) 
+y_avg_w(i) = y_avg_w(i) + w(i,j,k) * y(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+y_avg_v(i) = y_avg_v(i) + v(i,j,k) * y(i,j,k)
 end do
 end do
 end do
@@ -119,7 +140,7 @@ call ddz ( ddz_c_auto, c_auto )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-c_auto_avg_2(i) = c_auto_avg_2(i) + w(i,j,k) * c_auto(i,j,k)
+c_auto_avg_w(i) = c_auto_avg_w(i) + w(i,j,k) * c_auto(i,j,k)
 end do
 end do
 end do
@@ -128,7 +149,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-c_auto_avg_1(i) = c_auto_avg_1(i) + v(i,j,k) * c_auto(i,j,k)
+c_auto_avg_v(i) = c_auto_avg_v(i) + v(i,j,k) * c_auto(i,j,k)
 end do
 end do
 end do
@@ -137,7 +158,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-c_auto_avg_0(i) = c_auto_avg_0(i) + u(i,j,k) * c_auto(i,j,k)
+c_auto_avg_u(i) = c_auto_avg_u(i) + u(i,j,k) * c_auto(i,j,k)
 end do
 end do
 end do
@@ -146,34 +167,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-y_avg_2(i) = y_avg_2(i) + w(i,j,k) * y(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-y_avg_0(i) = y_avg_0(i) + u(i,j,k) * y(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-y_avg_1(i) = y_avg_1(i) + v(i,j,k) * y(i,j,k)
-end do
-end do
-end do
-
-
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-avg_0(i) = avg_0(i) + u(i,j,k) 
+y_avg_u(i) = y_avg_u(i) + u(i,j,k) * y(i,j,k)
 end do
 end do
 end do
@@ -183,7 +177,34 @@ call ddy ( ddy_c_auto, c_auto )
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-avg_2(i) = avg_2(i) + w(i,j,k) 
+y_avg_c_auto(i) = y_avg_c_auto(i) + c_auto(i,j,k) * y(i,j,k)
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_u(i) = avg_u(i) + u(i,j,k) 
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_v(i) = avg_v(i) + v(i,j,k) 
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_w(i) = avg_w(i) + w(i,j,k) 
 end do
 end do
 end do
@@ -211,7 +232,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_auto_avg_0(i) = fsd_auto_avg_0(i) + u(i,j,k) * fsd_auto(i,j,k)
+fsd_auto_avg_u(i) = fsd_auto_avg_u(i) + u(i,j,k) * fsd_auto(i,j,k)
 end do
 end do
 end do
@@ -226,15 +247,6 @@ end do
 end do
 
 
-do k = 1, nzp
-do j = 1, nyp
-do i = 1, nxp
-wrate_auto(i,j,k) = ( rxn_rate ( c_auto(i,j,k) ) )
-end do
-end do
-end do
-
-
 
 end subroutine m_terms_calculate_pass1
 
@@ -246,27 +258,92 @@ real*8 :: denum
 
 denum=real(nfiles*nyp*nzp)
 
-avg_y = avg_y / denum 
-avg_2 = avg_2 / denum 
-avg_1 = avg_1 / denum 
-avg_0 = avg_0 / denum 
+avg_u = avg_u / denum 
+avg_v = avg_v / denum 
+avg_w = avg_w / denum 
 avg_c_auto = avg_c_auto / denum 
+avg_y = avg_y / denum 
 avg_fsd_auto = avg_fsd_auto / denum 
-fsd_auto_avg_0 = fsd_auto_avg_0 / denum / avg_fsd_auto
-c_auto_avg_2 = c_auto_avg_2 / denum / avg_c_auto
-y_avg_2 = y_avg_2 / denum / avg_y
-c_auto_avg_1 = c_auto_avg_1 / denum / avg_c_auto
-y_avg_0 = y_avg_0 / denum / avg_y
-y_avg_1 = y_avg_1 / denum / avg_y
-c_auto_avg_0 = c_auto_avg_0 / denum / avg_c_auto
+c_auto_avg_u = c_auto_avg_u / denum / avg_c_auto
+c_auto_avg_v = c_auto_avg_v / denum / avg_c_auto
+c_auto_avg_w = c_auto_avg_w / denum / avg_c_auto
+y_avg_u = y_avg_u / denum / avg_y
+y_avg_v = y_avg_v / denum / avg_y
+y_avg_w = y_avg_w / denum / avg_y
+y_avg_c_auto = y_avg_c_auto / denum / avg_y
+fsd_auto_avg_u = fsd_auto_avg_u / denum / avg_fsd_auto
 
 
 end subroutine m_terms_average_pass1
 
 
+subroutine m_terms_calculate_pass2
+
+integer :: i, j, k
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+c_auto(i,j,k) = ( 1.0 - y(i,j,k) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+tflux_(i,j,k) = ( (u(i,j,k) - avg_u(i)) * (c_auto(i,j,k) - avg_c_auto(i)) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+avg_tflux(i) = avg_tflux(i) + tflux_(i,j,k) 
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+tflux_y(i,j,k) = ( (u(i,j,k) - y_avg_u(i)) * (c_auto(i,j,k) - y_avg_c_auto(i)) )
+end do
+end do
+end do
+
+
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+y_avg_tflux(i) = y_avg_tflux(i) + tflux_y(i,j,k) * y(i,j,k)
+end do
+end do
+end do
+
+
+
+end subroutine m_terms_calculate_pass2
+
+
+subroutine m_terms_average_pass2 (nfiles)
+
+integer :: nfiles
+real*8 :: denum
+
+denum=real(nfiles*nyp*nzp)
+
+avg_tflux = avg_tflux / denum 
+y_avg_tflux = y_avg_tflux / denum / avg_y
+
+
+end subroutine m_terms_average_pass2
+
+
 end module m_terms
 
-!  {'c_auto': '( 1.0 - y(i,j,k) )', 'fsd_auto': '( sqrt ( ( ( ((ddx_c_auto(i,j,k))*(ddx_c_auto(i,j,k))) + ((ddy_c_auto(i,j,k))*(ddy_c_auto(i,j,k))) ) + ((ddz_c_auto(i,j,k))*(ddz_c_auto(i,j,k))) ) ) )', 'wrate_auto': '( rxn_rate ( c_auto(i,j,k) ) )'}
-!  {'ddz_c_auto': ('ddz', 'c_auto'), 'ddy_c_auto': ('ddy', 'c_auto'), 'ddx_c_auto': ('ddx', 'c_auto')}
-!  {'fsd_auto': ['ddx_c_auto', 'ddy_c_auto', 'ddz_c_auto'], 'avg_fsd_auto': ['fsd_auto'], u'c_auto_avg_2': ['c_auto'], 'avg_c_auto': ['c_auto'], 'avg_1': [], 'ddz_c_auto': ['c_auto'], 'c_auto': [], u'y_avg_2': [], u'c_auto_avg_1': ['c_auto'], u'y_avg_0': [], u'y_avg_1': [], 'avg_0': [], 'ddy_c_auto': ['c_auto'], 'avg_2': [], u'c_auto_avg_0': ['c_auto'], 'avg_y': [], u'fsd_auto_avg_0': ['fsd_auto'], 'ddx_c_auto': ['c_auto'], 'wrate_auto': ['c_auto']}
-!  ['avg_1', 'c_auto', 'avg_c_auto', 'ddz_c_auto', u'c_auto_avg_2', u'c_auto_avg_1', u'c_auto_avg_0', u'y_avg_2', u'y_avg_0', u'y_avg_1', 'avg_0', 'ddy_c_auto', 'avg_2', 'avg_y', 'ddx_c_auto', 'fsd_auto', u'fsd_auto_avg_0', 'avg_fsd_auto', 'wrate_auto']
diff --git a/code/post.f90 b/code/post.f90
index b5e2e3c..79c6de5 100644
--- a/code/post.f90
+++ b/code/post.f90
@@ -106,11 +106,15 @@
          CALL CAL_FLUCTUATION
          CALL SAVE_SUM_FLUCTUATION
 
+      CALL m_terms_calculate_pass2
+
        ENDIF
 
       ENDDO secondloop
 
       CALL FLUCTUATION_AVG
+
+      CALL m_terms_average_pass2(countnum)
       
   !    CALL FINAL_AVG  ! edge-cold-bc-3
       CALL SAVE_AVG_RESULTS