Update of macros for 4th order differencing.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinBase/ADMMacros/trunk@53 b1d164ef-f17a-46e7-89d4-021c7118ef4e

23 files changed, 659 insertions, 564 deletions

diff --git a/src/macro/ADM_Derivative.h b/src/macro/ADM_Derivative.h
new file mode 100644
index 0000000..af8cc91
--- /dev/null
+++ b/src/macro/ADM_Derivative.h
@@ -0,0 +1,93 @@
+/*@@
+  @header   ADM_Derivative.h
+  @date     June 2002
+  @author   Denis Pollney
+  @desc
+            Derivative operators.
+  @enddesc
+@@*/
+
+#ifndef ADM_DERIVATIVE_H
+#define ADM_DERIVATIVE_H
+
+#include "ADM_Spacing.h"
+
+/*
+ *   2nd order operators.
+ */ 
+
+#define ADM_DX_2(var,i,j,k)   I2DX*(var(i+1,j,k) - var(i-1,j,k))
+#define ADM_DY_2(var,i,j,k)   I2DY*(var(i,j+1,k) - var(i,j-1,k))
+#define ADM_DZ_2(var,i,j,k)   I2DZ*(var(i,j,k+1) - var(i,j,k-1))
+
+#define ADM_DXX_2(var,i,j,k)  IDXX*(var(i+1,j,k) - 2.d0*var(i,j,k) \
+                                    + var(i-1,j,k))
+#define ADM_DYY_2(var,i,j,k)  IDYY*(var(i,j+1,k) - 2.d0*var(i,j,k) \
+                                    + var(i,j-1,k))
+#define ADM_DZZ_2(var,i,j,k)  IDZZ*(var(i,j,k+1) - 2.d0*var(i,j,k) \
+                                    + var(i,j,k-1))
+
+#define ADM_DXY_2(var,i,j,k)  IDXY*(var(i+1,j+1,k) - var(i+1,j-1,k) \
+                                - var(i-1,j+1,k) + var(i-1,j-1,k))
+#define ADM_DXZ_2(var,i,j,k)  IDXZ*(var(i+1,j,k+1) - var(i+1,j,k-1) \
+                                - var(i-1,j,k+1) + var(i-1,j,k-1))
+#define ADM_DYZ_2(var,i,j,k)  IDYZ*(var(i,j+1,k+1) - var(i,j+1,k-1) \
+                                - var(i,j-1,k+1) + var(i,j-1,k-1))
+/*
+ *   4th order operators.
+ */
+
+#define ADM_DX_4(var,i,j,k)   (I2DX/6.d0)*(-var(i+2,j,k) + 8.d0*var(i+1,j,k) \
+                                           - 8.d0*var(i-1,j,k) + var(i-2,j,k))
+#define ADM_DY_4(var,i,j,k)   (I2DY/6.d0)*(-var(i,j+2,k) + 8.d0*var(i,j+1,k) \
+                                           - 8.d0*var(i,j-1,k) + var(i,j-2,k))
+#define ADM_DZ_4(var,i,j,k)   (I2DZ/6.d0)*(-var(i,j,k+2) + 8.d0*var(i,j,k+1) \
+                                           - 8.d0*var(i,j,k-1) + var(i,j,k-2))
+
+#define ADM_DXX_4(var,i,j,k)  (IDXX/12.d0)*(-var(i+2,j,k) + 16.d0*var(i+1,j,k)\
+                        - 30.d0*var(i,j,k) + 16.d0*var(i-1,j,k) - var(i-2,j,k))
+#define ADM_DYY_4(var,i,j,k)  (IDYY/12.d0)*(-var(i,j+2,k) + 16.d0*var(i,j+1,k)\
+                        - 30.d0*var(i,j,k) + 16.d0*var(i,j-1,k) - var(i,j-2,k))
+#define ADM_DZZ_4(var,i,j,k)  (IDZZ/12.d0)*(-var(i,j,k+2) + 16.d0*var(i,j,k+1)\
+                        - 30.d0*var(i,j,k) + 16.d0*var(i,j,k-1) - var(i,j,k-2))
+
+#define ADM_DXY_4(var,i,j,k)  (IDXY/36.d0)* \
+  (var(i+2,j+2,k) - var(i+2,j-2,k) - var(i-2,j+2,k) + var(i-2,j-2,k) \
+ + 8.d0*(-var(i+2,j+1,k) + var(i+2,j-1,k) - var(i+1,j+2,k) + var(i+1,j-2,k) \
+         + var(i-2,j+1,k) - var(i-2,j-1,k) - var(i-1,j-2,k) + var(i-1,j+2,k)) \
+ + 64.d0*(var(i+1,j+1,k) - var(i+1,j-1,k) - var(i-1,j+1,k) + var(i-1,j-1,k)))
+#define ADM_DXZ_4(var,i,j,k)  (IDXZ/36.d0)* \
+  (var(i+2,j,k+2) - var(i+2,j,k-2) - var(i-2,j,k+2) + var(i-2,j,k-2) \
+ + 8.d0*(-var(i+2,j,k+1) + var(i+2,j,k-1) - var(i+1,j,k+2) + var(i+1,j,k-2) \
+         + var(i-2,j,k+1) - var(i-2,j,k-1) - var(i-1,j,k-2) + var(i-1,j,k+2)) \
+ + 64.d0*(var(i+1,j,k+1) - var(i+1,j,k-1) - var(i-1,j,k+1) + var(i-1,j,k-1)))
+#define ADM_DYZ_4(var,i,j,k)  (IDYZ/36.d0)* \
+  (var(i,j+2,k+2) - var(i,j+2,k-2) - var(i,j-2,k+2) + var(i,j-2,k-2) \
+ + 8.d0*(-var(i,j+2,k+1) + var(i,j+2,k-1) - var(i,j+1,k+2) + var(i,j+1,k-2) \
+         + var(i,j-2,k+1) - var(i,j-2,k-1) - var(i,j-1,k-2) + var(i,j-1,k+2)) \
+ + 64.d0*(var(i,j+1,k+1) - var(i,j+1,k-1) - var(i,j-1,k+1) + var(i,j-1,k-1)))
+
+
+#define ADM_ADV_DX_P1(var,i,j,k) betax(i,j,k)*IDX*(var(i+i,j,k) - var(i,j,k))
+#define ADM_ADV_DY_P1(var,i,j,k) betay(i,j,k)*IDY*(var(i,j+1,k) - var(i,j,k))
+#define ADM_ADV_DZ_P1(var,i,j,k) betaz(i,j,k)*IDZ*(var(i,j,k+1) - var(i,j,k))
+
+#define ADM_ADV_DX_M1(var,i,j,k) betax(i,j,k)*IDX*(var(i,j,k) - var(i-1,j,k))
+#define ADM_ADV_DY_M1(var,i,j,k) betay(i,j,k)*IDY*(var(i,j,k) - var(i,j-1,k))
+#define ADM_ADV_DZ_M1(var,i,j,k) betaz(i,j,k)*IDZ*(var(i,j,k) - var(i,j,k-1))
+
+#define ADM_ADV_DX_P2(var,i,j,k) betax(i,j,k)*I2DX*(-3.d0*var(i,j,k) \
+					    + 4.d0*var(i+1,j,k) - var(i+2,j,k))
+#define ADM_ADV_DY_P2(var,i,j,k) betay(i,j,k)*I2DY*(-3.d0*var(i,j,k) \
+					    + 4.d0*var(i,j+1,k) - var(i,j+2,k))
+#define ADM_ADV_DZ_P2(var,i,j,k) betaz(i,j,k)*I2DZ*(-3.d0*var(i,j,k) \
+					    + 4.d0*var(i,j,k+1) - var(i,j,k+2))
+
+#define ADM_ADV_DX_M2(var,i,j,k) betax(i,j,k)*I2DX*(3.d0*var(i,j,k) \
+					   - 4.d0*var(i-1,j,k) + var(i-2,j,k))
+#define ADM_ADV_DY_M2(var,i,j,k) betay(i,j,k)*I2DY*(3.d0*var(i,j,k) \
+					   - 4.d0*var(i,j-1,k) + var(i,j-2,k))
+#define ADM_ADV_DZ_M2(var,i,j,k) betaz(i,j,k)*I2DZ*(3.d0*var(i,j,k) \
+					   - 4.d0*var(i,j,k-1) + var(i,j,k-2))
+
+#endif
diff --git a/src/macro/ADM_Spacing.h b/src/macro/ADM_Spacing.h
new file mode 100644
index 0000000..234d554
--- /dev/null
+++ b/src/macro/ADM_Spacing.h
@@ -0,0 +1,34 @@
+/*@@
+  @header   ADM_Spacing.h
+  @date     June 2002
+  @author   Denis Pollney
+  @desc
+            Calculate various spacing dependent scalars.
+  @enddesc
+@@*/
+
+#ifndef ADM_SPACING_H
+#define ADM_SPACING_H
+
+       dt=cctk_delta_time
+       dx=cctk_delta_space(1)
+       dy=cctk_delta_space(2)
+       dz=cctk_delta_space(3)
+
+       idx=1.d0/dx
+       idy=1.d0/dy
+       idz=1.d0/dz
+
+       i2dx=0.5d0*idx
+       i2dy=0.5d0*idy
+       i2dz=0.5d0*idz
+
+       idxx=idx**2
+       idyy=idy**2
+       idzz=idz**2
+
+       idxy=i2dx*i2dy
+       idxz=i2dx*i2dz
+       idyz=i2dy*i2dz
+
+#endif
diff --git a/src/macro/ADM_Spacing_declare.h b/src/macro/ADM_Spacing_declare.h
new file mode 100644
index 0000000..21a1ee6
--- /dev/null
+++ b/src/macro/ADM_Spacing_declare.h
@@ -0,0 +1,13 @@
+/*@@
+  @header   ADM_Spacing.h
+  @date     June 2002
+  @author   Denis Pollney
+  @desc
+            Calculate various spacing dependent scalars.
+  @enddesc
+@@*/
+
+       CCTK_REAL :: dt, dx, dy, dz
+       CCTK_REAL :: idx, idy, idz, i2dx, i2dy, i2dz
+       CCTK_REAL :: idxx, idxy, idxz, idyy, idyz, idzz
+
diff --git a/src/macro/DA_guts.h b/src/macro/DA_guts.h
index 2b3ea63..cceb5c8 100644
--- a/src/macro/DA_guts.h
+++ b/src/macro/DA_guts.h
@@ -12,13 +12,17 @@
 
 #ifdef FCODE 
 
-      DA_OO2DX = 1D0/(2D0*DA_DX)
-      DA_OO2DY = 1D0/(2D0*DA_DY)
-      DA_OO2DZ = 1D0/(2D0*DA_DZ)
-
-      DA_DXDA = DA_OO2DX*(DA_A_IP - DA_A_IM)
-      DA_DYDA = DA_OO2DY*(DA_A_JP - DA_A_JM)
-      DA_DZDA = DA_OO2DZ*(DA_A_KP - DA_A_KM)
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DA_DXDA = ADM_DX_2(alp,i,j,k)
+        DA_DYDA = ADM_DY_2(alp,i,j,k)
+        DA_DZDA = ADM_DZ_2(alp,i,j,k)
+      else
+        DA_DXDA = ADM_DX_4(alp,i,j,k)
+        DA_DYDA = ADM_DY_4(alp,i,j,k)
+        DA_DZDA = ADM_DZ_4(alp,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DDA_guts.h b/src/macro/DDA_guts.h
index de4c849..7f115f0 100644
--- a/src/macro/DDA_guts.h
+++ b/src/macro/DDA_guts.h
@@ -12,20 +12,23 @@
 
 #ifdef FCODE 
 
-      DDA_OODX2 = 1D0/(DDA_DX*DDA_DX)
-      DDA_OODY2 = 1D0/(DDA_DY*DDA_DY)
-      DDA_OODZ2 = 1D0/(DDA_DZ*DDA_DZ)
-      DDA_OO4DXDY = 1D0/(4D0*DDA_DX*DDA_DY)
-      DDA_OO4DXDZ = 1D0/(4D0*DDA_DX*DDA_DZ)
-      DDA_OO4DYDZ = 1D0/(4D0*DDA_DY*DDA_DZ)
-
-      DDA_DXXDA = DDA_OODX2*(DDA_A_IP  - 2D0*DDA_A + DDA_A_IM)
-      DDA_DYYDA = DDA_OODY2*(DDA_A_JP  - 2D0*DDA_A + DDA_A_JM)
-      DDA_DZZDA = DDA_OODZ2*(DDA_A_KP  - 2D0*DDA_A + DDA_A_KM)
-
-      DDA_DXYDA = DDA_OO4DXDY*(DDA_A_IPJP-DDA_A_IPJM-DDA_A_IMJP+DDA_A_IMJM)
-      DDA_DXZDA = DDA_OO4DXDZ*(DDA_A_IPKP-DDA_A_IPKM-DDA_A_IMKP+DDA_A_IMKM)
-      DDA_DYZDA = DDA_OO4DYDZ*(DDA_A_JPKP-DDA_A_JPKM-DDA_A_JMKP+DDA_A_JMKM)
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DDA_DXXDA = ADM_DXX_2(alp,i,j,k)
+        DDA_DXYDA = ADM_DXY_2(alp,i,j,k)
+        DDA_DXZDA = ADM_DXZ_2(alp,i,j,k)
+        DDA_DYYDA = ADM_DYY_2(alp,i,j,k)
+        DDA_DYZDA = ADM_DYZ_2(alp,i,j,k)
+        DDA_DZZDA = ADM_DZZ_2(alp,i,j,k)
+      else
+        DDA_DXXDA = ADM_DXX_4(alp,i,j,k)
+        DDA_DXYDA = ADM_DXY_4(alp,i,j,k)
+        DDA_DXZDA = ADM_DXZ_4(alp,i,j,k)
+        DDA_DYYDA = ADM_DYY_4(alp,i,j,k)
+        DDA_DYZDA = ADM_DYZ_4(alp,i,j,k)
+        DDA_DZZDA = ADM_DZZ_4(alp,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DXDB_guts.h b/src/macro/DXDB_guts.h
index c47e81d..45a773c 100644
--- a/src/macro/DXDB_guts.h
+++ b/src/macro/DXDB_guts.h
@@ -17,11 +17,17 @@
 
 #ifdef FCODE 
 
-      DXDB_OO2DX = 1.0D0/(2.0D0*DXDB_DX)
-
-      DXDB_DXDBX = DXDB_OO2DX*(DXDB_BX_IP - DXDB_BX_IM)
-      DXDB_DXDBY = DXDB_OO2DX*(DXDB_BY_IP - DXDB_BY_IM)
-      DXDB_DXDBZ = DXDB_OO2DX*(DXDB_BZ_IP - DXDB_BZ_IM)
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DXDB_DXDBX = ADM_DX_2(betax,i,j,k)
+        DXDB_DXDBY = ADM_DX_2(betay,i,j,k)
+        DXDB_DXDBZ = ADM_DX_2(betaz,i,j,k)
+      else
+        DXDB_DXDBX = ADM_DX_4(betax,i,j,k)
+        DXDB_DXDBY = ADM_DX_4(betay,i,j,k)
+        DXDB_DXDBZ = ADM_DX_4(betaz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DXDCG_guts.h b/src/macro/DXDCG_guts.h
index 464fdb7..b922068 100644
--- a/src/macro/DXDCG_guts.h
+++ b/src/macro/DXDCG_guts.h
@@ -21,15 +21,23 @@
 
 #ifdef FCODE 
 
-      DXDCG_OO2DX = 1D0/(2D0*DXDCG_DX)
-    
-      DXDCG_DXDCGXX = DXDCG_OO2DX*(DXDCG_GXX_IP - DXDCG_GXX_IM)
-      DXDCG_DXDCGXY = DXDCG_OO2DX*(DXDCG_GXY_IP - DXDCG_GXY_IM)
-      DXDCG_DXDCGXZ = DXDCG_OO2DX*(DXDCG_GXZ_IP - DXDCG_GXZ_IM)
-      DXDCG_DXDCGYY = DXDCG_OO2DX*(DXDCG_GYY_IP - DXDCG_GYY_IM)
-      DXDCG_DXDCGYZ = DXDCG_OO2DX*(DXDCG_GYZ_IP - DXDCG_GYZ_IM)
-      DXDCG_DXDCGZZ = DXDCG_OO2DX*(DXDCG_GZZ_IP - DXDCG_GZZ_IM)
-
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DXDCG_DXDCGXX = ADM_DX_2(gxx,i,j,k)
+        DXDCG_DXDCGXY = ADM_DX_2(gxy,i,j,k)
+        DXDCG_DXDCGXZ = ADM_DX_2(gxz,i,j,k)
+        DXDCG_DXDCGYY = ADM_DX_2(gyy,i,j,k)
+        DXDCG_DXDCGYZ = ADM_DX_2(gyz,i,j,k)
+        DXDCG_DXDCGZZ = ADM_DX_2(gzz,i,j,k)
+      else
+        DXDCG_DXDCGXX = ADM_DX_4(gxx,i,j,k)
+        DXDCG_DXDCGXY = ADM_DX_4(gxy,i,j,k)
+        DXDCG_DXDCGXZ = ADM_DX_4(gxz,i,j,k)
+        DXDCG_DXDCGYY = ADM_DX_4(gyy,i,j,k)
+        DXDCG_DXDCGYZ = ADM_DX_4(gyz,i,j,k)
+        DXDCG_DXDCGZZ = ADM_DX_4(gzz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/DXDK_guts.h b/src/macro/DXDK_guts.h
index 6a6df3f..fd8af92 100644
--- a/src/macro/DXDK_guts.h
+++ b/src/macro/DXDK_guts.h
@@ -13,14 +13,23 @@
 
 #ifdef FCODE 
 
-      DXDK_OO2DX = 1D0/(2D0*DXDK_DX)
-    
-      DXDK_DXDKXX = DXDK_OO2DX*(DXDK_KXX_IP - DXDK_KXX_IM)
-      DXDK_DXDKXY = DXDK_OO2DX*(DXDK_KXY_IP - DXDK_KXY_IM)
-      DXDK_DXDKXZ = DXDK_OO2DX*(DXDK_KXZ_IP - DXDK_KXZ_IM)
-      DXDK_DXDKYY = DXDK_OO2DX*(DXDK_KYY_IP - DXDK_KYY_IM)
-      DXDK_DXDKYZ = DXDK_OO2DX*(DXDK_KYZ_IP - DXDK_KYZ_IM)
-      DXDK_DXDKZZ = DXDK_OO2DX*(DXDK_KZZ_IP - DXDK_KZZ_IM)
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DXDK_DXDKXX = ADM_DX_2(kxx,i,j,k)
+        DXDK_DXDKXY = ADM_DX_2(kxy,i,j,k)
+        DXDK_DXDKXZ = ADM_DX_2(kxz,i,j,k)
+        DXDK_DXDKYY = ADM_DX_2(kyy,i,j,k)
+        DXDK_DXDKYZ = ADM_DX_2(kyz,i,j,k)
+        DXDK_DXDKZZ = ADM_DX_2(kzz,i,j,k)
+      else
+        DXDK_DXDKXX = ADM_DX_4(kxx,i,j,k)
+        DXDK_DXDKXY = ADM_DX_4(kxy,i,j,k)
+        DXDK_DXDKXZ = ADM_DX_4(kxz,i,j,k)
+        DXDK_DXDKYY = ADM_DX_4(kyy,i,j,k)
+        DXDK_DXDKYZ = ADM_DX_4(kyz,i,j,k)
+        DXDK_DXDKZZ = ADM_DX_4(kzz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DXXDG_guts.h b/src/macro/DXXDG_guts.h
index 997a62b..89dcedf 100644
--- a/src/macro/DXXDG_guts.h
+++ b/src/macro/DXXDG_guts.h
@@ -20,33 +20,55 @@
 
 #ifdef FCODE 
 
-      DXXDG_OODX2 = 1D0/(dx*dx)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DXXDG_FAC  = 0
       ELSE
-        DXXDG_FAC  = DXDG_PSI4*(4*DXXDG_DXXDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DXDG_DXDPSI_O_PSI)
+        DXXDG_FAC  = DXDG_PSI4*(4*DXXDG_DXXDPSI_O_PSI \
+                                + 12*DXDG_DXDPSI_O_PSI*DXDG_DXDPSI_O_PSI)
       ENDIF
       
       /* Now calculate the second derivatives */
-      DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC+DXXDG_FAC*DXDG_GXX+DXDG_PSI4 \
-       *DXXDG_OODX2*(DXDCG_GXX_IP-2*DXDG_GXX+DXDCG_GXX_IM)
 
-      DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC+DXXDG_FAC*DXDG_GXY+DXDG_PSI4 \
-      *DXXDG_OODX2*(DXDCG_GXY_IP-2*DXDG_GXY+DXDCG_GXY_IM)
+      if (spatial_order.eq.2) then
+        DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC + DXXDG_FAC*DXDG_GXX \
+                      + DXDG_PSI4*ADM_DXX_2(gxx,i,j,k)
 
-      DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC+DXXDG_FAC*DXDG_GXZ+DXDG_PSI4 \
-       *DXXDG_OODX2*(DXDCG_GXZ_IP-2*DXDG_GXZ+DXDCG_GXZ_IM)
+        DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC + DXXDG_FAC*DXDG_GXY \
+                      + DXDG_PSI4*ADM_DXX_2(gxy,i,j,k)
 
-      DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC+DXXDG_FAC*DXDG_GYY+DXDG_PSI4 \
-       *DXXDG_OODX2*(DXDCG_GYY_IP-2*DXDG_GYY+DXDCG_GYY_IM)
+        DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC + DXXDG_FAC*DXDG_GXZ \
+                      + DXDG_PSI4*ADM_DXX_2(gxz,i,j,k)
 
-      DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC+DXXDG_FAC*DXDG_GYZ+DXDG_PSI4 \
-       *DXXDG_OODX2*(DXDCG_GYZ_IP-2*DXDG_GYZ+DXDCG_GYZ_IM)
+        DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC + DXXDG_FAC*DXDG_GYY \
+                      + DXDG_PSI4*ADM_DXX_2(gyy,i,j,k)
 
-      DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC+DXXDG_FAC*DXDG_GZZ+DXDG_PSI4 \
-       *DXXDG_OODX2*(DXDCG_GZZ_IP-2*DXDG_GZZ+DXDCG_GZZ_IM)
+        DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC + DXXDG_FAC*DXDG_GYZ \
+                      + DXDG_PSI4*ADM_DXX_2(gyz,i,j,k)
+
+        DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC + DXXDG_FAC*DXDG_GZZ \
+                      + DXDG_PSI4*ADM_DXX_2(gzz,i,j,k)
+      else
+        DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC + DXXDG_FAC*DXDG_GXX \
+                      + DXDG_PSI4*ADM_DXX_4(gxx,i,j,k)
+
+        DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC + DXXDG_FAC*DXDG_GXY \
+                      + DXDG_PSI4*ADM_DXX_4(gxy,i,j,k)
+
+        DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC + DXXDG_FAC*DXDG_GXZ \
+                      + DXDG_PSI4*ADM_DXX_4(gxz,i,j,k)
+
+        DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC + DXXDG_FAC*DXDG_GYY \
+                      + DXDG_PSI4*ADM_DXX_4(gyy,i,j,k)
+
+        DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC + DXXDG_FAC*DXDG_GYZ \
+                      + DXDG_PSI4*ADM_DXX_4(gyz,i,j,k)
+
+        DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC + DXXDG_FAC*DXDG_GZZ \
+                      + DXDG_PSI4*ADM_DXX_4(gzz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DXYDG_guts.h b/src/macro/DXYDG_guts.h
index de1744e..f43de6c 100644
--- a/src/macro/DXYDG_guts.h
+++ b/src/macro/DXYDG_guts.h
@@ -22,39 +22,54 @@
 
 #ifdef FCODE 
 
-      DXYDG_OO4DXDY = 1D0/(4D0*dx*dy)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DXYDG_FAC  = 0
       ELSE
-        DXYDG_FAC  = DXDG_PSI4*(4*DXYDG_DXYDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DYDG_DYDPSI_O_PSI)
+        DXYDG_FAC  = DXDG_PSI4*(4*DXYDG_DXYDPSI_O_PSI + \
+                     12*DXDG_DXDPSI_O_PSI*DYDG_DYDPSI_O_PSI)
       ENDIF
       
       /* Now calculate the second deriatives */
-      DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC+DXDCG_DXDCGXX*DYDG_FAC+DXYDG_FAC*DXDG_GXX \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GXX_IPJP-DXYDG_GXX_IPJM-DXYDG_GXX_IMJP+DXYDG_GXX_IMJM);
+      if (spatial_order.eq.2) then
+        DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXY_2(gxx,i,j,k)
 
-      DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC+DXDCG_DXDCGXY*DYDG_FAC+DXYDG_FAC*DXDG_GXY \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GXY_IPJP-DXYDG_GXY_IPJM-DXYDG_GXY_IMJP+DXYDG_GXY_IMJM);
+        DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXY_2(gxy,i,j,k)
 
-      DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC+DXDCG_DXDCGXZ*DYDG_FAC+DXYDG_FAC*DXDG_GXZ \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GXZ_IPJP-DXYDG_GXZ_IPJM-DXYDG_GXZ_IMJP+DXYDG_GXZ_IMJM);
+        DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXY_2(gxz,i,j,k)
 
-      DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC+DXDCG_DXDCGYY*DYDG_FAC+DXYDG_FAC*DXDG_GYY \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GYY_IPJP-DXYDG_GYY_IPJM-DXYDG_GYY_IMJP+DXYDG_GYY_IMJM);
+        DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXY_2(gyy,i,j,k)
 
-      DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC+DXDCG_DXDCGYZ*DYDG_FAC+DXYDG_FAC*DXDG_GYZ \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GYZ_IPJP-DXYDG_GYZ_IPJM-DXYDG_GYZ_IMJP+DXYDG_GYZ_IMJM);
+        DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXY_2(gyz,i,j,k)
 
-      DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC+DXDCG_DXDCGZZ*DYDG_FAC+DXYDG_FAC*DXDG_GZZ \
-     +DXDG_PSI4*DXYDG_OO4DXDY* \
-     (DXYDG_GZZ_IPJP-DXYDG_GZZ_IPJM-DXYDG_GZZ_IMJP+DXYDG_GZZ_IMJM);
+        DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXY_2(gzz,i,j,k)
+      else
+        DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXY_4(gxx,i,j,k)
+
+        DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXY_4(gxy,i,j,k)
+
+        DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXY_4(gxz,i,j,k)
+
+        DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXY_4(gyy,i,j,k)
+
+        DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXY_4(gyz,i,j,k)
+
+        DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DYDG_FAC \
+                      + DXYDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXY_4(gzz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DXZDG_guts.h b/src/macro/DXZDG_guts.h
index 4978507..cdee51c 100644
--- a/src/macro/DXZDG_guts.h
+++ b/src/macro/DXZDG_guts.h
@@ -22,40 +22,56 @@
 
 #ifdef FCODE 
 
-      DXZDG_OO4DXDZ = 1D0/(4D0*dx*dz)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DXZDG_FAC  = 0
       ELSE
-        DXZDG_FAC  =  DXDG_PSI4*(4*DXZDG_DXZDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
+        DXZDG_FAC  =  DXDG_PSI4*(4*DXZDG_DXZDPSI_O_PSI \
+                   + 12*DXDG_DXDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
       ENDIF
 
       /* Now calculate the second deriatives */
-      DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC+DXDCG_DXDCGXX*DZDG_FAC+DXZDG_FAC*DXDG_GXX+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXX_IPKP-DXZDG_GXX_IPKM-DXZDG_GXX_IMKP+ \
-     DXZDG_GXX_IMKM)
+      if (spatial_order.eq.2) then 
+        DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXZ_2(gxx,i,j,k)
 
-      DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC+DXDCG_DXDCGXY*DZDG_FAC+DXZDG_FAC*DXDG_GXY+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXY_IPKP-DXZDG_GXY_IPKM-DXZDG_GXY_IMKP+ \
-     DXZDG_GXY_IMKM)
+        DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXZ_2(gxy,i,j,k)
 
-      DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC+DXDCG_DXDCGXZ*DZDG_FAC+DXZDG_FAC*DXDG_GXZ+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXZ_IPKP-DXZDG_GXZ_IPKM-DXZDG_GXZ_IMKP+ \
-     DXZDG_GXZ_IMKM)
+        DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXZ_2(gxz,i,j,k)
 
-      DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC+DXDCG_DXDCGYY*DZDG_FAC+DXZDG_FAC*DXDG_GYY+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GYY_IPKP-DXZDG_GYY_IPKM-DXZDG_GYY_IMKP+ \
-     DXZDG_GYY_IMKM)
+        DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXZ_2(gyy,i,j,k)
 
