Adding possibility of using upwind and second order upwind for calculation

of advection terms in the Lie derivatives.


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

2 files changed, 704 insertions, 90 deletions

diff --git a/src/macro/LIEG_guts.h b/src/macro/LIEG_guts.h
index d188c3d..4bdafcb 100644
--- a/src/macro/LIEG_guts.h
+++ b/src/macro/LIEG_guts.h
@@ -1,10 +1,14 @@
 /*@@
   @header   LIEG_guts.h
-  @date     Jun 98
-  @author   Gabrielle Allen
+  @date     July 2000
+  @author   Gabrielle Allen, Miguel Alcubierre
   @desc
+
   Macro to calculate the Lie derivative of the lower 
-  conformal metric
+  conformal metric.  Notice that the advection term
+  can be calculated in several different ways.
+
+  IMPORTANT:  THE C VERSION ONLY HAS CENTERED DIFFERENCES!
 
   @enddesc
 @@*/
@@ -17,63 +21,365 @@
 
 #ifdef FCODE
 
-      LIEG_LGXX = DXDCG_DXDCGXX*LIEG_BX + DYDCG_DYDCGXX*LIEG_BY +
-     & DZDCG_DZDCGXX*LIEG_BZ + 2*(DXDB_DXDBX*LIEG_GXX + DXDB_DXDBY*LIEG_GXY +
-     & DXDB_DXDBZ*LIEG_GXZ)
+c     Advection.
+
+      LIEG_LGXX = 0.0D0
+      LIEG_LGYY = 0.0D0
+      LIEG_LGZZ = 0.0D0
+      LIEG_LGXY = 0.0D0
+      LIEG_LGXZ = 0.0D0
+      LIEG_LGYZ = 0.0D0
+
+      if (CCTK_Equals(advection,'center').eq.1) then
+
+         LIEG_LGXX = DXDCG_DXDCGXX*LIEG_BX + DYDCG_DYDCGXX*LIEG_BY
+     &      + DZDCG_DZDCGXX*LIEG_BZ
+
+         LIEG_LGYY = DXDCG_DXDCGYY*LIEG_BX + DYDCG_DYDCGYY*LIEG_BY
+     &      + DZDCG_DZDCGYY*LIEG_BZ
+
+         LIEG_LGZZ = DXDCG_DXDCGZZ*LIEG_BX + DYDCG_DYDCGZZ*LIEG_BY
+     &      + DZDCG_DZDCGZZ*LIEG_BZ
+
+         LIEG_LGXY = DXDCG_DXDCGXY*LIEG_BX + DYDCG_DYDCGXY*LIEG_BY
+     &      + DZDCG_DZDCGXY*LIEG_BZ
+
+         LIEG_LGXZ = DXDCG_DXDCGXZ*LIEG_BX + DYDCG_DYDCGXZ*LIEG_BY
+     &      + DZDCG_DZDCGXZ*LIEG_BZ
+
+         LIEG_LGYZ = DXDCG_DXDCGYZ*LIEG_BX + DYDCG_DYDCGYZ*LIEG_BY
+     &      + DZDCG_DZDCGYZ*LIEG_BZ
+
+      else if ((CCTK_Equals(advection,'upwind1').eq.1).or.
+     &   (i.eq.2).or.(i.eq.cctk_lsh(1)-1).or.
+     &   (j.eq.2).or.(j.eq.cctk_lsh(2)-1).or.
+     &   (k.eq.2).or.(k.eq.cctk_lsh(3)-1)) then
+
+         if (LIEG_BX.gt.0.0D0) then
+
+            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
+
+         else
+
+            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
+
+         end if
+
+         if (LIEG_BY.gt.0.0D0) then
+
+            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
+
+         else
+
+            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
+
+         end if
+
+         if (LIEG_BZ.gt.0.0D0) then
+
+            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
+
+         else
+
+            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
+
+         end if
+
+      else if (CCTK_Equals(advection,'upwind2').eq.1) then
 
-      LIEG_LGXY = DXDCG_DXDCGXY*LIEG_BX + DYDCG_DYDCGXY*LIEG_BY +
-     & DZDCG_DZDCGXY*LIEG_BZ + DYDB_DYDBX*LIEG_GXX + (DXDB_DXDBX +
-     & DYDB_DYDBY)*LIEG_GXY + DYDB_DYDBZ*LIEG_GXZ + DXDB_DXDBY*LIEG_GYY +
-     & DXDB_DXDBZ*LIEG_GYZ
+         if (LIEG_BX.gt.0.0D0) then
 
-      LIEG_LGXZ = DXDCG_DXDCGXZ*LIEG_BX + DYDCG_DYDCGXZ*LIEG_BY +
-     & DZDCG_DZDCGXZ*LIEG_BZ + DZDB_DZDBX*LIEG_GXX + DZDB_DZDBY*LIEG_GXY +
-     & (DXDB_DXDBX + DZDB_DZDBZ)*LIEG_GXZ + DXDB_DXDBY*LIEG_GYZ +
-     & DXDB_DXDBZ*LIEG_GZZ
+            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 = DXDCG_DXDCGYY*LIEG_BX +DYDCG_DYDCGYY*LIEG_BY +
-     & DZDCG_DZDCGYY*LIEG_BZ + 2*(DYDB_DYDBX*LIEG_GXY + DYDB_DYDBY*LIEG_GYY +
-     & DYDB_DYDBZ*LIEG_GYZ)
+            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_LGYZ = DXDCG_DXDCGYZ*LIEG_BX + DYDCG_DYDCGYZ*LIEG_BY +
-     & DZDCG_DZDCGYZ*LIEG_BZ + DZDB_DZDBX*LIEG_GXY + DYDB_DYDBX*LIEG_GXZ +
-     & DZDB_DZDBY*LIEG_GYY + (DYDB_DYDBY + DZDB_DZDBZ)*LIEG_GYZ +
-     & DYDB_DYDBZ*LIEG_GZZ
+            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_LGZZ = DXDCG_DXDCGZZ*LIEG_BX + DYDCG_DYDCGZZ*LIEG_BY +
-     &DZDCG_DZDCGZZ*LIEG_BZ + 2*(DZDB_DZDBX*LIEG_GXZ + DZDB_DZDBY*LIEG_GYZ +
-     &DZDB_DZDBZ*LIEG_GZZ)
+            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))
+
+         else
+
+            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))
+
+         end if
+
+         if (LIEG_BY.gt.0.0D0) then
+
+            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))
+
+         else
+
+            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))
+
+         end if
+
+         if (LIEG_BZ.gt.0.0D0) then
+
+            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))
+
+         else
+
+            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))
+
+         end if
+
+      end if
+
+c     Extra terms in the Lie derivative.
+
+      LIEG_LGXX = LIEG_LGXX + 2.0D0*(DXDB_DXDBX*LIEG_GXX
+     &   + DXDB_DXDBY*LIEG_GXY + DXDB_DXDBZ*LIEG_GXZ)
+
+      LIEG_LGYY = LIEG_LGYY + 2.0D0*(DYDB_DYDBX*LIEG_GXY
+     &   + DYDB_DYDBY*LIEG_GYY + DYDB_DYDBZ*LIEG_GYZ)
+
+      LIEG_LGZZ = LIEG_LGZZ + 2.0D0*(DZDB_DZDBX*LIEG_GXZ
+     &   + DZDB_DZDBY*LIEG_GYZ + DZDB_DZDBZ*LIEG_GZZ)
+
+      LIEG_LGXY = LIEG_LGXY + DYDB_DYDBX*LIEG_GXX + DXDB_DXDBY*LIEG_GYY
+     &   + (DXDB_DXDBX + DYDB_DYDBY)*LIEG_GXY
+     &   + DYDB_DYDBZ*LIEG_GXZ + DXDB_DXDBZ*LIEG_GYZ
+
+      LIEG_LGXZ = LIEG_LGXZ + DZDB_DZDBX*LIEG_GXX + DXDB_DXDBZ*LIEG_GZZ
+     &   + (DXDB_DXDBX + DZDB_DZDBZ)*LIEG_GXZ
+     &   + DZDB_DZDBY*LIEG_GXY + DXDB_DXDBY*LIEG_GYZ
+
+      LIEG_LGYZ = LIEG_LGYZ + DZDB_DZDBY*LIEG_GYY + DYDB_DYDBZ*LIEG_GZZ
+     &   + (DYDB_DYDBY + DZDB_DZDBZ)*LIEG_GYZ
+     &   + DZDB_DZDBX*LIEG_GXY + DYDB_DYDBX*LIEG_GXZ
 
 #endif
 
 #ifdef CCODE
 
-      LIEG_LGXX = DXDCG_DXDCGXX*LIEG_BX + DYDCG_DYDCGXX*LIEG_BY +
-      DZDCG_DZDCGXX*LIEG_BZ + 2*(DXDB_DXDBX*LIEG_GXX + DXDB_DXDBY*LIEG_GXY +
-      DXDB_DXDBZ*LIEG_GXZ) ;
+/* Advection */
+
+      LIEG_LGXX = DXDCG_DXDCGXX*LIEG_BX + DYDCG_DYDCGXX*LIEG_BY
+         + DZDCG_DZDCGXX*LIEG_BZ;
+
+      LIEG_LGYY = DXDCG_DXDCGYY*LIEG_BX + DYDCG_DYDCGYY*LIEG_BY
+         + DZDCG_DZDCGYY*LIEG_BZ;
+
+      LIEG_LGZZ = DXDCG_DXDCGZZ*LIEG_BX + DYDCG_DYDCGZZ*LIEG_BY
+         + DZDCG_DZDCGZZ*LIEG_BZ;
+
+      LIEG_LGXY = DXDCG_DXDCGXY*LIEG_BX + DYDCG_DYDCGXY*LIEG_BY
+         + DZDCG_DZDCGXY*LIEG_BZ;
+
+      LIEG_LGXZ = DXDCG_DXDCGXZ*LIEG_BX + DYDCG_DYDCGXZ*LIEG_BY
+         + DZDCG_DZDCGXZ*LIEG_BZ;
+
+      LIEG_LGYZ = DXDCG_DXDCGYZ*LIEG_BX + DYDCG_DYDCGYZ*LIEG_BY
+         + DZDCG_DZDCGYZ*LIEG_BZ;
+
+/* Extra terms in the Lie derivative */
+
+      LIEG_LGXX = LIEG_LGXX + 2*(DXDB_DXDBX*LIEG_GXX
+         + DXDB_DXDBY*LIEG_GXY + DXDB_DXDBZ*LIEG_GXZ);
 
-      LIEG_LGXY = DXDCG_DXDCGXY*LIEG_BX + DYDCG_DYDCGXY*LIEG_BY +
-      DZDCG_DZDCGXY*LIEG_BZ + DYDB_DYDBX*LIEG_GXX + (DXDB_DXDBX +
-      DYDB_DYDBY)*LIEG_GXY + DYDB_DYDBZ*LIEG_GXZ+ DXDB_DXDBY*LIEG_GYY +
-      DXDB_DXDBZ*LIEG_GYZ ;
+      LIEG_LGYY = LIEG_LGYY + 2*(DYDB_DYDBX*LIEG_GXY
+         + DYDB_DYDBY*LIEG_GYY + DYDB_DYDBZ*LIEG_GYZ);
 
-      LIEG_LGXZ = DXDCG_DXDCGXZ*LIEG_BX + DYDCG_DYDCGXZ*LIEG_BY +
-      DZDCG_DZDCGXZ*LIEG_BZ + DZDB_DZDBX*LIEG_GXX + DZDB_DZDBY*LIEG_GXY +
-      (DXDB_DXDBX + DZDB_DZDBZ)*LIEG_GXZ+ DXDB_DXDBY*LIEG_GYZ +
-      DXDB_DXDBZ*LIEG_GZZ ;
+      LIEG_LGZZ = LIEG_LGZZ + 2*(DZDB_DZDBX*LIEG_GXZ
+         + DZDB_DZDBY*LIEG_GYZ + DZDB_DZDBZ*LIEG_GZZ);
 
