diff --git a/code/code_gen/code_gen.py b/code/code_gen/code_gen.py
index f79e7b8..687e063 100644
--- a/code/code_gen/code_gen.py
+++ b/code/code_gen/code_gen.py
@@ -43,12 +43,12 @@ calc_grammar = """
              | "sqrt" -> sqrt
              | "rxn_rate" -> rxn_rate
 
-    ?derivative: "dx"  -> dx
-               | "d2x" -> d2x
-               | "dy" -> dy
-               | "d2y" -> d2y
-               | "dz" -> dz
-               | "d2z" -> d2z
+    ?derivative: "ddx"  -> ddx
+               | "dd2x" -> dd2x
+               | "ddy" -> ddy
+               | "dd2y" -> dd2y
+               | "ddz" -> ddz
+               | "dd2z" -> dd2z
 
     %import common.CNAME -> NAME
     %import common.NUMBER
@@ -71,6 +71,7 @@ end do
 end do
 """
 
+real_array_diff  = "call {0[0]} ( {0[0]}_{0[1]}, {0[1]} )"
 
 @v_args(inline=True)    # Affects the signatures of the methods
 class CalculateTree(Transformer):
@@ -195,8 +196,16 @@ class CalculateTree(Transformer):
 
         f_code = ""
 
+        code_dict = {}
+
         for tup in self.derived.iteritems():
-            f_code = f_code + real_array_loop.format(tup) + "\n"
+            code_dict[tup[0]] = real_array_loop.format(tup) + "\n"
+
+        for tup in self.derivatives.iteritems():
+            code_dict[tup[0]] = real_array_diff.format(tup[1]) + "\n"
+
+        for var in self.sort_vars():
+            f_code = f_code + code_dict[var]
 
         return f_code
 
@@ -211,6 +220,19 @@ class CalculateTree(Transformer):
 
         return md
 
+    def sort_vars (self):
+        order = []
+        remain = set(self.derived.iterkeys()) | set(self.derivatives.iterkeys())
+
+        while len(remain) > 0:
+            for v in remain:
+                if len(set(self.dependency[v]) & remain) == 0:
+                    order.append(v)
+                    remain.remove(v)
+                    break
+
+        return order
+
 
 tf=CalculateTree()
 
@@ -245,6 +267,8 @@ def test():
     print "! ", tf.derivatives
     print "! ", tf.dependency
 
+    print "! ", tf.sort_vars()
+
 if __name__ == '__main__':
     test()
     # main()
diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index de2b056..2bb19a6 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -2,4 +2,4 @@
 
 c_auto = 1.0 - y
 wrate_auto = rxn_rate(c_auto)
-fsd_auto = sqrt (sqr(dx(c_auto)) + sqr(dy(c_auto)) + sqr(dz(c_auto)))
+fsd_auto = sqrt (sqr(ddx(c_auto)) + sqr(ddy(c_auto)) + sqr(ddz(c_auto)))
diff --git a/code/m_calculate.f90 b/code/m_calculate.f90
index ba3cf17..ec5ceda 100644
--- a/code/m_calculate.f90
+++ b/code/m_calculate.f90
@@ -24,7 +24,7 @@ contains
 
      integer :: ierr
 
-     allocate(work(nxp, nyp, nzp, 3), stat=ierr)
+     ! allocate(work(nxp, nyp, nzp, 3), stat=ierr)
 
      allocate(xsrc(nyp, nxp), stat=ierr)
      allocate(xdst(nyp, nxp), stat=ierr)
@@ -37,7 +37,7 @@ contains
 
      subroutine m_calculate_finalize
 
-     deallocate(work)
+     ! deallocate(work)
 
      deallocate(xsrc)
      deallocate(xdst)
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 85ca4df..3f83f8f 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -9,9 +9,9 @@ implicit none
 real*8, allocatable, dimension(:,:,:) :: c_auto
 real*8, allocatable, dimension(:,:,:) :: fsd_auto
 real*8, allocatable, dimension(:,:,:) :: wrate_auto
-real*8, allocatable, dimension(:,:,:) :: dx_c_auto
-real*8, allocatable, dimension(:,:,:) :: dy_c_auto
-real*8, allocatable, dimension(:,:,:) :: dz_c_auto
+real*8, allocatable, dimension(:,:,:) :: ddz_c_auto
+real*8, allocatable, dimension(:,:,:) :: ddy_c_auto
+real*8, allocatable, dimension(:,:,:) :: ddx_c_auto
 
 
 contains