-      DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC+DXDCG_DXDCGYZ*DZDG_FAC+DXZDG_FAC*DXDG_GYZ+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GYZ_IPKP-DXZDG_GYZ_IPKM-DXZDG_GYZ_IMKP+ \
-     DXZDG_GYZ_IMKM)
+        DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXZ_2(gyz,i,j,k)
 
-      DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC+DXDCG_DXDCGZZ*DZDG_FAC+DXZDG_FAC*DXDG_GZZ+ \
-     DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GZZ_IPKP-DXZDG_GZZ_IPKM-DXZDG_GZZ_IMKP+ \
-     DXZDG_GZZ_IMKM)
+        DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXZ_2(gzz,i,j,k)
 
+      else
+
+        DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXZ_4(gxx,i,j,k)
+
+        DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXZ_4(gxy,i,j,k)
+
+        DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXZ_4(gxz,i,j,k)
+
+        DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXZ_4(gyy,i,j,k)
+
+        DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXZ_4(gyz,i,j,k)
+
+        DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DZDG_FAC \
+                      + DXZDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXZ_4(gzz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DYDB_guts.h b/src/macro/DYDB_guts.h
index 54cba6c..bd92f57 100644
--- a/src/macro/DYDB_guts.h
+++ b/src/macro/DYDB_guts.h
@@ -17,11 +17,17 @@
 
 #ifdef FCODE 
 
-      DYDB_OO2DY = 1.0D0/(2.0D0*DYDB_DY)
-
-      DYDB_DYDBX = DYDB_OO2DY*(DYDB_BX_JP - DYDB_BX_JM)
-      DYDB_DYDBY = DYDB_OO2DY*(DYDB_BY_JP - DYDB_BY_JM)
-      DYDB_DYDBZ = DYDB_OO2DY*(DYDB_BZ_JP - DYDB_BZ_JM)
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then
+        DYDB_DYDBX = ADM_DY_2(betax,i,j,k)
+        DYDB_DYDBY = ADM_DY_2(betay,i,j,k)
+        DYDB_DYDBZ = ADM_DY_2(betaz,i,j,k)
+      else
+        DYDB_DYDBX = ADM_DY_4(betax,i,j,k)
+        DYDB_DYDBY = ADM_DY_4(betay,i,j,k)
+        DYDB_DYDBZ = ADM_DY_4(betaz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DYDCG_guts.h b/src/macro/DYDCG_guts.h
index 01bced2..3c328d7 100644
--- a/src/macro/DYDCG_guts.h
+++ b/src/macro/DYDCG_guts.h
@@ -21,15 +21,23 @@
 
 #ifdef FCODE 
 
-      DYDCG_OO2DY = 1D0/(2D0*DYDCG_DY)
-    
-      DYDCG_DYDCGXX = DYDCG_OO2DY*(DYDCG_GXX_JP - DYDCG_GXX_JM)
-      DYDCG_DYDCGXY = DYDCG_OO2DY*(DYDCG_GXY_JP - DYDCG_GXY_JM)
-      DYDCG_DYDCGXZ = DYDCG_OO2DY*(DYDCG_GXZ_JP - DYDCG_GXZ_JM)
-      DYDCG_DYDCGYY = DYDCG_OO2DY*(DYDCG_GYY_JP - DYDCG_GYY_JM)
-      DYDCG_DYDCGYZ = DYDCG_OO2DY*(DYDCG_GYZ_JP - DYDCG_GYZ_JM)
-      DYDCG_DYDCGZZ = DYDCG_OO2DY*(DYDCG_GZZ_JP - DYDCG_GZZ_JM)
-
+#include "ADM_Derivative.h"
+
+      if (spatial_order.eq.2) then 
+        DYDCG_DYDCGXX = ADM_DY_2(gxx,i,j,k)
+        DYDCG_DYDCGXY = ADM_DY_2(gxy,i,j,k)
+        DYDCG_DYDCGXZ = ADM_DY_2(gxz,i,j,k)
+        DYDCG_DYDCGYY = ADM_DY_2(gyy,i,j,k)
+        DYDCG_DYDCGYZ = ADM_DY_2(gyz,i,j,k)
+        DYDCG_DYDCGZZ = ADM_DY_2(gzz,i,j,k)
+      else
+        DYDCG_DYDCGXX = ADM_DY_4(gxx,i,j,k)
+        DYDCG_DYDCGXY = ADM_DY_4(gxy,i,j,k)
+        DYDCG_DYDCGXZ = ADM_DY_4(gxz,i,j,k)
+        DYDCG_DYDCGYY = ADM_DY_4(gyy,i,j,k)
+        DYDCG_DYDCGYZ = ADM_DY_4(gyz,i,j,k)
+        DYDCG_DYDCGZZ = ADM_DY_4(gzz,i,j,k)
+      end if
 #endif
 
 
diff --git a/src/macro/DYDK_guts.h b/src/macro/DYDK_guts.h
index 26364b3..7dbfc0c 100644
--- a/src/macro/DYDK_guts.h
+++ b/src/macro/DYDK_guts.h
@@ -13,15 +13,23 @@
 
 #ifdef FCODE 
 
-      DYDK_OO2DY = 1D0/(2D0*DYDK_DY)
-    
-      DYDK_DYDKXX = DYDK_OO2DY*(DYDK_KXX_JP - DYDK_KXX_JM)
-      DYDK_DYDKXY = DYDK_OO2DY*(DYDK_KXY_JP - DYDK_KXY_JM)
-      DYDK_DYDKXZ = DYDK_OO2DY*(DYDK_KXZ_JP - DYDK_KXZ_JM)
-      DYDK_DYDKYY = DYDK_OO2DY*(DYDK_KYY_JP - DYDK_KYY_JM)
-      DYDK_DYDKYZ = DYDK_OO2DY*(DYDK_KYZ_JP - DYDK_KYZ_JM)
-      DYDK_DYDKZZ = DYDK_OO2DY*(DYDK_KZZ_JP - DYDK_KZZ_JM)
-
+#include "ADM_Derivative.h"
+  
+      if (spatial_order.eq.2) then 
+        DYDK_DYDKXX = ADM_DY_2(kxx,i,j,k)
+        DYDK_DYDKXY = ADM_DY_2(kxy,i,j,k)
+        DYDK_DYDKXZ = ADM_DY_2(kxz,i,j,k)
+        DYDK_DYDKYY = ADM_DY_2(kyy,i,j,k)
+        DYDK_DYDKYZ = ADM_DY_2(kyz,i,j,k)
+        DYDK_DYDKZZ = ADM_DY_2(kzz,i,j,k)
+      else
+        DYDK_DYDKXX = ADM_DY_4(kxx,i,j,k)
+        DYDK_DYDKXY = ADM_DY_4(kxy,i,j,k)
+        DYDK_DYDKXZ = ADM_DY_4(kxz,i,j,k)
+        DYDK_DYDKYY = ADM_DY_4(kyy,i,j,k)
+        DYDK_DYDKYZ = ADM_DY_4(kyz,i,j,k)
+        DYDK_DYDKZZ = ADM_DY_4(kzz,i,j,k)
+      end if
 #endif
 
 
diff --git a/src/macro/DYYDG_guts.h b/src/macro/DYYDG_guts.h
index 1896229..2dc495e 100644
--- a/src/macro/DYYDG_guts.h
+++ b/src/macro/DYYDG_guts.h
@@ -20,33 +20,54 @@
 
 #ifdef FCODE 
 
-      DYYDG_OODY2 = 1D0/(dy*dy)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DYYDG_FAC  = 0
       ELSE
-        DYYDG_FAC  = DYDG_PSI4*(4*DYYDG_DYYDPSI_O_PSI + 12*DYDG_DYDPSI_O_PSI*DYDG_DYDPSI_O_PSI)
+        DYYDG_FAC  = DYDG_PSI4*(4*DYYDG_DYYDPSI_O_PSI \
+                                + 12*DYDG_DYDPSI_O_PSI*DYDG_DYDPSI_O_PSI)
       ENDIF
 
       /* Now calculate the second deriatives */
-      DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC+DYYDG_FAC*DYDG_GXX+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GXX_JP-2*DYDG_GXX+DYDCG_GXX_JM)
+      if (spatial_order.eq.2) then 
+        DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC + DYYDG_FAC*DYDG_GXX \
+                      + DYDG_PSI4*ADM_DYY_2(gxx,i,j,k)
 
-      DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC+DYYDG_FAC*DYDG_GXY+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GXY_JP-2*DYDG_GXY+DYDCG_GXY_JM)
+        DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC + DYYDG_FAC*DYDG_GXY \
+                      + DYDG_PSI4*ADM_DYY_2(gxy,i,j,k)
 
-      DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC+DYYDG_FAC*DYDG_GXZ+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GXZ_JP-2*DYDG_GXZ+DYDCG_GXZ_JM)
+        DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC + DYYDG_FAC*DYDG_GXZ \
+                      + DYDG_PSI4*ADM_DYY_2(gxz,i,j,k)
 
-      DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC+DYYDG_FAC*DYDG_GYY+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GYY_JP-2*DYDG_GYY+DYDCG_GYY_JM)
+        DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC + DYYDG_FAC*DYDG_GYY \
+                      + DYDG_PSI4*ADM_DYY_2(gyy,i,j,k)
 
-      DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC+DYYDG_FAC*DYDG_GYZ+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GYZ_JP-2*DYDG_GYZ+DYDCG_GYZ_JM)
+        DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC + DYYDG_FAC*DYDG_GYZ \
+                      + DYDG_PSI4*ADM_DYY_2(gyz,i,j,k)
 
-      DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC+DYYDG_FAC*DYDG_GZZ+DYDG_PSI4 \
-     *DYYDG_OODY2*(DYDCG_GZZ_JP-2*DYDG_GZZ+DYDCG_GZZ_JM)
+        DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC + DYYDG_FAC*DYDG_GZZ \
+                      + DYDG_PSI4*ADM_DYY_2(gzz,i,j,k)
+      else
+        DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC + DYYDG_FAC*DYDG_GXX \
+                      + DYDG_PSI4*ADM_DYY_4(gxx,i,j,k)
+
+        DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC + DYYDG_FAC*DYDG_GXY \
+                      + DYDG_PSI4*ADM_DYY_4(gxy,i,j,k)
+
+        DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC + DYYDG_FAC*DYDG_GXZ \
+                      + DYDG_PSI4*ADM_DYY_4(gxz,i,j,k)
+
+        DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC + DYYDG_FAC*DYDG_GYY \
+                      + DYDG_PSI4*ADM_DYY_4(gyy,i,j,k)
+
+        DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC + DYYDG_FAC*DYDG_GYZ \
+                      + DYDG_PSI4*ADM_DYY_4(gyz,i,j,k)
+
+        DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC + DYYDG_FAC*DYDG_GZZ \
+                      + DYDG_PSI4*ADM_DYY_4(gzz,i,j,k)
+      end if
 
 #endif
 
diff --git a/src/macro/DYZDG_guts.h b/src/macro/DYZDG_guts.h
index 80347ff..65a513f 100644
--- a/src/macro/DYZDG_guts.h
+++ b/src/macro/DYZDG_guts.h
@@ -22,40 +22,54 @@
 
 #ifdef FCODE 
 
-      DYZDG_OO4DYDZ = 1D0/(4D0*dy*dz)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DYZDG_FAC  = 0
       ELSE
-        DYZDG_FAC  = DYDG_PSI4*(4*DYZDG_DYZDPSI_O_PSI + 12*DYDG_DYDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
+        DYZDG_FAC  = DYDG_PSI4*(4*DYZDG_DYZDPSI_O_PSI \
+                                + 12*DYDG_DYDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
       ENDIF
       
       /* Now calculate the second deriatives */
-      DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC+DYDCG_DYDCGXX*DZDG_FAC+DYZDG_FAC*DYDG_GXX \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXX_JPKP-DYZDG_GXX_JPKM-DYZDG_GXX_JMKP \
-        +DYZDG_GXX_JMKM)
+      if (spatial_order.eq.2) then
+        DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC + DYDCG_DYDCGXX*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXX + DYDG_PSI4*ADM_DYZ_2(gxx,i,j,k)
 
-      DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC+DYDCG_DYDCGXY*DZDG_FAC+DYZDG_FAC*DYDG_GXY \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXY_JPKP-DYZDG_GXY_JPKM-DYZDG_GXY_JMKP \
-        +DYZDG_GXY_JMKM)
+        DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC + DYDCG_DYDCGXY*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXY + DYDG_PSI4*ADM_DYZ_2(gxy,i,j,k)
 
-      DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC+DYDCG_DYDCGXZ*DZDG_FAC+DYZDG_FAC*DYDG_GXZ \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXZ_JPKP-DYZDG_GXZ_JPKM-DYZDG_GXZ_JMKP \
-        +DYZDG_GXZ_JMKM)
+        DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC + DYDCG_DYDCGXZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXZ + DYDG_PSI4*ADM_DYZ_2(gxz,i,j,k)
 
-      DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC+DYDCG_DYDCGYY*DZDG_FAC+DYZDG_FAC*DYDG_GYY \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GYY_JPKP-DYZDG_GYY_JPKM-DYZDG_GYY_JMKP \
-        +DYZDG_GYY_JMKM)
+        DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC + DYDCG_DYDCGYY*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GYY + DYDG_PSI4*ADM_DYZ_2(gyy,i,j,k)
 
-      DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC+DYDCG_DYDCGYZ*DZDG_FAC+DYZDG_FAC*DYDG_GYZ \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GYZ_JPKP-DYZDG_GYZ_JPKM-DYZDG_GYZ_JMKP \
-        +DYZDG_GYZ_JMKM)
+        DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC + DYDCG_DYDCGYZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GYZ + DYDG_PSI4*ADM_DYZ_2(gyz,i,j,k)
 
-      DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC+DYDCG_DYDCGZZ*DZDG_FAC+DYZDG_FAC*DYDG_GZZ \
-        +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GZZ_JPKP-DYZDG_GZZ_JPKM-DYZDG_GZZ_JMKP \
-        +DYZDG_GZZ_JMKM)
+        DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC + DYDCG_DYDCGZZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GZZ + DYDG_PSI4*ADM_DYZ_2(gzz,i,j,k)
+      else 
+        DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC + DYDCG_DYDCGXX*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXX + DYDG_PSI4*ADM_DYZ_4(gxx,i,j,k)
 
