diff --git a/code/code_gen/code_gen.py b/code/code_gen/code_gen.py
index 0cc9864..680741a 100644
--- a/code/code_gen/code_gen.py
+++ b/code/code_gen/code_gen.py
@@ -170,7 +170,8 @@ def test(terms_raw, report=False, latex=False):
             # print r"\end{equation}"
         print("}")
     else:
-        print(mod_form.format( cdict ))
+        from jinja2 import Template
+        print(Template(mod_form).render(**cdict))
 
 
 
diff --git a/code/code_gen/post.py b/code/code_gen/post.py
index 7488960..ec4fd96 100644
--- a/code/code_gen/post.py
+++ b/code/code_gen/post.py
@@ -1,5 +1,6 @@
 from lark import Lark, Visitor, Transformer, v_args, Token
 import warnings
+from jinja2 import Template
 
 
 class CollectDefinitions(Visitor):
@@ -355,90 +356,90 @@ integer(kind=MPI_OFFSET_KIND)  :: offset
     # Subarray version
     fmt_decl_subarray='''
 ! - file_handles and mpi_infos
-integer(kind=MPI_INTEGER_KIND) :: {field_name}_fh
-integer(kind=MPI_INTEGER_KIND) :: {field_name}_info
-integer(kind=MPI_INTEGER_KIND) :: {field_name}_filetype
+integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_fh
+integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_info
+integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_filetype
 '''
 
     fmt_init_subarray='''
-! init subarray datatype for {field_name}
+! init subarray datatype for {{ field_name }}
 block
     integer(4) :: sizes(3), subsizes(3), starts(3)
-    call MPI_INFO_CREATE({field_name}_info, mpi_err)
-    call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{field_name}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{field_name}_info,{field_name}_fh,mpi_err)
+    call MPI_INFO_CREATE({{ field_name }}_info, mpi_err)
+    call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{{ field_name }}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{{ field_name }}_info,{{ field_name }}_fh,mpi_err)
     sizes = (/ nxp, nyp, nzp /)
-    subsizes = (/ {len_xpts}, {ye} - {ys} + 1, {ze} - {zs} + 1 /)
-    starts = (/ {xs} - 1, {ys} - 1, {zs} - 1 /)
-    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts, MPI_ORDER_FORTRAN, MPI_REAL8, {field_name}_filetype, mpi_err)
-    call MPI_TYPE_COMMIT({field_name}_filetype, mpi_err)
+    subsizes = (/ {{ len_xpts }}, {{ ye }} - {{ ys }} + 1, {{ ze }} - {{ zs }} + 1 /)
+    starts = (/ {{ xs }} - 1, {{ ys }} - 1, {{ zs }} - 1 /)
+    call MPI_TYPE_CREATE_SUBARRAY(3, sizes, subsizes, starts, MPI_ORDER_FORTRAN, MPI_REAL8, {{ field_name }}_filetype, mpi_err)
+    call MPI_TYPE_COMMIT({{ field_name }}_filetype, mpi_err)
 end block
 '''
 
     fmt_final_subarray='''
 ! finalize
-call MPI_FILE_CLOSE({field_name}_fh, mpi_err)
-call MPI_INFO_FREE({field_name}_info, mpi_err)
-call MPI_TYPE_FREE({field_name}_filetype, mpi_err)
+call MPI_FILE_CLOSE({{ field_name }}_fh, mpi_err)
+call MPI_INFO_FREE({{ field_name }}_info, mpi_err)
+call MPI_TYPE_FREE({{ field_name }}_filetype, mpi_err)
 '''
 
     fmt_calc_subarray='''
 ! write to file via MPI Subarray
-count = ({len_xpts}) * ({ye} - {ys} + 1) * ({ze} - {zs} + 1)
+count = ({{ len_xpts }}) * ({{ ye }} - {{ ys }} + 1) * ({{ ze }} - {{ zs }} + 1)
 offset = export_offset(fidx) * count * 8
-call MPI_FILE_WRITE_AT({field_name}_fh, offset, {work_array}, 1, {field_name}_filetype, mpi_status, mpi_err)
+call MPI_FILE_WRITE_AT({{ field_name }}_fh, offset, {{ work_array }}, 1, {{ field_name }}_filetype, mpi_status, mpi_err)
 '''
 
     # Legacy copy version (fallback)
     fmt_decl_legacy='''
 ! - file_handles and mpi_infos
-integer(kind=MPI_INTEGER_KIND) :: {field_name}_fh
-integer(kind=MPI_INTEGER_KIND) :: {field_name}_info
+integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_fh
+integer(kind=MPI_INTEGER_KIND) :: {{ field_name }}_info
 
 ! - buffer
-real(real64), allocatable, dimension(:,:,:) :: {field_name}_export_array
-integer, allocatable, dimension(:) :: {field_name}_xpts
+real(real64), allocatable, dimension(:,:,:) :: {{ field_name }}_export_array
+integer, allocatable, dimension(:) :: {{ field_name }}_xpts
 '''
 
     fmt_init_legacy='''
 ! init
-call MPI_INFO_CREATE({field_name}_info, mpi_err)
-call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{field_name}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{field_name}_info,{field_name}_fh,mpi_err)
-allocate({field_name}_export_array(1:{len_xpts},{ys}:{ye},{zs}:{ze}), stat=ierr)
+call MPI_INFO_CREATE({{ field_name }}_info, mpi_err)
+call MPI_FILE_OPEN(MPI_COMM_TASK,'export-{{ field_name }}.dat',MPI_MODE_WRONLY+MPI_MODE_CREATE,{{ field_name }}_info,{{ field_name }}_fh,mpi_err)
+allocate({{ field_name }}_export_array(1:{{ len_xpts }},{{ ys }}:{{ ye }},{{ zs }}:{{ ze }}), stat=ierr)
 if (ierr /= 0) then
-    write(0,*) 'Error: allocation of {field_name}_export_array failed on process', myid
+    write(0,*) 'Error: allocation of {{ field_name }}_export_array failed on process', myid
     call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
 end if
-{field_name}_export_array = 0.
-allocate({field_name}_xpts(1:{len_xpts}), stat=ierr)
+{{ field_name }}_export_array = 0.
+allocate({{ field_name }}_xpts(1:{{ len_xpts }}), stat=ierr)
 if (ierr /= 0) then
-    write(0,*) 'Error: allocation of {field_name}_xpts failed on process', myid
+    write(0,*) 'Error: allocation of {{ field_name }}_xpts failed on process', myid
     call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err)
 end if
-{xpts_init}
+{{ xpts_init }}
 '''
 
     fmt_final_legacy='''
 ! finalize
-call MPI_FILE_CLOSE({field_name}_fh, mpi_err)
-call MPI_INFO_FREE({field_name}_info, mpi_err)
-deallocate({field_name}_export_array)
-deallocate({field_name}_xpts)
+call MPI_FILE_CLOSE({{ field_name }}_fh, mpi_err)
+call MPI_INFO_FREE({{ field_name }}_info, mpi_err)
+deallocate({{ field_name }}_export_array)
+deallocate({{ field_name }}_xpts)
 '''
 
     fmt_calc_legacy='''
 ! copy to array for export
-do k = {zs}, {ze}
-do j = {ys}, {ye}
-do i = 1, {len_xpts}
-{field_name}_export_array(i,j,k) = {work_array}({field_name}_xpts(i),j,k)
+do k = {{ zs }}, {{ ze }}
+do j = {{ ys }}, {{ ye }}
+do i = 1, {{ len_xpts }}
+{{ field_name }}_export_array(i,j,k) = {{ work_array }}({{ field_name }}_xpts(i),j,k)
 end do
 end do
 end do
 
 ! write to file
-count = ({len_xpts}) * ({ye} - {ys} + 1) * ({ze} - {zs} + 1)
+count = ({{ len_xpts }}) * ({{ ye }} - {{ ys }} + 1) * ({{ ze }} - {{ zs }} + 1)
 offset = export_offset(fidx) * count * 8
-call MPI_FILE_WRITE_AT({field_name}_fh, offset, {field_name}_export_array, count, MPI_REAL8, mpi_status, mpi_err)
+call MPI_FILE_WRITE_AT({{ field_name }}_fh, offset, {{ field_name }}_export_array, count, MPI_REAL8, mpi_status, mpi_err)
 '''
 
     def __init__ (self, name, attr, parent):