-      LIEG_LGYY = DXDCG_DXDCGYY*LIEG_BX + DYDCG_DYDCGYY*LIEG_BY +
-      DZDCG_DZDCGYY*LIEG_BZ + 2*(DYDB_DYDBX*LIEG_GXY + DYDB_DYDBY*LIEG_GYY +
-      DYDB_DYDBZ*LIEG_GYZ) ;
+      LIEG_LGXY = LIEG_LGXY + DYDB_DYDBX*LIEG_GXX + DXDB_DXDBY*LIEG_GYY
+         + (DXDB_DXDBX + DYDB_DYDBY)*LIEG_GXY
+         + DYDB_DYDBZ*LIEG_GXZ + DXDB_DXDBZ*LIEG_GYZ;
 
-      LIEG_LGYZ = DXDCG_DXDCGYZ*LIEG_BX + DYDCG_DYDCGYZ*LIEG_BY +
-      DZDCG_DZDCGYZ*LIEG_BZ + DZDB_DZDBX*LIEG_GXY + DYDB_DYDBX*LIEG_GXZ+
-      DZDB_DZDBY*LIEG_GYY + (DYDB_DYDBY + DZDB_DZDBZ)*LIEG_GYZ +
-      DYDB_DYDBZ*LIEG_GZZ ;
+      LIEG_LGXZ = LIEG_LGXZ + DZDB_DZDBX*LIEG_GXX + DXDB_DXDBZ*LIEG_GZZ
+         + (DXDB_DXDBX + DZDB_DZDBZ)*LIEG_GXZ
+         + DZDB_DZDBY*LIEG_GXY + DXDB_DXDBY*LIEG_GYZ;
 
-      LIEG_LGZZ = DXDCG_DXDCGZZ*LIEG_BX + DYDCG_DYDCGZZ*LIEG_BY +
-      DZDCG_DZDCGZZ*LIEG_BZ + 2*(DZDB_DZDBX*LIEG_GXZ + DZDB_DZDBY*LIEG_GYZ +
-      DZDB_DZDBZ*LIEG_GZZ) ;
+      LIEG_LGYZ = LIEG_LGYZ + DZDB_DZDBY*LIEG_GYY + DYDB_DYDBZ*LIEG_GZZ
+         + (DYDB_DYDBY + DZDB_DZDBZ)*LIEG_GYZ
+         + DZDB_DZDBX*LIEG_GXY + DYDB_DYDBX*LIEG_GXZ;
 
 #endif
 
diff --git a/src/macro/LIEK_guts.h b/src/macro/LIEK_guts.h
index 2c66513..29f634f 100644
--- a/src/macro/LIEK_guts.h
+++ b/src/macro/LIEK_guts.h
@@ -1,10 +1,15 @@
 /*@@
   @header   LIEK_guts.h
-  @date     Jul 98
-  @author   Gabrielle Allen
+  @date     Jul 2000
+  @author   Gabrielle Allen + Miguel Alcubierre
   @desc
+
   Macro to calculate the Lie derivative of the lower 
-  extrinsic curvature along the shift vector
+  extrinsic curvature along the shift vector.   Notice
+  that the advection term can be calculated in several
+  different ways.
+
+  IMPORTANT:  THE C VERSION ONLY HAS CENTERED DIFFERENCES!
 
   @enddesc
 @@*/
@@ -17,63 +22,366 @@
 
 #ifdef FCODE
 