+        DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC + DYDCG_DYDCGXY*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXY + DYDG_PSI4*ADM_DYZ_4(gxy,i,j,k)
+
+        DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC + DYDCG_DYDCGXZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GXZ + DYDG_PSI4*ADM_DYZ_4(gxz,i,j,k)
+
+        DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC + DYDCG_DYDCGYY*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GYY + DYDG_PSI4*ADM_DYZ_4(gyy,i,j,k)
+
+        DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC + DYDCG_DYDCGYZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GYZ + DYDG_PSI4*ADM_DYZ_4(gyz,i,j,k)
+
+        DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC + DYDCG_DYDCGZZ*DZDG_FAC \
+                      + DYZDG_FAC*DYDG_GZZ + DYDG_PSI4*ADM_DYZ_4(gzz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/DZDB_guts.h b/src/macro/DZDB_guts.h
index 5797db7..4b01986 100644
--- a/src/macro/DZDB_guts.h
+++ b/src/macro/DZDB_guts.h
@@ -17,12 +17,16 @@
 
 #ifdef FCODE 
 
-      DZDB_OO2DZ = 1.0D0/(2.0D0*DZDB_DZ)
-
-      DZDB_DZDBX = DZDB_OO2DZ*(DZDB_BX_KP - DZDB_BX_KM)
-      DZDB_DZDBY = DZDB_OO2DZ*(DZDB_BY_KP - DZDB_BY_KM)
-      DZDB_DZDBZ = DZDB_OO2DZ*(DZDB_BZ_KP - DZDB_BZ_KM)
-
+#include "ADM_Derivative.h"
+      if (spatial_order.eq.2) then
+        DZDB_DZDBX = ADM_DZ_2(betax,i,j,k)
+        DZDB_DZDBY = ADM_DZ_2(betay,i,j,k)
+        DZDB_DZDBZ = ADM_DZ_2(betaz,i,j,k)
+      else
+        DZDB_DZDBX = ADM_DZ_4(betax,i,j,k)
+        DZDB_DZDBY = ADM_DZ_4(betay,i,j,k)
+        DZDB_DZDBZ = ADM_DZ_4(betaz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/DZDCG_guts.h b/src/macro/DZDCG_guts.h
index 82f3a09..08664d0 100644
--- a/src/macro/DZDCG_guts.h
+++ b/src/macro/DZDCG_guts.h
@@ -21,15 +21,22 @@
 
 #ifdef FCODE 
 
-      DZDCG_OO2DZ = 1D0/(2D0*DZDCG_DZ)
-    
-      DZDCG_DZDCGXX = DZDCG_OO2DZ*(DZDCG_GXX_KP - DZDCG_GXX_KM)
-      DZDCG_DZDCGXY = DZDCG_OO2DZ*(DZDCG_GXY_KP - DZDCG_GXY_KM)
-      DZDCG_DZDCGXZ = DZDCG_OO2DZ*(DZDCG_GXZ_KP - DZDCG_GXZ_KM)
-      DZDCG_DZDCGYY = DZDCG_OO2DZ*(DZDCG_GYY_KP - DZDCG_GYY_KM)
-      DZDCG_DZDCGYZ = DZDCG_OO2DZ*(DZDCG_GYZ_KP - DZDCG_GYZ_KM)
-      DZDCG_DZDCGZZ = DZDCG_OO2DZ*(DZDCG_GZZ_KP - DZDCG_GZZ_KM)
-
+#include "ADM_Derivative.h"
+      if (spatial_order.eq.2) then     
+        DZDCG_DZDCGXX = ADM_DZ_2(gxx,i,j,k)
+        DZDCG_DZDCGXY = ADM_DZ_2(gxy,i,j,k)
+        DZDCG_DZDCGXZ = ADM_DZ_2(gxz,i,j,k)
+        DZDCG_DZDCGYY = ADM_DZ_2(gyy,i,j,k)
+        DZDCG_DZDCGYZ = ADM_DZ_2(gyz,i,j,k)
+        DZDCG_DZDCGZZ = ADM_DZ_2(gzz,i,j,k)
+      else
+        DZDCG_DZDCGXX = ADM_DZ_4(gxx,i,j,k)
+        DZDCG_DZDCGXY = ADM_DZ_4(gxy,i,j,k)
+        DZDCG_DZDCGXZ = ADM_DZ_4(gxz,i,j,k)
+        DZDCG_DZDCGYY = ADM_DZ_4(gyy,i,j,k)
+        DZDCG_DZDCGYZ = ADM_DZ_4(gyz,i,j,k)
+        DZDCG_DZDCGZZ = ADM_DZ_4(gzz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/DZDK_guts.h b/src/macro/DZDK_guts.h
index 1e205ab..cc8e5e6 100644
--- a/src/macro/DZDK_guts.h
+++ b/src/macro/DZDK_guts.h
@@ -13,15 +13,23 @@
 
 #ifdef FCODE 
 
-      DZDK_OO2DZ = 1D0/(2D0*DZDK_DZ)
-    
-      DZDK_DZDKXX = DZDK_OO2DZ*(DZDK_KXX_KP - DZDK_KXX_KM)
-      DZDK_DZDKXY = DZDK_OO2DZ*(DZDK_KXY_KP - DZDK_KXY_KM)
-      DZDK_DZDKXZ = DZDK_OO2DZ*(DZDK_KXZ_KP - DZDK_KXZ_KM)
-      DZDK_DZDKYY = DZDK_OO2DZ*(DZDK_KYY_KP - DZDK_KYY_KM)
-      DZDK_DZDKYZ = DZDK_OO2DZ*(DZDK_KYZ_KP - DZDK_KYZ_KM)
-      DZDK_DZDKZZ = DZDK_OO2DZ*(DZDK_KZZ_KP - DZDK_KZZ_KM)
-
+#include "ADM_Derivative.h"
+  
+      if (spatial_order.eq.2) then   
+        DZDK_DZDKXX = ADM_DZ_2(kxx,i,j,k)
+        DZDK_DZDKXY = ADM_DZ_2(kxy,i,j,k)
+        DZDK_DZDKXZ = ADM_DZ_2(kxz,i,j,k)
+        DZDK_DZDKYY = ADM_DZ_2(kyy,i,j,k)
+        DZDK_DZDKYZ = ADM_DZ_2(kyz,i,j,k)
+        DZDK_DZDKZZ = ADM_DZ_2(kzz,i,j,k)
+      else
+        DZDK_DZDKXX = ADM_DZ_4(kxx,i,j,k)
+        DZDK_DZDKXY = ADM_DZ_4(kxy,i,j,k)
+        DZDK_DZDKXZ = ADM_DZ_4(kxz,i,j,k)
+        DZDK_DZDKYY = ADM_DZ_4(kyy,i,j,k)
+        DZDK_DZDKYZ = ADM_DZ_4(kyz,i,j,k)
+        DZDK_DZDKZZ = ADM_DZ_4(kzz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/DZZDG_guts.h b/src/macro/DZZDG_guts.h
index 601589b..853c12d 100644
--- a/src/macro/DZZDG_guts.h
+++ b/src/macro/DZZDG_guts.h
@@ -20,33 +20,53 @@
 
 #ifdef FCODE 
 
-      DZZDG_OODZ2 = 1D0/(dz*dz)
+#include "ADM_Derivative.h"
 
       /* Factor involving 2nd derivative of conformal factor */ 
       IF (conformal_state .eq. 0) THEN
         DZZDG_FAC  = 0
       ELSE
-        DZZDG_FAC  = DZDG_PSI4*(4*DZZDG_DZZDPSI_O_PSI + 12*DZDG_DZDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
+        DZZDG_FAC  = DZDG_PSI4*(4*DZZDG_DZZDPSI_O_PSI \
+                   + 12*DZDG_DZDPSI_O_PSI*DZDG_DZDPSI_O_PSI)
       ENDIF
 
-      DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC+DZZDG_FAC*DZDG_GXX+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GXX_KP-2*DZDG_GXX+DZDCG_GXX_KM)
+      if (spatial_order.eq.2) then
+        DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC + DZZDG_FAC*DZDG_GXX \
+        + DZDG_PSI4*ADM_DZZ_2(gxx,i,j,k)
 
-      DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC+DZZDG_FAC*DZDG_GXY+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GXY_KP-2*DZDG_GXY+DZDCG_GXY_KM)
+        DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC + DZZDG_FAC*DZDG_GXY \
+        + DZDG_PSI4*ADM_DZZ_2(gxy,i,j,k)
 
-      DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC+DZZDG_FAC*DZDG_GXZ+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GXZ_KP-2*DZDG_GXZ+DZDCG_GXZ_KM)
+        DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC + DZZDG_FAC*DZDG_GXZ \
+        + DZDG_PSI4*ADM_DZZ_2(gxz,i,j,k)
 
-      DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC+DZZDG_FAC*DZDG_GYY+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GYY_KP-2*DZDG_GYY+DZDCG_GYY_KM)
+        DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC + DZZDG_FAC*DZDG_GYY \
+        + DZDG_PSI4*ADM_DZZ_2(gyy,i,j,k)
 
-      DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC+DZZDG_FAC*DZDG_GYZ+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GYZ_KP-2*DZDG_GYZ+DZDCG_GYZ_KM)
+        DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC + DZZDG_FAC*DZDG_GYZ \
+        + DZDG_PSI4*ADM_DZZ_2(gyz,i,j,k)
 
-      DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC+DZZDG_FAC*DZDG_GZZ+DZDG_PSI4 \
-     *DZZDG_OODZ2*(DZDCG_GZZ_KP-2*DZDG_GZZ+DZDCG_GZZ_KM)
+        DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC + DZZDG_FAC*DZDG_GZZ \
+        + DZDG_PSI4*ADM_DZZ_2(gzz,i,j,k)
+      else
+        DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC + DZZDG_FAC*DZDG_GXX \
+        + DZDG_PSI4*ADM_DZZ_4(gxx,i,j,k)
 
+        DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC + DZZDG_FAC*DZDG_GXY \
+        + DZDG_PSI4*ADM_DZZ_4(gxy,i,j,k)
+
+        DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC + DZZDG_FAC*DZDG_GXZ \
+        + DZDG_PSI4*ADM_DZZ_4(gxz,i,j,k)
+
+        DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC + DZZDG_FAC*DZDG_GYY \
+        + DZDG_PSI4*ADM_DZZ_4(gyy,i,j,k)
+
+        DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC + DZZDG_FAC*DZDG_GYZ \
+        + DZDG_PSI4*ADM_DZZ_4(gyz,i,j,k)
+
+        DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC + DZZDG_FAC*DZDG_GZZ \
+        + DZDG_PSI4*ADM_DZZ_4(gzz,i,j,k)
+      end if
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/HAMADM_guts.h b/src/macro/HAMADM_guts.h
index 523d0d5..9c87ccb 100644
--- a/src/macro/HAMADM_guts.h
+++ b/src/macro/HAMADM_guts.h
@@ -21,14 +21,12 @@
 #ifdef FCODE 
 
       HAMADM_HAMADM = TRRICCI_TRRICCI - TRKK_TRKK + TRK_TRK**2
-      HAMADM_HAMADMABS = abs(TRRICCI_TRRICCI) + abs(TRKK_TRKK) + TRK_TRK**2
 
 #endif
 
 #ifdef CCODE
 
       HAMADM_HAMADM = TRRICCI_TRRICCI - TRKK_TRKK + TRK_TRK**2
-      HAMADM_HAMADMABS = fabs(TRRICCI_TRRICCI) + fabs(TRKK_TRKK) + TRK_TRK**2
 
 #endif
 
diff --git a/src/macro/LIEG_guts.h b/src/macro/LIEG_guts.h
index f36c812..3afd945 100644
--- a/src/macro/LIEG_guts.h
+++ b/src/macro/LIEG_guts.h
@@ -21,6 +21,8 @@
 
 #ifdef FCODE
 
+#include "ADM_Derivative.h"
+
 /*  Advection. */
 
       LIEG_LGXX = 0.0D0
@@ -43,63 +45,39 @@
 
       else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BX \
-              *(gxx(i+1,j,k) - gxx(i,j,k))/dx
-         LIEG_LGYY = LIEG_LGYY + LIEG_BX \
-              *(gyy(i+1,j,k) - gyy(i,j,k))/dx
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BX \
-              *(gzz(i+1,j,k) - gzz(i,j,k))/dx
-         LIEG_LGXY = LIEG_LGXY + LIEG_BX \
-              *(gxy(i+1,j,k) - gxy(i,j,k))/dx
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BX \
-              *(gxz(i+1,j,k) - gxz(i,j,k))/dx
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BX \
-              *(gyz(i+1,j,k) - gyz(i,j,k))/dx
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_P1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_P1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_P1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_P1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_P1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_P1(gyz,i,j,k)
 
       else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BX \
-              *(gxx(i,j,k) - gxx(i-1,j,k))/dx
-         LIEG_LGYY = LIEG_LGYY + LIEG_BX \
-              *(gyy(i,j,k) - gyy(i-1,j,k))/dx
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BX \
-              *(gzz(i,j,k) - gzz(i-1,j,k))/dx
-         LIEG_LGXY = LIEG_LGXY + LIEG_BX \
-              *(gxy(i,j,k) - gxy(i-1,j,k))/dx
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BX \
-              *(gxz(i,j,k) - gxz(i-1,j,k))/dx
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BX \
-              *(gyz(i,j,k) - gyz(i-1,j,k))/dx
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_M1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_M1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_M1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_M1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_M1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_M1(gyz,i,j,k)
 
       else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i+1,j,k) + gxx(i+2,j,k))
-         LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i+1,j,k) + gyy(i+2,j,k))
-         LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i+1,j,k) + gzz(i+2,j,k))
-         LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i+1,j,k) + gxy(i+2,j,k))
-         LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i+1,j,k) + gxz(i+2,j,k))
-         LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i+1,j,k) + gyz(i+2,j,k))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_P2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_P2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_P2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_P2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_P2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_P2(gyz,i,j,k)
 
       else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i-1,j,k) + gxx(i-2,j,k))
-         LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i-1,j,k) + gyy(i-2,j,k))
-         LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i-1,j,k) + gzz(i-2,j,k))
-         LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i-1,j,k) + gxy(i-2,j,k))
-         LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i-1,j,k) + gxz(i-2,j,k))
-         LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BX/dx \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i-1,j,k) + gyz(i-2,j,k))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_M2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_M2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_M2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_M2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_M2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_M2(gyz,i,j,k)
 
       end if
 
@@ -117,63 +95,39 @@
 
       else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BY \
-              *(gxx(i,j+1,k) - gxx(i,j,k))/dy
-         LIEG_LGYY = LIEG_LGYY + LIEG_BY \
-              *(gyy(i,j+1,k) - gyy(i,j,k))/dy
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BY \
-              *(gzz(i,j+1,k) - gzz(i,j,k))/dy
-         LIEG_LGXY = LIEG_LGXY + LIEG_BY \
-              *(gxy(i,j+1,k) - gxy(i,j,k))/dy
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BY \
-              *(gxz(i,j+1,k) - gxz(i,j,k))/dy
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BY \
-              *(gyz(i,j+1,k) - gyz(i,j,k))/dy
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_P1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_P1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_P1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_P1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_P1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_P1(gyz,i,j,k)
 
       else if (admmacros_advectiony .eq. -1) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BY \
-              *(gxx(i,j,k) - gxx(i,j-1,k))/dy
-         LIEG_LGYY = LIEG_LGYY + LIEG_BY \
-              *(gyy(i,j,k) - gyy(i,j-1,k))/dy
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BY \
-              *(gzz(i,j,k) - gzz(i,j-1,k))/dy
-         LIEG_LGXY = LIEG_LGXY + LIEG_BY \
-              *(gxy(i,j,k) - gxy(i,j-1,k))/dy
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BY \
-              *(gxz(i,j,k) - gxz(i,j-1,k))/dy
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BY \
-              *(gyz(i,j,k) - gyz(i,j-1,k))/dy
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_M1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_M1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_M1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_M1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_M1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_M1(gyz,i,j,k)
 
       else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j+1,k) + gxx(i,j+2,k))
-         LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j+1,k) + gyy(i,j+2,k))
-         LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j+1,k) + gzz(i,j+2,k))
-         LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j+1,k) + gxy(i,j+2,k))
-         LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j+1,k) + gxz(i,j+2,k))
-         LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j+1,k) + gyz(i,j+2,k))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_P2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_P2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_P2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_P2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_P2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_P2(gyz,i,j,k)
 
       else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j-1,k) + gxx(i,j-2,k))
-         LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j-1,k) + gyy(i,j-2,k))
-         LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j-1,k) + gzz(i,j-2,k))
-         LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j-1,k) + gxy(i,j-2,k))
-         LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j-1,k) + gxz(i,j-2,k))
-         LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BY/dy \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j-1,k) + gyz(i,j-2,k))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_M2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_M2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_M2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_M2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_M2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_M2(gyz,i,j,k)
 
       end if
 
@@ -191,73 +145,39 @@
 
       else if ((admmacros_advectionz .eq. 1)) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BZ \
-              *(gxx(i,j,k+1) - gxx(i,j,k))/dz
-         LIEG_LGYY = LIEG_LGYY + LIEG_BZ \
-              *(gyy(i,j,k+1) - gyy(i,j,k))/dz
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BZ \
-              *(gzz(i,j,k+1) - gzz(i,j,k))/dz
-         LIEG_LGXY = LIEG_LGXY + LIEG_BZ \
-              *(gxy(i,j,k+1) - gxy(i,j,k))/dz
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BZ \
-              *(gxz(i,j,k+1) - gxz(i,j,k))/dz
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BZ \
-              *(gyz(i,j,k+1) - gyz(i,j,k))/dz
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_P1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_P1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_P1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_P1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_P1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_P1(gyz,i,j,k)
 
       else if ((admmacros_advectionz .eq. -1)) then /* UPWIND1 */
 
-         LIEG_LGXX = LIEG_LGXX + LIEG_BZ \
-              *(gxx(i,j,k) - gxx(i,j,k-1))/dz
-         LIEG_LGYY = LIEG_LGYY + LIEG_BZ \
-              *(gyy(i,j,k) - gyy(i,j,k-1))/dz
-         LIEG_LGZZ = LIEG_LGZZ + LIEG_BZ \
-              *(gzz(i,j,k) - gzz(i,j,k-1))/dz
-         LIEG_LGXY = LIEG_LGXY + LIEG_BZ \
-              *(gxy(i,j,k) - gxy(i,j,k-1))/dz
-         LIEG_LGXZ = LIEG_LGXZ + LIEG_BZ \
-              *(gxz(i,j,k) - gxz(i,j,k-1))/dz
-         LIEG_LGYZ = LIEG_LGYZ + LIEG_BZ \
-              *(gyz(i,j,k) - gyz(i,j,k-1))/dz
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_M1(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_M1(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_M1(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_M1(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_M1(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_M1(gyz,i,j,k)
 
       else if (admmacros_advectionz .eq. 2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j,k+1) + gxx(i,j,k+2))
-
-         LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j,k+1) + gyy(i,j,k+2))
-
-         LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j,k+1) + gzz(i,j,k+2))
-
-         LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j,k+1) + gxy(i,j,k+2))
-
-         LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j,k+1) + gxz(i,j,k+2))
-
-         LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j,k+1) + gyz(i,j,k+2))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_P2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_P2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_P2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_P2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_P2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_P2(gyz,i,j,k)
 
       else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */
 
-         LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j,k-1) + gxx(i,j,k-2))
-
-         LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j,k-1) + gyy(i,j,k-2))
-
-         LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j,k-1) + gzz(i,j,k-2))
-
-         LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j,k-1) + gxy(i,j,k-2))
-
-         LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j,k-1) + gxz(i,j,k-2))
-
-         LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BZ/dz \
-              *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j,k-1) + gyz(i,j,k-2))
+         LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_M2(gxx,i,j,k)
+         LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_M2(gyy,i,j,k)
+         LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_M2(gzz,i,j,k)
+         LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_M2(gxy,i,j,k)
+         LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_M2(gxz,i,j,k)
+         LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_M2(gyz,i,j,k)
 
       end if 
 
diff --git a/src/macro/LIEK_guts.h b/src/macro/LIEK_guts.h
index b18b3f9..bdb6141 100644
--- a/src/macro/LIEK_guts.h
+++ b/src/macro/LIEK_guts.h
@@ -22,6 +22,8 @@
 
 #ifdef FCODE
 
+#include "ADM_Derivative.h"
+
 /*     Advection. */
 
       LIEK_LKXX = 0.0D0
@@ -34,96 +36,47 @@
       if (admmacros_advectionx .eq. 0) then /* CENTER X */
          
          LIEK_LKXX = LIEK_LKXX + DXDK_DXDKXX*LIEK_BX
-
          LIEK_LKYY = LIEK_LKYY + DXDK_DXDKYY*LIEK_BX
-
          LIEK_LKZZ = LIEK_LKZZ + DXDK_DXDKZZ*LIEK_BX
-
          LIEK_LKXY = LIEK_LKXY + DXDK_DXDKXY*LIEK_BX
-
          LIEK_LKXZ = LIEK_LKXZ + DXDK_DXDKXZ*LIEK_BX
-
          LIEK_LKYZ = LIEK_LKYZ + DXDK_DXDKYZ*LIEK_BX
 
       else if (admmacros_advectionx .eq. 1) then /* UPWIND1 X */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BX \
-         *(kxx(i+1,j,k) - kxx(i,j,k))/dx
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BX \
-         *(kyy(i+1,j,k) - kyy(i,j,k))/dx
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BX \
-         *(kzz(i+1,j,k) - kzz(i,j,k))/dx
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BX \
-         *(kxy(i+1,j,k) - kxy(i,j,k))/dx
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BX \
-         *(kxz(i+1,j,k) - kxz(i,j,k))/dx
-
-         LIEK_LKYZ = LIEK_LKYZ + LIEK_BX \
-         *(kyz(i+1,j,k) - kyz(i,j,k))/dx
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_P1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_P1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_P1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_P1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_P1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_P1(kyz,i,j,k)
 
       else if (admmacros_advectionx .eq. -1) then /* UPWIND1 X */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BX \
-         *(kxx(i,j,k) - kxx(i-1,j,k))/dx
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BX \
-         *(kyy(i,j,k) - kyy(i-1,j,k))/dx
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BX \
-         *(kzz(i,j,k) - kzz(i-1,j,k))/dx
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BX \
-         *(kxy(i,j,k) - kxy(i-1,j,k))/dx
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BX \
-         *(kxz(i,j,k) - kxz(i-1,j,k))/dx
-
-         LIEK_LKYZ = LIEK_LKYZ + LIEK_BX \
-         *(kyz(i,j,k) - kyz(i-1,j,k))/dx
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_M1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_M1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_M1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_M1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_M1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_M1(kyz,i,j,k)
 
       else if (admmacros_advectionx .eq. 2) then /* UPWIND2 X */
 
-         LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i+1,j,k) + kxx(i+2,j,k))
-
-         LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i+1,j,k) + kyy(i+2,j,k))
-
-         LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i+1,j,k) + kzz(i+2,j,k))
-
-         LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i+1,j,k) + kxy(i+2,j,k))
-
-         LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i+1,j,k) + kxz(i+2,j,k))
-
-         LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i+1,j,k) + kyz(i+2,j,k))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_P2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_P2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_P2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_P2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_P2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_P2(kyz,i,j,k)
 
       else if (admmacros_advectionx .eq. -2) then /* UPWIND2 X */
 
-         LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i-1,j,k) + kxx(i-2,j,k))
-
-         LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i-1,j,k) + kyy(i-2,j,k))
-
-         LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i-1,j,k) + kzz(i-2,j,k))
-
-         LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i-1,j,k) + kxy(i-2,j,k))
-
-         LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i-1,j,k) + kxz(i-2,j,k))
-
-         LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BX/dx \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i-1,j,k) + kyz(i-2,j,k))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_M2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_M2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_M2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_M2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_M2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_M2(kyz,i,j,k)
 
       end if
 
@@ -131,96 +84,47 @@
       if (admmacros_advectiony .eq. 0) then /* CENTER Y */
 
          LIEK_LKXX = LIEK_LKXX + DYDK_DYDKXX*LIEK_BY
-
          LIEK_LKYY = LIEK_LKYY + DYDK_DYDKYY*LIEK_BY
-
          LIEK_LKZZ = LIEK_LKZZ + DYDK_DYDKZZ*LIEK_BY
-
          LIEK_LKXY = LIEK_LKXY + DYDK_DYDKXY*LIEK_BY
-
          LIEK_LKXZ = LIEK_LKXZ + DYDK_DYDKXZ*LIEK_BY
-
          LIEK_LKYZ = LIEK_LKYZ + DYDK_DYDKYZ*LIEK_BY
 
       else if (admmacros_advectiony .eq. 1) then /* UPWIND1 Y */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BY \
-         *(kxx(i,j+1,k) - kxx(i,j,k))/dy
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BY \
-         *(kyy(i,j+1,k) - kyy(i,j,k))/dy
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BY \
-         *(kzz(i,j+1,k) - kzz(i,j,k))/dy
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BY \
-         *(kxy(i,j+1,k) - kxy(i,j,k))/dy
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BY \
-         *(kxz(i,j+1,k) - kxz(i,j,k))/dy
-
-         LIEK_LKYZ = LIEK_LKYZ + LIEK_BY \
-         *(kyz(i,j+1,k) - kyz(i,j,k))/dy
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P1(kyz,i,j,k)
 
       else if (admmacros_advectiony .eq. -1) then /* UPWIND1 Y */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BY \
-         *(kxx(i,j,k) - kxx(i,j-1,k))/dy
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BY \
-         *(kyy(i,j,k) - kyy(i,j-1,k))/dy
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BY \
-         *(kzz(i,j,k) - kzz(i,j-1,k))/dy
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BY \
-         *(kxy(i,j,k) - kxy(i,j-1,k))/dy
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BY \
-         *(kxz(i,j,k) - kxz(i,j-1,k))/dy
-
-         LIEK_LKYZ = LIEK_LKYZ + LIEK_BY \
-         *(kyz(i,j,k) - kyz(i,j-1,k))/dy
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M1(kyz,i,j,k)
 
       else if (admmacros_advectiony .eq. 2) then /* UPWIND2 Y */
 
-         LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j+1,k) + kxx(i,j+2,k))
-
-         LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j+1,k) + kyy(i,j+2,k))
-
-         LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j+1,k) + kzz(i,j+2,k))
-
-         LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j+1,k) + kxy(i,j+2,k))
-
-         LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j+1,k) + kxz(i,j+2,k))
-
-         LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j+1,k) + kyz(i,j+2,k))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P2(kyz,i,j,k)
 
       else if (admmacros_advectiony .eq. -2) then /* UPWIND2 Y */
 
-         LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j-1,k) + kxx(i,j-2,k))
-
-         LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j-1,k) + kyy(i,j-2,k))
-
-         LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j-1,k) + kzz(i,j-2,k))
-
-         LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j-1,k) + kxy(i,j-2,k))
-
-         LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j-1,k) + kxz(i,j-2,k))
-
-         LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BY/dy \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j-1,k) + kyz(i,j-2,k))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M2(kyz,i,j,k)
 
       end if
 
@@ -228,96 +132,50 @@
       if (admmacros_advectionz .eq. 0) then /* CENTER Z */
 
          LIEK_LKXX = LIEK_LKXX + DZDK_DZDKXX*LIEK_BZ
-
          LIEK_LKYY = LIEK_LKYY + DZDK_DZDKYY*LIEK_BZ
-
          LIEK_LKZZ = LIEK_LKZZ + DZDK_DZDKZZ*LIEK_BZ
-
          LIEK_LKXY = LIEK_LKXY + DZDK_DZDKXY*LIEK_BZ
-
          LIEK_LKXZ = LIEK_LKXZ + DZDK_DZDKXZ*LIEK_BZ
-
          LIEK_LKYZ = LIEK_LKYZ + DZDK_DZDKYZ*LIEK_BZ
 
       else if (admmacros_advectionz .eq. 1) then /* UPWIND1 Z */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BZ \
-         *(kxx(i,j,k+1) - kxx(i,j,k))/dz
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BZ \
-         *(kyy(i,j,k+1) - kyy(i,j,k))/dz
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BZ \
-         *(kzz(i,j,k+1) - kzz(i,j,k))/dz
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BZ \
-         *(kxy(i,j,k+1) - kxy(i,j,k))/dz
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BZ \
-         *(kxz(i,j,k+1) - kxz(i,j,k))/dz
-
-         LIEK_LKYZ = LIEK_LKYZ + LIEK_BZ \
-         *(kyz(i,j,k+1) - kyz(i,j,k))/dz
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P1(kyz,i,j,k)
 
       else if (admmacros_advectionz .eq. -1) then /* UPWIND1 Z */
 
-         LIEK_LKXX = LIEK_LKXX + LIEK_BZ \
-         *(kxx(i,j,k) - kxx(i,j,k-1))/dz
-
-         LIEK_LKYY = LIEK_LKYY + LIEK_BZ \
-         *(kyy(i,j,k) - kyy(i,j,k-1))/dz
-
-         LIEK_LKZZ = LIEK_LKZZ + LIEK_BZ \
-         *(kzz(i,j,k) - kzz(i,j,k-1))/dz
-
-         LIEK_LKXY = LIEK_LKXY + LIEK_BZ \
-         *(kxy(i,j,k) - kxy(i,j,k-1))/dz
-
-         LIEK_LKXZ = LIEK_LKXZ + LIEK_BZ \
-         *(kxz(i,j,k) - kxz(i,j,k-1))/dz
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M1(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M1(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M1(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M1(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M1(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M1(kyz,i,j,k)
 
       else if (admmacros_advectionz .eq. 2) then /* UPWIND2 Z */
 
-         LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j,k+1) + kxx(i,j,k+2))
-
-         LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j,k+1) + kyy(i,j,k+2))
-
-         LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j,k+1) + kzz(i,j,k+2))
-
-         LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j,k+1) + kxy(i,j,k+2))
-
-         LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j,k+1) + kxz(i,j,k+2))
-
-         LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j,k+1) + kyz(i,j,k+2))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P2(kyz,i,j,k)
 
       else if (admmacros_advectionz .eq. -2) then /* UPWIND2 Z */
 
-         LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j,k-1) + kxx(i,j,k-2))
-
-         LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j,k-1) + kyy(i,j,k-2))
+         LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M2(kxx,i,j,k)
+         LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M2(kyy,i,j,k)
+         LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M2(kzz,i,j,k)
+         LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M2(kxy,i,j,k)
+         LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M2(kxz,i,j,k)
+         LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M2(kyz,i,j,k)
 
-         LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j,k-1) + kzz(i,j,k-2))
-
-         LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j,k-1) + kxy(i,j,k-2))
-
-         LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j,k-1) + kxz(i,j,k-2))
-
-         LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BZ/dz \
-         *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j,k-1) + kyz(i,j,k-2))
       end if
 
-
       /*    Extra terms in the Lie derivative. */
 
       LIEK_LKXX = LIEK_LKXX + 2.0D0*(DXDB_DXDBX*LIEK_KXX \