@@ -457,18 +458,18 @@ call MPI_FILE_WRITE_AT({field_name}_fh, offset, {field_name}_export_array, count
 
         self.params.setdefault("field_name", self.name)
 
-        self.params.setdefault("len_xpts", "({xe} - {xs} + 1)".format(**self.params))
+        self.params.setdefault("len_xpts", f"({self.params['xe']} - {self.params['xs']} + 1)")
         
         fmt_xpts_init = '''
-do i = {xs}, {xe}
-{field_name}_xpts(i-{xs}+1) = i
+do i = {{ xs }}, {{ xe }}
+{{ field_name }}_xpts(i-{{ xs }}+1) = i
 end do
 '''
-        self.params.setdefault("xpts_init", fmt_xpts_init.format(**self.params))
+        self.params.setdefault("xpts_init", Template(fmt_xpts_init).render(**self.params))
 
         try:
         # Sampling at listed x coordinates
-            fmt_xpts_init_list = "{field_name}_xpts = (/ {list_xpts} /)"
+            fmt_xpts_init_list = "{{ field_name }}_xpts = (/ {{ list_xpts }} /)"
 
             import sys
             raw_xpts = self.params["xpts"]
@@ -476,7 +477,7 @@ end do
             len_xpts = len(int_xpts)
             self.params["len_xpts"] = len_xpts
             self.params["list_xpts"] = ",".join(map(str, int_xpts))
-            self.params["xpts_init"] = fmt_xpts_init_list.format(**self.params)
+            self.params["xpts_init"] = Template(fmt_xpts_init_list).render(**self.params)
         except KeyError:
             pass
 
@@ -486,27 +487,27 @@ end do
     def code (self):
         self.params["work_array"] = self.parent.array
         if self.use_subarray:
-            return FieldExporter.fmt_calc_subarray.format(**self.params)
+            return Template(FieldExporter.fmt_calc_subarray).render(**self.params)
         else:
-            return FieldExporter.fmt_calc_legacy.format(**self.params)
+            return Template(FieldExporter.fmt_calc_legacy).render(**self.params)
 
     def code_decl (self):
         if self.use_subarray:
-            return FieldExporter.fmt_decl_subarray.format(**self.params)
+            return Template(FieldExporter.fmt_decl_subarray).render(**self.params)
         else:
-            return FieldExporter.fmt_decl_legacy.format(**self.params)
+            return Template(FieldExporter.fmt_decl_legacy).render(**self.params)
 
     def code_alloc (self):
         if self.use_subarray:
-            return FieldExporter.fmt_init_subarray.format(**self.params)
+            return Template(FieldExporter.fmt_init_subarray).render(**self.params)
         else:
-            return FieldExporter.fmt_init_legacy.format(**self.params)
+            return Template(FieldExporter.fmt_init_legacy).render(**self.params)
 
     def code_free (self):
         if self.use_subarray:
-            return FieldExporter.fmt_final_subarray.format(**self.params)
+            return Template(FieldExporter.fmt_final_subarray).render(**self.params)
         else:
-            return FieldExporter.fmt_final_legacy.format(**self.params)
+            return Template(FieldExporter.fmt_final_legacy).render(**self.params)
 
 
 
@@ -545,11 +546,11 @@ class Field (FieldBase):
 
     def code (self, alloc=None):
         real_array_loop  = """
-! {1}
+! {{ comment }}
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-{0[0]}(i,j,k) = {0[1]}
+{{ array }}(i,j,k) = {{ rhs }}
 end do
 end do
 end do
@@ -557,7 +558,7 @@ end do
         self.array = alloc[self.name] if alloc else self.name
 
         rhs = ExpToCode(self.fdict).transform(self.exp)
-        calculation_code = real_array_loop.format((self.array, rhs), self.comment)
+        calculation_code = Template(real_array_loop).render(comment=self.comment, array=self.array, rhs=rhs)
 
         export_code = ( self.exporter.code() if self.export_on() else "")
 
@@ -589,11 +590,11 @@ class FluctuationField (FieldBase):
 
     def code (self, alloc=None):
         real_array_loop  = """
-! {1}
+! {{ comment }}
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-{0[0]}(i,j,k) = {0[1]}
+{{ array }}(i,j,k) = {{ rhs }}
 end do
 end do
 end do
@@ -605,7 +606,7 @@ end do
         if self.field.is_fluctuation():
             rhs = rhs.format(self.w)
 
-        return real_array_loop.format((self.array, rhs), self.comment)
+        return Template(real_array_loop).render(comment=self.comment, array=self.array, rhs=rhs)
 
     @classmethod
     def id (cls, w, field):
@@ -702,7 +703,7 @@ class AveragedField (FieldBase):
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-{0}(i) = {0}(i) + {1}
+{{ name }}(i) = {{ name }}(i) + {{ arrname }}
 end do
 end do
 end do
@@ -711,20 +712,19 @@ end do
         if self.weighted is not None:
             arrname = arrname + " * " + self.w.array + "(i,j,k)"
 
-        return avg_array_sum.format(self.name, arrname)
+        return Template(avg_array_sum).render(name=self.name, arrname=arrname)
 
 
     def code_avg (self):
 
         avg_array_divide  = """
-call MPI_ALLREDUCE(MPI_IN_PLACE, {0}, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
+call MPI_ALLREDUCE(MPI_IN_PLACE, {{ name }}, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
 
-{0} = {0} {1} / denum
+{{ name }} = {{ name }} {{ dWeight }} / denum
 """
-        meanw = "/ avg_{}"
-        dWeight = (meanw.format(self.weighted) if self.weighted else "")
+        dWeight = (f"/ avg_{self.weighted}" if self.weighted else "")
 
-        return avg_array_divide.format(self.name, dWeight)
+        return Template(avg_array_divide).render(name=self.name, dWeight=dWeight)
 
     def isWeighted (self):
         return self.weighted is not None
@@ -1027,29 +1027,34 @@ integer :: i
 open (200, file="qEdge_X.dat")
 write (200,*) output_header
 do i=1,nxp
- write (200,'({0}e20.10)') real(i)*hxp, {1}
+ write (200,'({{ num_args }}e20.10)') real(i)*hxp, {{ formatted_avglist }}
 end do
 close (200)
 
 open (200, file="d1.dat")
-{2}
+{{ deriv1_lines }}
 close (200)
 
 open (200, file="d2.dat")
-{3}
+{{ deriv2_lines }}
 close (200)
 '''
 
         avgarr = "{}(i)"
         deriv1_avgarr = """call ddx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
         deriv2_avgarr = """call d2dx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
-        write_avg = avg_array_write.format(
-                len(self.averaged)+1,
-                ", ".join(map(avgarr.format, avglist)),
-                "\n".join(map(deriv1_avgarr.format, avglist)),
-                "\n".join(map(deriv2_avgarr.format, avglist)),
-                )
 
+        num_args = len(self.averaged) + 1
+        formatted_avglist = ", ".join(map(avgarr.format, avglist))
+        deriv1_lines = "\n".join(map(deriv1_avgarr.format, avglist))
+        deriv2_lines = "\n".join(map(deriv2_avgarr.format, avglist))
+
+        write_avg = Template(avg_array_write).render(
+            num_args=num_args,
+            formatted_avglist=formatted_avglist,
+            deriv1_lines=deriv1_lines,
+            deriv2_lines=deriv2_lines
+        )
         return write_avg
 
 
diff --git a/code/code_gen/resources/m_template.py b/code/code_gen/resources/m_template.py
index 572c575..e5134d3 100644
--- a/code/code_gen/resources/m_template.py
+++ b/code/code_gen/resources/m_template.py
@@ -1,4 +1,4 @@
-mod_form='''module m_{0[module_name]}
+mod_form='''module m_{{ module_name }}
 
 use, intrinsic :: iso_fortran_env, only: real64
 use m_openmpi
@@ -8,77 +8,77 @@ use m_calculate
 
 implicit none
 
-{0[module_data]}
+{{ module_data }}
 
 contains
 
-subroutine m_{0[module_name]}_init
+subroutine m_{{ module_name }}_init
 
 integer :: ierr
 integer :: i, j, k
 
-{0[module_init]}
+{{ module_init }}
 
-end subroutine m_{0[module_name]}_init
+end subroutine m_{{ module_name }}_init
 
 
-subroutine m_{0[module_name]}_finalize
+subroutine m_{{ module_name }}_finalize
 
-{0[module_finalize]}
+{{ module_finalize }}
 
-end subroutine m_{0[module_name]}_finalize
+end subroutine m_{{ module_name }}_finalize
 
 
-subroutine m_{0[module_name]}_calculate_pass1(fidx)
+subroutine m_{{ module_name }}_calculate_pass1(fidx)
 
 integer :: fidx
 integer :: i, j, k
 
-{0[module_pass1]}
+{{ module_pass1 }}
 
-end subroutine m_{0[module_name]}_calculate_pass1
+end subroutine m_{{ module_name }}_calculate_pass1
 
 
-subroutine m_{0[module_name]}_average_pass1 (nfiles)
+subroutine m_{{ module_name }}_average_pass1 (nfiles)
 
 integer :: nfiles
 real*8 :: denum
 
 denum=real(nfiles*nyp*nzp)
 
-{0[module_pass1_avg]}
+{{ module_pass1_avg }}
 
-end subroutine m_{0[module_name]}_average_pass1
+end subroutine m_{{ module_name }}_average_pass1
 
 
-subroutine m_{0[module_name]}_calculate_pass2(fidx)
+subroutine m_{{ module_name }}_calculate_pass2(fidx)
 
 integer :: fidx
 integer :: i, j, k
 
-{0[module_pass2]}
+{{ module_pass2 }}
 
-end subroutine m_{0[module_name]}_calculate_pass2
+end subroutine m_{{ module_name }}_calculate_pass2
 
 
-subroutine m_{0[module_name]}_average_pass2 (nfiles)
+subroutine m_{{ module_name }}_average_pass2 (nfiles)
 
 integer :: nfiles
 real*8 :: denum
 
 denum=real(nfiles*nyp*nzp)
 
-{0[module_pass2_avg]}
+{{ module_pass2_avg }}
 
-end subroutine m_{0[module_name]}_average_pass2
+end subroutine m_{{ module_name }}_average_pass2
 
 
-subroutine m_{0[module_name]}_write_result
+subroutine m_{{ module_name }}_write_result
 
-{0[module_write_result]}
+{{ module_write_result }}
 
-end subroutine m_{0[module_name]}_write_result
+end subroutine m_{{ module_name }}_write_result
 
 
-end module m_{0[module_name]}
+end module m_{{ module_name }}
 '''
diff --git a/code/termspecs/termspec.banded_fsd-0.05-0.95-0.1 b/code/termspecs/termspec.banded_fsd-0.05-0.95-0.1
index 67b6c1a..2496e3d 100644
--- a/code/termspecs/termspec.banded_fsd-0.05-0.95-0.1
+++ b/code/termspecs/termspec.banded_fsd-0.05-0.95-0.1
@@ -18,8 +18,6 @@ mask9 = $threshold_min_max(c, 0.85, 0.95)
 
 fsd_orig = sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c)))
 
-fsd = fsd_orig * mask
-
 fsd1 = fsd_orig * mask1
 fsd2 = fsd_orig * mask2
 fsd3 = fsd_orig * mask3