-      LIEK_LKXX = DXDK_DXDKXX*LIEK_BX + DYDK_DYDKXX*LIEK_BY +
-     &DZDK_DZDKXX*LIEK_BZ + 2*(DXDB_DXDBX*LIEK_KXX + DXDB_DXDBY*LIEK_KXY +
-     &DXDB_DXDBZ*LIEK_KXZ)
+c     Advection.
+
+      LIEK_LKXX = 0.0D0
+      LIEK_LKYY = 0.0D0
+      LIEK_LKZZ = 0.0D0
+      LIEK_LKXY = 0.0D0
+      LIEK_LKXZ = 0.0D0
+      LIEK_LKYZ = 0.0D0
+
+      if (CCTK_Equals(advection,'center').eq.1) then
+
+         LIEK_LKXX = DXDK_DXDKXX*LIEK_BX + DYDK_DYDKXX*LIEK_BY
+     &      + DZDK_DZDKXX*LIEK_BZ
+
+         LIEK_LKYY = DXDK_DXDKYY*LIEK_BX + DYDK_DYDKYY*LIEK_BY
+     &      + DZDK_DZDKYY*LIEK_BZ
+
+         LIEK_LKZZ = DXDK_DXDKZZ*LIEK_BX + DYDK_DYDKZZ*LIEK_BY
+     &      + DZDK_DZDKZZ*LIEK_BZ
+
+         LIEK_LKXY = DXDK_DXDKXY*LIEK_BX + DYDK_DYDKXY*LIEK_BY
+     &      + DZDK_DZDKXY*LIEK_BZ
+
+         LIEK_LKXZ = DXDK_DXDKXZ*LIEK_BX + DYDK_DYDKXZ*LIEK_BY
+     &      + DZDK_DZDKXZ*LIEK_BZ
+
+         LIEK_LKYZ = DXDK_DXDKYZ*LIEK_BX + DYDK_DYDKYZ*LIEK_BY
+     &      + DZDK_DZDKYZ*LIEK_BZ
+
+      else if ((CCTK_Equals(advection,'upwind1').eq.1).or.
+     &   (i.eq.2).or.(i.eq.cctk_lsh(1)-1).or.
+     &   (j.eq.2).or.(j.eq.cctk_lsh(2)-1).or.
+     &   (k.eq.2).or.(k.eq.cctk_lsh(3)-1)) then
+
+         if (LIEK_BX.gt.0.0D0) then
+
+            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
+
+         else
+
+            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
+
+         end if
+
+         if (LIEK_BY.gt.0.0D0) then
+
+            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
+
+         else
+
+            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
+
+         end if
+
+         if (LIEK_BZ.gt.0.0D0) then
+
+            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
+
+         else
+
+            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_LKYZ = LIEK_LKYZ + LIEK_BZ
+     &         *(kyz(i,j,k) - kyz(i,j,k-1))/dz
+
+         end if
+
+      else if (CCTK_Equals(advection,'upwind2').eq.1) then
 
-      LIEK_LKXY = DXDK_DXDKXY*LIEK_BX + DYDK_DYDKXY*LIEK_BY +
-     &DZDK_DZDKXY*LIEK_BZ + DYDB_DYDBX*LIEK_KXX + (DXDB_DXDBX +
-     &DYDB_DYDBY)*LIEK_KXY + DYDB_DYDBZ*LIEK_KXZ + DXDB_DXDBY*LIEK_KYY +
-     &DXDB_DXDBZ*LIEK_KYZ
+         if (LIEK_BX.gt.0.0D0) then
 
