diff --git a/code/code_gen/code_gen.py b/code/code_gen/code_gen.py
index c20e9d0..0b562ab 100644
--- a/code/code_gen/code_gen.py
+++ b/code/code_gen/code_gen.py
@@ -30,14 +30,16 @@ calc_grammar = """
          | "-" atom         -> neg
          | NAME             -> var
          | "(" sum ")"
-         | function "(" sum ")" -> fcall
+         | mathfunc "(" sum ")" -> fcall
+         | derivative "(" NAME ")" -> dnx
 
-    ?function: "log" -> log
+    ?mathfunc: "log" -> log
              | "exp" -> exp
-             | "dx"  -> dx
-             | "d2x" -> d2x
              | "sqrt" -> sqrt
 
+    ?derivative: "dx"  -> dx
+               | "d2x" -> d2x
+
     %import common.CNAME -> NAME
     %import common.NUMBER
     %import common.WS
@@ -62,13 +64,14 @@ end do
 
 @v_args(inline=True)    # Affects the signatures of the methods
 class CalculateTree(Transformer):
-    from operator import add, sub, mul, truediv as div, neg
+
     number = float
 
     def __init__(self):
         self.primary = []
         self.derived = {}
         self.averaged = {}
+        self.derivatives = {}
 
     def varlist(self, *args):
         for arg in args:
@@ -80,7 +83,14 @@ class CalculateTree(Transformer):
         return "{} = {}".format(name, value)
 
     def var(self, name):
-        return name+"(i,j,k)"
+        return name + "(i,j,k)"
+
+    def dnx (self, a, b):
+        signature = "{}_{}".format(a, b)
+
+        self.derivatives[signature] = (a, b)
+
+        return "{}_{}".format(a, b) + "(i,j,k)"
 
     def fcall (self, a, b):
         return "( {} ( {} ) )".format(a, b)
@@ -107,13 +117,22 @@ class CalculateTree(Transformer):
     sqrt = lambda self : "sqrt"
 
     def array_decl (self):
-        return "\n".join(map(real_array_decl.format, self.derived.iterkeys()))
+        return "\n".join(
+                map(real_array_decl.format, self.derived.iterkeys())
+                + map(real_array_decl.format, self.derivatives.iterkeys())
+                )
 
     def array_init (self):
-        return "\n".join(map(real_array_alloc.format, self.derived.iterkeys()))
+        return "\n".join(
+                map(real_array_alloc.format, self.derived.iterkeys())
+                + map(real_array_alloc.format, self.derivatives.iterkeys())
+                )
 
     def array_final (self):
-        return "\n".join(map(real_array_free.format, self.derived.iterkeys()))
+        return "\n".join(
+                map(real_array_free.format, self.derived.iterkeys())
+                + map(real_array_free.format, self.derivatives.iterkeys())
+                )
 
     def array_ci (self):
         return "\n".join(map(real_array_loop.format, self.derived.iteritems()))
@@ -134,6 +153,7 @@ tf=CalculateTree()
 calc_parser = Lark(calc_grammar, parser='lalr' , transformer=tf)
 calc = calc_parser.parse
 
+import sys
 
 def main():
     while True:
@@ -146,7 +166,6 @@ def main():
 
 def test():
 
-
     with open("resources/m_template.f90") as template_file:
         mod_form = template_file.read()
 
diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input
index 44e4e78..d1cc4e7 100644
--- a/code/code_gen/terms.input
+++ b/code/code_gen/terms.input
@@ -2,3 +2,4 @@
 
 c_auto = 1 - y
 wrate_auto = 21000 * (1 - c_auto) * exp (- 26.7 / (1 + 3 * c_auto))
+dcdx_auto = dx (c_auto)
diff --git a/code/m_terms.f90 b/code/m_terms.f90
index 4920bc9..a320a47 100644
--- a/code/m_terms.f90
+++ b/code/m_terms.f90
@@ -6,7 +6,9 @@ use m_arrays
 implicit none
 
 real*8, allocatable, dimension(:,:,:) :: c_auto
+real*8, allocatable, dimension(:,:,:) :: dcdx_auto
 real*8, allocatable, dimension(:,:,:) :: wrate_auto
+real*8, allocatable, dimension(:,:,:) :: dx_c_auto
 
 contains
 
@@ -15,7 +17,9 @@ subroutine m_terms_init
 integer :: ierr
 
 allocate(c_auto(nxp,nyp,nzp), stat=ierr) ; c_auto = 0.
+allocate(dcdx_auto(nxp,nyp,nzp), stat=ierr) ; dcdx_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.
 
 end subroutine m_terms_init
 
@@ -23,7 +27,9 @@ end subroutine m_terms_init
 subroutine m_terms_finalize
 
 deallocate(c_auto)
+deallocate(dcdx_auto)
 deallocate(wrate_auto)
+deallocate(dx_c_auto)
 
 end subroutine m_terms_finalize
 
@@ -42,6 +48,15 @@ end do
 end do
 
 
+do k = 1, nzp
+do j = 1, nyp
+do i = 1, nxp
+dcdx_auto(i,j,k) = dx_c_auto(i,j,k)
+end do
+end do
+end do
+
+
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp