To use this macro you need to set the following variables (which are already defined

in the declare part)

admmacros_advectionx
admmacros_advectiony
admmacros_advectionz

to have the values

0 centered differencing
-1,1 first order upwind for the two different directions
-2,2 second order upwind for the two different directions


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

2 files changed, 155 insertions, 201 deletions

diff --git a/src/macro/LIEG_declare.h b/src/macro/LIEG_declare.h
index b4765d9..4737407 100644
--- a/src/macro/LIEG_declare.h
+++ b/src/macro/LIEG_declare.h
@@ -59,6 +59,10 @@
       CCTK_REAL LIEG_LGYZ
       CCTK_REAL LIEG_LGZZ
 
+      integer admmacros_advectionx
+      integer admmacros_advectiony
+      integer admmacros_advectionz
+
 #endif
 
 #ifdef CCODE
diff --git a/src/macro/LIEG_guts.h b/src/macro/LIEG_guts.h
index 6376815..544f067 100644
--- a/src/macro/LIEG_guts.h
+++ b/src/macro/LIEG_guts.h
@@ -30,286 +30,236 @@
       LIEG_LGXZ = 0.0D0
       LIEG_LGYZ = 0.0D0
 
-      if (CCTK_Equals(advection,'center').eq.1) then
+      /* Do the x-direction */
 
-         LIEG_LGXX = DXDCG_DXDCGXX*LIEG_BX + DYDCG_DYDCGXX*LIEG_BY
-     &      + DZDCG_DZDCGXX*LIEG_BZ
+      if (admmacros_advectionx .eq. 0) then /* CENTER */
 
-         LIEG_LGYY = DXDCG_DXDCGYY*LIEG_BX + DYDCG_DYDCGYY*LIEG_BY
-     &      + DZDCG_DZDCGYY*LIEG_BZ
+         LIEG_LGXX = LIEG_LGXX + DXDCG_DXDCGXX*LIEG_BX 
+         LIEG_LGYY = LIEG_LGYY + DXDCG_DXDCGYY*LIEG_BX 
+         LIEG_LGZZ = LIEG_LGZZ + DXDCG_DXDCGZZ*LIEG_BX 
+         LIEG_LGXY = LIEG_LGXY + DXDCG_DXDCGXY*LIEG_BX
+         LIEG_LGXZ = LIEG_LGXZ + DXDCG_DXDCGXZ*LIEG_BX 
+         LIEG_LGYZ = LIEG_LGYZ + DXDCG_DXDCGYZ*LIEG_BX 
 
-         LIEG_LGZZ = DXDCG_DXDCGZZ*LIEG_BX + DYDCG_DYDCGZZ*LIEG_BY
-     &      + DZDCG_DZDCGZZ*LIEG_BZ
+      else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */
 
-         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
+         LIEG_LGXX = LIEG_LGXX + LIEG_BX
      &         *(gxx(i+1,j,k) - gxx(i,j,k))/dx
-
-            LIEG_LGYY = LIEG_LGYY + LIEG_BX
+         LIEG_LGYY = LIEG_LGYY + LIEG_BX
      &         *(gyy(i+1,j,k) - gyy(i,j,k))/dx
-
-            LIEG_LGZZ = LIEG_LGZZ + LIEG_BX
+         LIEG_LGZZ = LIEG_LGZZ + LIEG_BX
      &         *(gzz(i+1,j,k) - gzz(i,j,k))/dx
-
-            LIEG_LGXY = LIEG_LGXY + LIEG_BX
+         LIEG_LGXY = LIEG_LGXY + LIEG_BX
      &         *(gxy(i+1,j,k) - gxy(i,j,k))/dx
-
-            LIEG_LGXZ = LIEG_LGXZ + LIEG_BX
+         LIEG_LGXZ = LIEG_LGXZ + LIEG_BX
      &         *(gxz(i+1,j,k) - gxz(i,j,k))/dx
-
-            LIEG_LGYZ = LIEG_LGYZ + LIEG_BX
+         LIEG_LGYZ = LIEG_LGYZ + LIEG_BX
      &         *(gyz(i+1,j,k) - gyz(i,j,k))/dx
 
-         else
+      else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */
 
-            LIEG_LGXX = LIEG_LGXX + LIEG_BX
+         LIEG_LGXX = LIEG_LGXX + LIEG_BX
      &         *(gxx(i,j,k) - gxx(i-1,j,k))/dx
-
-            LIEG_LGYY = LIEG_LGYY + LIEG_BX
+         LIEG_LGYY = LIEG_LGYY + LIEG_BX
      &         *(gyy(i,j,k) - gyy(i-1,j,k))/dx
-
-            LIEG_LGZZ = LIEG_LGZZ + LIEG_BX
+         LIEG_LGZZ = LIEG_LGZZ + LIEG_BX
      &         *(gzz(i,j,k) - gzz(i-1,j,k))/dx
-
-            LIEG_LGXY = LIEG_LGXY + LIEG_BX
+         LIEG_LGXY = LIEG_LGXY + LIEG_BX
      &         *(gxy(i,j,k) - gxy(i-1,j,k))/dx
-
-            LIEG_LGXZ = LIEG_LGXZ + LIEG_BX
+         LIEG_LGXZ = LIEG_LGXZ + LIEG_BX
      &         *(gxz(i,j,k) - gxz(i-1,j,k))/dx
-
-            LIEG_LGYZ = LIEG_LGYZ + LIEG_BX
+         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
+      else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */
 
-         if (LIEG_BX.gt.0.0D0) then
-
-            LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BX/dx
+         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
+         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
+         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
+         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
+         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
+         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
+      else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */
 
-            LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BX/dx
+         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
+         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
+         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
+         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
+         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
+         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
+      end if
 
-         if (LIEG_BY.gt.0.0D0) then
+        
+      /* Do the y-direction */
 
-            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))
+      if (admmacros_advectiony .eq. 0) then /* CENTER */
 
-            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_LGXX = LIEG_LGXX + DYDCG_DYDCGXX*LIEG_BY
+         LIEG_LGYY = LIEG_LGYY + DYDCG_DYDCGYY*LIEG_BY 
+         LIEG_LGZZ = LIEG_LGZZ + DYDCG_DYDCGZZ*LIEG_BY 
+         LIEG_LGXY = LIEG_LGXY + DYDCG_DYDCGXY*LIEG_BY
+         LIEG_LGXZ = LIEG_LGXZ + DYDCG_DYDCGXZ*LIEG_BY 
+         LIEG_LGYZ = LIEG_LGYZ + DYDCG_DYDCGYZ*LIEG_BY 
 
-            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))
+      else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */
 
-            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_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_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))
+      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
+
+      else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */
 
-            LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BY/dy
+         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
+      else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */
 
-            LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BY/dy
+         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
+         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
+         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
+         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
+         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
+         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
+      end if
+
+
+      /* Do the z-direction */
 
-         if (LIEG_BZ.gt.0.0D0) then
+      if (admmacros_advectionz .eq. 0) then /* CENTER */
 
-            LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BZ/dz
+         LIEG_LGXX = LIEG_LGXX + DZDCG_DZDCGXX*LIEG_BZ 
+         LIEG_LGYY = LIEG_LGYY + DZDCG_DZDCGYY*LIEG_BZ
+         LIEG_LGZZ = LIEG_LGZZ + DZDCG_DZDCGZZ*LIEG_BZ 
+         LIEG_LGXY = LIEG_LGXY + DZDCG_DZDCGXY*LIEG_BZ
+         LIEG_LGXZ = LIEG_LGXZ + DZDCG_DZDCGXZ*LIEG_BZ 
+         LIEG_LGYZ = LIEG_LGYZ + DZDCG_DZDCGYZ*LIEG_BZ 
+
+      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
+
+      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
+
+      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
+         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
+         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
+         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
+         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
+         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
+      else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */
 
-            LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BZ/dz
+         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
+         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
+         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
+         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
+         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
+         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
+      end if 
 
 /*  Extra terms in the Lie derivative.  */