@@ -23,9 +23,9 @@ integer :: ierr
 allocate(c_auto(nxp,nyp,nzp), stat=ierr) ; c_auto = 0.
 allocate(fsd_auto(nxp,nyp,nzp), stat=ierr) ; fsd_auto = 0.
 allocate(wrate_auto(nxp,nyp,nzp), stat=ierr) ; wrate_auto = 0.
-allocate(dx_c_auto(nxp,nyp,nzp), stat=ierr) ; dx_c_auto = 0.
-allocate(dy_c_auto(nxp,nyp,nzp), stat=ierr) ; dy_c_auto = 0.
-allocate(dz_c_auto(nxp,nyp,nzp), stat=ierr) ; dz_c_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.
 
 
 end subroutine m_terms_init
@@ -36,9 +36,9 @@ subroutine m_terms_finalize
 deallocate(c_auto)
 deallocate(fsd_auto)
 deallocate(wrate_auto)
-deallocate(dx_c_auto)
-deallocate(dy_c_auto)
-deallocate(dz_c_auto)
+deallocate(ddz_c_auto)
+deallocate(ddy_c_auto)
+deallocate(ddx_c_auto)
 
 
 end subroutine m_terms_finalize
@@ -57,11 +57,14 @@ end do
 end do
 end do
 
+call ddy ( ddy_c_auto, c_auto )
+call ddz ( ddz_c_auto, c_auto )
+call ddx ( ddx_c_auto, c_auto )
 
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-fsd_auto(i,j,k) = ( sqrt ( ( ( ((dx_c_auto(i,j,k))*(dx_c_auto(i,j,k))) + ((dy_c_auto(i,j,k))*(dy_c_auto(i,j,k))) ) + ((dz_c_auto(i,j,k))*(dz_c_auto(i,j,k))) ) ) )
+fsd_auto(i,j,k) = ( 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))) ) ) )
 end do
 end do
 end do
@@ -82,6 +85,7 @@ end subroutine m_terms_calculate_instant
 
 end module m_terms
 
-!  {'c_auto': '( 1.0 - y(i,j,k) )', 'fsd_auto': '( sqrt ( ( ( ((dx_c_auto(i,j,k))*(dx_c_auto(i,j,k))) + ((dy_c_auto(i,j,k))*(dy_c_auto(i,j,k))) ) + ((dz_c_auto(i,j,k))*(dz_c_auto(i,j,k))) ) ) )', 'wrate_auto': '( rxn_rate ( c_auto(i,j,k) ) )'}
-!  {'dx_c_auto': ('dx', 'c_auto'), 'dy_c_auto': ('dy', 'c_auto'), 'dz_c_auto': ('dz', 'c_auto')}
-!  {'dx_c_auto': ['c_auto'], 'fsd_auto': ['dx_c_auto', 'dy_c_auto', 'dz_c_auto'], 'dy_c_auto': ['c_auto'], 'c_auto': [], 'wrate_auto': ['c_auto'], 'dz_c_auto': ['c_auto']}
+!  {'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'], 'ddz_c_auto': ['c_auto'], 'c_auto': [], 'ddy_c_auto': ['c_auto'], 'ddx_c_auto': ['c_auto'], 'wrate_auto': ['c_auto']}
+!  ['c_auto', 'ddy_c_auto', 'ddz_c_auto', 'ddx_c_auto', 'fsd_auto', 'wrate_auto']
diff --git a/code/post.f90 b/code/post.f90
index ab40afb..4bdf649 100644
--- a/code/post.f90
+++ b/code/post.f90
@@ -1474,6 +1474,7 @@
       INTEGER :: ierr
 
        CALL m_arrays_init
+       CALL m_calculate_init
        CALL m_terms_init
 
        ALLOCATE(old_scalar(2,nxp,nyp,nzp),STAT=ierr)  ; old_scalar=0.          ! Main variables
@@ -1564,6 +1565,7 @@
      SUBROUTINE DEALLOCATES_CLOSE
 
       CALL m_arrays_finalize
+      CALL m_calculate_finalize
       CALL m_terms_finalize
 
       DEALLOCATE(old_scalar)