-      LIEK_LKXZ = DXDK_DXDKXZ*LIEK_BX + DYDK_DYDKXZ*LIEK_BY +
-     &DZDK_DZDKXZ*LIEK_BZ + DZDB_DZDBX*LIEK_KXX + DZDB_DZDBY*LIEK_KXY +
-     &(DXDB_DXDBX + DZDB_DZDBZ)*LIEK_KXZ + DXDB_DXDBY*LIEK_KYZ +
-     &DXDB_DXDBZ*LIEK_KZZ
+            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 = DXDK_DXDKYY*LIEK_BX + DYDK_DYDKYY*LIEK_BY +
-     &DZDK_DZDKYY*LIEK_BZ + 2*(DYDB_DYDBX*LIEK_KXY + DYDB_DYDBY*LIEK_KYY +
-     &DYDB_DYDBZ*LIEK_KYZ)
+            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_LKYZ = DXDK_DXDKYZ*LIEK_BX + DYDK_DYDKYZ*LIEK_BY +
-     &DZDK_DZDKYZ*LIEK_BZ + DZDB_DZDBX*LIEK_KXY + DYDB_DYDBX*LIEK_KXZ +
-     &DZDB_DZDBY*LIEK_KYY + (DYDB_DYDBY + DZDB_DZDBZ)*LIEK_KYZ +
-     &DYDB_DYDBZ*LIEK_KZZ
+            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))
+
+         else
+
+            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))
+
+         end if
+
+         if (LIEK_BY.gt.0.0D0) then
+
+            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))
+
+         else
+
+            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))
+
+         end if
+
+         if (LIEK_BZ.gt.0.0D0) then
+
+            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))
+
+         else
+
+            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))
+
+         end if
+
+      end if
+
+c     Extra terms in the Lie derivative.
+
+      LIEK_LKXX = LIEK_LKXX + 2.0D0*(DXDB_DXDBX*LIEK_KXX
+     &   + DXDB_DXDBY*LIEK_KXY + DXDB_DXDBZ*LIEK_KXZ)
+
+      LIEK_LKYY = LIEK_LKYY + 2.0D0*(DYDB_DYDBX*LIEK_KXY
+     &   + DYDB_DYDBY*LIEK_KYY + DYDB_DYDBZ*LIEK_KYZ)
+
+      LIEK_LKZZ = LIEK_LKZZ + 2.0D0*(DZDB_DZDBX*LIEK_KXZ
+     &   + DZDB_DZDBY*LIEK_KYZ + DZDB_DZDBZ*LIEK_KZZ)
+
+      LIEK_LKXY = LIEK_LKXY + DYDB_DYDBX*LIEK_KXX + DXDB_DXDBY*LIEK_KYY
+     &   + (DXDB_DXDBX + DYDB_DYDBY)*LIEK_KXY
+     &   + DYDB_DYDBZ*LIEK_KXZ + DXDB_DXDBZ*LIEK_KYZ
+
+      LIEK_LKXZ = LIEK_LKXZ + DZDB_DZDBX*LIEK_KXX + DXDB_DXDBZ*LIEK_KZZ
+     &   + (DXDB_DXDBX + DZDB_DZDBZ)*LIEK_KXZ
+     &   + DZDB_DZDBY*LIEK_KXY + DXDB_DXDBY*LIEK_KYZ
+
+      LIEK_LKYZ = LIEK_LKYZ + DZDB_DZDBY*LIEK_KYY + DYDB_DYDBZ*LIEK_KZZ
+     &   + (DYDB_DYDBY + DZDB_DZDBZ)*LIEK_KYZ
+     &   + DZDB_DZDBX*LIEK_KXY + DYDB_DYDBX*LIEK_KXZ
 
-      LIEK_LKZZ = DXDK_DXDKZZ*LIEK_BX + DYDK_DYDKZZ*LIEK_BY +
-     &DZDK_DZDKZZ*LIEK_BZ + 2*(DZDB_DZDBX*LIEK_KXZ + DZDB_DZDBY*LIEK_KYZ +
-     &DZDB_DZDBZ*LIEK_KZZ)
 
 #endif
 
 #ifdef CCODE
 
-      LIEK_LKXX = DXDK_DXDKXX*LIEK_BX + DYDK_DYDKXX*LIEK_BY +
-      DZDK_DZDKXX*LIEK_BZ + 2*(DXDB_DXDBX*LIEK_KXX + DXDB_DXDBY*LIEK_KXY +
-      DXDB_DXDBZ*LIEK_KXZ) ;
+/* Advection */
+
+      LIEK_LKXX = DXDK_DXDKXX*LIEK_BX + DYDK_DYDKXX*LIEK_BY
+         + DZDK_DZDKXX*LIEK_BZ
+
+      LIEK_LKYY = DXDK_DXDKYY*LIEK_BX + DYDK_DYDKYY*LIEK_BY
+         + DZDK_DZDKYY*LIEK_BZ
+
+      LIEK_LKZZ = DXDK_DXDKZZ*LIEK_BX + DYDK_DYDKZZ*LIEK_BY
+         + DZDK_DZDKZZ*LIEK_BZ
+
+      LIEK_LKXY = DXDK_DXDKXY*LIEK_BX + DYDK_DYDKXY*LIEK_BY
+         + DZDK_DZDKXY*LIEK_BZ
+
+      LIEK_LKXZ = DXDK_DXDKXZ*LIEK_BX + DYDK_DYDKXZ*LIEK_BY
+         + DZDK_DZDKXZ*LIEK_BZ
+
+      LIEK_LKYZ = DXDK_DXDKYZ*LIEK_BX + DYDK_DYDKYZ*LIEK_BY
+         + DZDK_DZDKYZ*LIEK_BZ
+
+/* Extra terms in the Lie derivative */
+
+      LIEK_LKXX = LIEK_LKXX + 2*(DXDB_DXDBX*LIEK_KXX
+         + DXDB_DXDBY*LIEK_KXY + DXDB_DXDBZ*LIEK_KXZ);
 
-      LIEK_LKXY = DXDK_DXDKXY*LIEK_BX + DYDK_DYDKXY*LIEK_BY +
-      DZDK_DZDKXY*LIEK_BZ + DYDB_DYDBX*LIEK_KXX + (DXDB_DXDBX +
-      DYDB_DYDBY)*LIEK_KXY + DYDB_DYDBZ*LIEK_KXZ + DXDB_DXDBY*LIEK_KYY +
-      DXDB_DXDBZ*LIEK_KYZ ;
+      LIEK_LKYY = LIEK_LKYY + 2*(DYDB_DYDBX*LIEK_KXY
+         + DYDB_DYDBY*LIEK_KYY + DYDB_DYDBZ*LIEK_KYZ);
 
-      LIEK_LKXZ = DXDK_DXDKXZ*LIEK_BX + DYDK_DYDKXZ*LIEK_BY +
-      DZDK_DZDKXZ*LIEK_BZ + DZDB_DZDBX*LIEK_KXX + DZDB_DZDBY*LIEK_KXY +
-      (DXDB_DXDBX + DZDB_DZDBZ)*LIEK_KXZ + DXDB_DXDBY*LIEK_KYZ +
-      DXDB_DXDBZ*LIEK_KZZ ;
+      LIEK_LKZZ = LIEK_LKZZ + 2*(DZDB_DZDBX*LIEK_KXZ
+         + DZDB_DZDBY*LIEK_KYZ + DZDB_DZDBZ*LIEK_KZZ);
 
-      LIEK_LKYY = DXDK_DXDKYY*LIEK_BX + DYDK_DYDKYY*LIEK_BY +
-      DZDK_DZDKYY*LIEK_BZ + 2*(DYDB_DYDBX*LIEK_KXY + DYDB_DYDBY*LIEK_KYY +
-      DYDB_DYDBZ*LIEK_KYZ) ;
+      LIEK_LKXY = LIEK_LKXY + DYDB_DYDBX*LIEK_KXX + DXDB_DXDBY*LIEK_KYY
+         + (DXDB_DXDBX + DYDB_DYDBY)*LIEK_KXY
+         + DYDB_DYDBZ*LIEK_KXZ + DXDB_DXDBZ*LIEK_KYZ;
 
-      LIEK_LKYZ = DXDK_DXDKYZ*LIEK_BX + DYDK_DYDKYZ*LIEK_BY +
-      DZDK_DZDKYZ*LIEK_BZ + DZDB_DZDBX*LIEK_KXY + DYDB_DYDBX*LIEK_KXZ +
-      DZDB_DZDBY*LIEK_KYY + (DYDB_DYDBY + DZDB_DZDBZ)*LIEK_KYZ +
-      DYDB_DYDBZ*LIEK_KZZ ;
+      LIEK_LKXZ = LIEK_LKXZ + DZDB_DZDBX*LIEK_KXX + DXDB_DXDBZ*LIEK_KZZ
+         + (DXDB_DXDBX + DZDB_DZDBZ)*LIEK_KXZ
+         + DZDB_DZDBY*LIEK_KXY + DXDB_DXDBY*LIEK_KYZ;
 
-      LIEK_LKZZ = DXDK_DXDKZZ*LIEK_BX + DYDK_DYDKZZ*LIEK_BY +
-      DZDK_DZDKZZ*LIEK_BZ + 2*(DZDB_DZDBX*LIEK_KXZ + DZDB_DZDBY*LIEK_KYZ +
-      DZDB_DZDBZ*LIEK_KZZ) ;
+      LIEK_LKYZ = LIEK_LKYZ + DZDB_DZDBY*LIEK_KYY + DYDB_DYDBZ*LIEK_KZZ
+         + (DYDB_DYDBY + DZDB_DZDBZ)*LIEK_KYZ
+         + DZDB_DZDBX*LIEK_KXY + DYDB_DYDBX*LIEK_KXZ;
 
 #endif