Modified the include files to handle the extra index in the level set

function after changing to vector groups to be able to evolve more
than one level set function at a time.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@125 2a26948c-0e4f-0410-aee8-f1d3e353619c

8 files changed, 460 insertions, 455 deletions

diff --git a/src/include/centered_second.h b/src/include/centered_second.h
index e0d552e..fbdba04 100644
--- a/src/include/centered_second.h
+++ b/src/include/centered_second.h
@@ -7,45 +7,47 @@ idx = half / cctk_delta_space(1)
 idy = half / cctk_delta_space(2)
 idz = half / cctk_delta_space(3)
 
-do k = kzl, kzr
-  do j = jyl, jyr
-    do i = ixl, ixr
-      if ( eh_mask(i,j,k) .ge. 0 ) then
-        if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),1) ) ) then
-          dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
-        else if ( btest(eh_mask(i,j,k),0) ) then
-          dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                                four * f(i+1,j,k) - f(i+2,j,k) )
+do l = 1, eh_number_level_sets
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
+        if ( eh_mask(i,j,k,l) .ge. 0 ) then
+          if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+            dfx(i,j,k,l) = idx * ( f(i+1,j,k,l) - f(i-1,j,k,l) )
+          else if ( btest(eh_mask(i,j,k,l),0) ) then
+            dfx(i,j,k,l) = idx * ( -three * f(i,j,k,l) + &
+                                  four * f(i+1,j,k,l) - f(i+2,j,k,l) )
+          else
+            dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - &
+                                 four * f(i-1,j,k,l) + f(i-2,j,k,l) )
+          end if
+          if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+            dfy(i,j,k,l) = idy * ( f(i,j+1,k,l) - f(i,j-1,k,l) )
+          else if ( btest(eh_mask(i,j,k,l),2) ) then
+            dfy(i,j,k,l) = idy * ( -three * f(i,j,k,l) + &
+                                  four * f(i,j+1,k,l) - f(i,j+2,k,l) )
+          else
+            dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - &
+                                  four * f(i,j-1,k,l) + f(i,j-2,k,l) )
+          end if
+          if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+            dfz(i,j,k,l) = idz * ( f(i,j,k+1,l) - f(i,j,k-1,l) )
+          else if ( btest(eh_mask(i,j,k,l),4) ) then
+            dfz(i,j,k,l) = idz * ( -three * f(i,j,k,l) + &
+                                   four * f(i,j,k+1,l) - f(i,j,k+2,l) )
+          else
+            dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - &
+                                  four * f(i,j,k-1,l) + f(i,j,k-2,l) )
+          end if
         else
-          dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                               four * f(i-1,j,k) + f(i-2,j,k) )
+          dfx(i,j,k,l) = zero
+          dfy(i,j,k,l) = zero
+          dfz(i,j,k,l) = zero
         end if
-        if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),3) ) ) then
-          dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
-        else if ( btest(eh_mask(i,j,k),2) ) then
-          dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                                four * f(i,j+1,k) - f(i,j+2,k) )
-        else
-          dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                                four * f(i,j-1,k) + f(i,j-2,k) )
-        end if
-        if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),5) ) ) then
-          dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
-        else if ( btest(eh_mask(i,j,k),4) ) then
-          dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                                 four * f(i,j,k+1) - f(i,j,k+2) )
-        else
-          dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                                four * f(i,j,k-1) + f(i,j,k-2) )
-        end if
-      else
-        dfx(i,j,k) = zero
-        dfy(i,j,k) = zero
-        dfz(i,j,k) = zero
-      end if
+      end do
     end do
   end do
 end do
diff --git a/src/include/centered_second2.h b/src/include/centered_second2.h
index c29fd41..e7aca4b 100644
--- a/src/include/centered_second2.h
+++ b/src/include/centered_second2.h
@@ -1,39 +1,39 @@
 ! Calculation of centered differences.
 ! $Header$
 
-if ( eh_mask(i,j,k) .ge. 0 ) then
-  if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),1) ) ) then
-    dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
-  else if ( btest(eh_mask(i,j,k),0) ) then
-    dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                          four * f(i+1,j,k) - f(i+2,j,k) )
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+    dfx(i,j,k,l) = idx * ( f(i+1,j,k,l) - f(i-1,j,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),0) ) then
+    dfx(i,j,k,l) = idx * ( -three * f(i,j,k,l) + &
+                          four * f(i+1,j,k,l) - f(i+2,j,k,l) )
   else
-    dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                         four * f(i-1,j,k) + f(i-2,j,k) )
+    dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - &
+                         four * f(i-1,j,k,l) + f(i-2,j,k,l) )
   end if
-  if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),3) ) ) then
-    dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
-  else if ( btest(eh_mask(i,j,k),2) ) then
-    dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                          four * f(i,j+1,k) - f(i,j+2,k) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+    dfy(i,j,k,l) = idy * ( f(i,j+1,k,l) - f(i,j-1,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),2) ) then
+    dfy(i,j,k,l) = idy * ( -three * f(i,j,k,l) + &
+                          four * f(i,j+1,k,l) - f(i,j+2,k,l) )
   else
-    dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                          four * f(i,j-1,k) + f(i,j-2,k) )
+    dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - &
+                          four * f(i,j-1,k,l) + f(i,j-2,k,l) )
   end if
-  if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),5) ) ) then
-    dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
-  else if ( btest(eh_mask(i,j,k),4) ) then
-    dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                           four * f(i,j,k+1) - f(i,j,k+2) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+    dfz(i,j,k,l) = idz * ( f(i,j,k+1,l) - f(i,j,k-1,l) )
+  else if ( btest(eh_mask(i,j,k,l),4) ) then
+    dfz(i,j,k,l) = idz * ( -three * f(i,j,k,l) + &
+                           four * f(i,j,k+1,l) - f(i,j,k+2,l) )
   else
-    dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                          four * f(i,j,k-1) + f(i,j,k-2) )
+    dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - &
+                          four * f(i,j,k-1,l) + f(i,j,k-2,l) )
   end if
 else
-  dfx(i,j,k) = zero
-  dfy(i,j,k) = zero
-  dfz(i,j,k) = zero
+  dfx(i,j,k,l) = zero
+  dfy(i,j,k,l) = zero
+  dfz(i,j,k,l) = zero
 end if
diff --git a/src/include/metric.h b/src/include/metric.h
index f3a22ff..d4890df 100644
--- a/src/include/metric.h
+++ b/src/include/metric.h
@@ -1,7 +1,7 @@
 ! Calculation of inverse three metric, lapse squared and shift
 ! $Header$
 
-if ( eh_mask(i,j,k) .ge. 0 ) then
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
   alp2 = alp(i,j,k)**2
 
   tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
diff --git a/src/include/upwind_characteristic_second2.h b/src/include/upwind_characteristic_second2.h
index 5dacaa4..04b1d06 100644
--- a/src/include/upwind_characteristic_second2.h
+++ b/src/include/upwind_characteristic_second2.h
@@ -2,51 +2,51 @@
 ! $Header$
 
 ! If this is an active cell we have to compute its derivative.
-if ( eh_mask(i,j,k) .ge. 0 ) then
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
 
   ! If it is not a boundary cell in the x-direction then calculate as a
   ! starting point the centered derivative in the x-direction.
-  if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),1) ) ) then
-    dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+    dfx(i,j,k,l) = idx * ( f(i+1,j,k,l) - f(i-1,j,k,l) )
 
 !   Otherwise we will have to calculate the appropriate one sided derivative.
-  else if ( btest(eh_mask(i,j,k),0) ) then
-    dfx(i,j,k) = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+  else if ( btest(eh_mask(i,j,k,l),0) ) then
+    dfx(i,j,k,l) = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
   else
-    dfx(i,j,k) = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
+    dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
   end if
 
   ! Do the same for the y direction.
-  if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),3) ) ) then
-    dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
-  else if ( btest(eh_mask(i,j,k),2) ) then
-    dfy(i,j,k) = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+    dfy(i,j,k,l) = idy * ( f(i,j+1,k,l) - f(i,j-1,k,l) )
+  else if ( btest(eh_mask(i,j,k,l),2) ) then
+    dfy(i,j,k,l) = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
   else
-    dfy(i,j,k) = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
+    dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
   end if
 
   ! And for the z direction.
-  if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),5) ) ) then
-    dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
-  else if ( btest(eh_mask(i,j,k),4) ) then
-    dfz(i,j,k) = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+    dfz(i,j,k,l) = idz * ( f(i,j,k+1,l) - f(i,j,k-1,l) )
+  else if ( btest(eh_mask(i,j,k,l),4) ) then
+    dfz(i,j,k,l) = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
   else
-    dfz(i,j,k) = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
+    dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
   end if
 
 ! Calculate the characteristic direction using these preliminary
 ! derivatives and multiply it with the timestep.
 
-  dfup(1) = g3xx * dfx(i,j,k) + g3xy * dfy(i,j,k) + g3xz * dfz(i,j,k)
-  dfup(2) = g3xy * dfx(i,j,k) + g3yy * dfy(i,j,k) + g3yz * dfz(i,j,k)
-  dfup(3) = g3xz * dfx(i,j,k) + g3yz * dfy(i,j,k) + g3zz * dfz(i,j,k)
+  dfup(1) = g3xx * dfx(i,j,k,l) + g3xy * dfy(i,j,k,l) + g3xz * dfz(i,j,k,l)
+  dfup(2) = g3xy * dfx(i,j,k,l) + g3yy * dfy(i,j,k,l) + g3yz * dfz(i,j,k,l)
+  dfup(3) = g3xz * dfx(i,j,k,l) + g3yz * dfy(i,j,k,l) + g3zz * dfz(i,j,k,l)
 
-  cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k) + &
-                                  dfup(2) * dfy(i,j,k) + &
-                                  dfup(3) * dfz(i,j,k) ) )
+  cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k,l) + &
+                                  dfup(2) * dfy(i,j,k,l) + &
+                                  dfup(3) * dfz(i,j,k,l) ) )
 
   cdx(1) = ( -betax(i,j,k) + ssign * alp2 * dfup(1) * cfactor ) * &
                                                       cctk_delta_time
@@ -55,8 +55,8 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
   cdx(3) = ( -betaz(i,j,k) + ssign * alp2 * dfup(3) * cfactor ) * &
                                                       cctk_delta_time
 
-  if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),1) ) ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
     ! Look at the characteristic direction to figure out in which direction
     ! the stencil should be chosen.
     if ( cdx(1) .gt. zero ) then           ! Choose the stencil to the left.
@@ -65,15 +65,15 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
       ! safe to assume that i-2, i-1,i and i+1 have good values.
       ! Otherwise the neigbour must be a boundary cell, so we can only use a
       ! centered derivative, which we have already computed.
-      if ( eh_mask(i-1,j,k) .eq. 0 ) then
+      if ( eh_mask(i-1,j,k,l) .eq. 0 ) then
 
         ! Choose the stencil, that gives the smalles second derivative.
         ! If it happens to be the centered, do not do anything since we have
         ! already calculated the centered derivative.
-!        if ( f(i-2,j,k) - two * f(i-1,j,k) + f(i,j,k) .le. &
-!             f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) then
-          dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                               four * f(i-1,j,k) + f(i-2,j,k) )
+!        if ( f(i-2,j,k,l) - two * f(i-1,j,k,l) + f(i,j,k,l) .le. &
+!             f(i-1,j,k,l) - two * f(i,j,k,l) + f(i+1,j,k,l) ) then
+          dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - &
+                               four * f(i-1,j,k,l) + f(i-2,j,k,l) )
 !        end if
       end if
 
@@ -83,67 +83,67 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
       ! safe to assume that i-1, i,i+1 and i+2 have good values.
       ! Otherwise the neigbour must be a boundary cell, so we can only use a
       ! centered derivative, which we have already computed.
-      if ( eh_mask(i+1,j,k) .eq. 0 ) then
+      if ( eh_mask(i+1,j,k,l) .eq. 0 ) then
 
         ! Choose the stencil, that gives the smalles second derivative.
         ! If it happens to be the centered, do not do anything since we have
         ! already calculated the centered derivative.
-!        if ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) .ge. &
-!             f(i,j,k) - two * f(i+1,j,k) + f(i+2,j,k) ) then
-          dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                               four * f(i+1,j,k) - f(i+2,j,k) )
+!        if ( f(i-1,j,k,l) - two * f(i,j,k,l) + f(i+1,j,k,l) .ge. &
+!             f(i,j,k,l) - two * f(i+1,j,k,l) + f(i+2,j,k,l) ) then
+          dfx(i,j,k,l) = idx * ( -three * f(i,j,k,l) + &
+                               four * f(i+1,j,k,l) - f(i+2,j,k,l) )
 !        end if
       end if
     end if
   end if
   
   ! Do the same for the y direction.
-  if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),3) ) ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
 
     if ( cdx(2) .gt. zero ) then           ! Choose the stencil to the left.
-      if ( eh_mask(i,j-1,k) .eq. 0 ) then
-!        if ( f(i,j-2,k) - two * f(i,j-1,k) + f(i,j,k) .le. &
-!             f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) then
-          dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                               four * f(i,j-1,k) + f(i,j-2,k) )
+      if ( eh_mask(i,j-1,k,l) .eq. 0 ) then
+!        if ( f(i,j-2,k,l) - two * f(i,j-1,k,l) + f(i,j,k,l) .le. &
+!             f(i,j-1,k,l) - two * f(i,j,k,l) + f(i,j+1,k,l) ) then
+          dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - &
+                               four * f(i,j-1,k,l) + f(i,j-2,k,l) )
 !        end if
       end if
     else if ( cdx(2) .lt. zero ) then      ! Choose the stencil to the left.
-      if ( eh_mask(i,j+1,k) .eq. 0 ) then
-!        if ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) .ge. &
-!             f(i,j,k) - two * f(i,j+1,k) + f(i,j+2,k) ) then
-          dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                               four * f(i,j+1,k) - f(i,j+2,k) )
+      if ( eh_mask(i,j+1,k,l) .eq. 0 ) then
+!        if ( f(i,j-1,k,l) - two * f(i,j,k,l) + f(i,j+1,k,l) .ge. &
+!             f(i,j,k,l) - two * f(i,j+1,k,l) + f(i,j+2,k,l) ) then
+          dfy(i,j,k,l) = idy * ( -three * f(i,j,k,l) + &
+                               four * f(i,j+1,k,l) - f(i,j+2,k,l) )
 !        end if
       end if
     end if
   end if
 
   ! And for the z direction.
-  if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),5) ) ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
 
     if ( cdx(3) .gt. zero ) then
-      if ( eh_mask(i,j,k-1) .eq. 0 ) then
-!        if ( f(i,j,k-2) - two * f(i,j,k-1) + f(i,j,k) .le. &
-!             f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) then
-          dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                               four * f(i,j,k-1) + f(i,j,k-2) )
+      if ( eh_mask(i,j,k-1,l) .eq. 0 ) then
+!        if ( f(i,j,k-2,l) - two * f(i,j,k-1,l) + f(i,j,k,l) .le. &
+!             f(i,j,k-1,l) - two * f(i,j,k,l) + f(i,j,k+1,l) ) then
+          dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - &
+                               four * f(i,j,k-1,l) + f(i,j,k-2,l) )
 !        end if
       end if
     else if ( cdx(3) .lt. zero ) then
-      if ( eh_mask(i,j,k+1) .eq. 0 ) then
-!        if ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) .ge. &
-!             f(i,j,k) - two * f(i,j,k+1) + f(i,j,k+2) ) then
-          dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                               four * f(i,j,k+1) - f(i,j,k+2) )
+      if ( eh_mask(i,j,k+1,l) .eq. 0 ) then
+!        if ( f(i,j,k-1,l) - two * f(i,j,k,l) + f(i,j,k+1,l) .ge. &
+!             f(i,j,k,l) - two * f(i,j,k+1,l) + f(i,j,k+2,l) ) then
+          dfz(i,j,k,l) = idz * ( -three * f(i,j,k,l) + &
+                               four * f(i,j,k+1,l) - f(i,j,k+2,l) )
 !        end if
       end if
     end if
   end if
 else
-  dfx(i,j,k) = zero
-  dfy(i,j,k) = zero
-  dfz(i,j,k) = zero
+  dfx(i,j,k,l) = zero
+  dfy(i,j,k,l) = zero
+  dfz(i,j,k,l) = zero
 end if
diff --git a/src/include/upwind_first.h b/src/include/upwind_first.h
index 6c99d29..1f8c309 100644
--- a/src/include/upwind_first.h
+++ b/src/include/upwind_first.h
@@ -5,100 +5,101 @@ idx = one / cctk_delta_space(1)
 idy = one / cctk_delta_space(2)
 idz = one / cctk_delta_space(3)
 
-do k = 1, nz
-  do j = 1, ny
-    do i = 1, nx
-      if ( eh_mask(i,j,k) .ge. 0 ) then
-        if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),1) ) ) then
-          al = idx * ( f(i,j,k) - f(i-1,j,k) )
-          ar = idx * ( f(i+1,j,k) - f(i,j,k) )
-        else if ( btest(eh_mask(i,j,k),0) ) then
-          al = zero
-          ar = idx * ( f(i+1,j,k) - f(i,j,k) )
-        else
-          al = idx * ( f(i,j,k) - f(i-1,j,k) )
-          ar = zero
-        end if
-        if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),3) ) ) then
-          bl = idy * ( f(i,j,k) - f(i,j-1,k) )
-          br = idy * ( f(i,j+1,k) - f(i,j,k) )
-        else if ( btest(eh_mask(i,j,k),2) ) then
-          bl = zero
-          br = idy * ( f(i,j+1,k) - f(i,j,k) )
-        else
-          bl = idy * ( f(i,j,k) - f(i,j-1,k) )
-          br = zero
-        end if
-        if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),5) ) ) then
-          cl = idz * ( f(i,j,k) - f(i,j,k-1) )
-          cr = idz * ( f(i,j,k+1) - f(i,j,k) )
-        else if ( btest(eh_mask(i,j,k),4) ) then
-          cl = zero
-          cr = idz * ( f(i,j,k+1) - f(i,j,k) )
-        else
-          cl = idz * ( f(i,j,k) - f(i,j,k-1) )
-          cr = zero
-        end if
-
-        alminus = half*(al-abs(al))
-        alplus = half*(al+abs(al))
-        arminus = half*(ar-abs(ar))
-        arplus = half*(ar+abs(ar))
-
-        blminus = half*(bl-abs(bl))
-        blplus = half*(bl+abs(bl))
-        brminus = half*(br-abs(br))
-        brplus = half*(br+abs(br))
-
-        clminus = half*(cl-abs(cl))
-        clplus = half*(cl+abs(cl))
-        crminus = half*(cr-abs(cr))
-        crplus = half*(cr+abs(cr))
-
-        if ( f(i,j,k) .gt. 0 ) then
-          if ( abs(alplus) .gt. abs(arminus) ) then
-            dfx(i,j,k) = alplus
+do l = 1, eh_number_level_sets
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        if ( eh_mask(i,j,k,l) .ge. 0 ) then
+          if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+            al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
+            ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
+          else if ( btest(eh_mask(i,j,k,l),0) ) then
+            al = zero
+            ar = idx * ( f(i+1,j,k,l) - f(i,j,k,l) )
           else
-            dfx(i,j,k) = arminus
+            al = idx * ( f(i,j,k,l) - f(i-1,j,k,l) )
+            ar = zero
           end if
-          if ( abs(blplus) .gt. abs(brminus) ) then
-            dfy(i,j,k) = blplus
+          if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+            bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
+            br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
+          else if ( btest(eh_mask(i,j,k,l),2) ) then
+            bl = zero
+            br = idy * ( f(i,j+1,k,l) - f(i,j,k,l) )
           else
-            dfy(i,j,k) = brminus
+            bl = idy * ( f(i,j,k,l) - f(i,j-1,k,l) )
+            br = zero
           end if
-          if ( abs(clplus) .gt. abs(crminus) ) then
-            dfz(i,j,k) = clplus
+          if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+            cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
+            cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
+          else if ( btest(eh_mask(i,j,k,l),4) ) then
+            cl = zero
+            cr = idz * ( f(i,j,k+1,l) - f(i,j,k,l) )
           else
-            dfz(i,j,k) = crminus
+            cl = idz * ( f(i,j,k,l) - f(i,j,k-1,l) )
+            cr = zero
           end if
-        endif
-
-        if ( f(i,j,k) .lt. 0 ) then
-          if ( abs(alminus) .gt. abs(arplus) ) then
-            dfx(i,j,k) = alminus
-          else
-            dfx(i,j,k) = arplus
-          end if
-          if ( abs(blminus) .gt. abs(brplus) ) then
-            dfy(i,j,k) = blminus
-          else
-            dfy(i,j,k) = brplus
-          end if
-          if ( abs(clminus) .gt. abs(crplus) ) then
-            dfz(i,j,k) = clminus
-          else
-            dfz(i,j,k) = crplus
+  
+          alminus = half*(al-abs(al))
+          alplus = half*(al+abs(al))
+          arminus = half*(ar-abs(ar))
+          arplus = half*(ar+abs(ar))
+  
+          blminus = half*(bl-abs(bl))
+          blplus = half*(bl+abs(bl))
+          brminus = half*(br-abs(br))
+          brplus = half*(br+abs(br))
+  
+          clminus = half*(cl-abs(cl))
+          clplus = half*(cl+abs(cl))
+          crminus = half*(cr-abs(cr))
+          crplus = half*(cr+abs(cr))
+  
+          if ( f(i,j,k,l) .gt. 0 ) then
+            if ( abs(alplus) .gt. abs(arminus) ) then
+              dfx(i,j,k,l) = alplus
+            else
+              dfx(i,j,k,l) = arminus
+            end if
+            if ( abs(blplus) .gt. abs(brminus) ) then
+              dfy(i,j,k,l) = blplus
+            else
+              dfy(i,j,k,l) = brminus
+            end if
+            if ( abs(clplus) .gt. abs(crminus) ) then
+              dfz(i,j,k,l) = clplus
+            else
+              dfz(i,j,k,l) = crminus
+            end if
+          endif
+  
+          if ( f(i,j,k,l) .lt. 0 ) then
+            if ( abs(alminus) .gt. abs(arplus) ) then
+              dfx(i,j,k,l) = alminus
+            else
+              dfx(i,j,k,l) = arplus
+            end if
+            if ( abs(blminus) .gt. abs(brplus) ) then
+              dfy(i,j,k,l) = blminus
+            else
+              dfy(i,j,k,l) = brplus
+            end if
+            if ( abs(clminus) .gt. abs(crplus) ) then
+              dfz(i,j,k,l) = clminus
+            else
+              dfz(i,j,k,l) = crplus
+            end if
           end if
+        else
+          dfx(i,j,k,l) = zero
+          dfy(i,j,k,l) = zero
+          dfz(i,j,k,l) = zero
         end if
-      else
-        dfx(i,j,k) = zero
-        dfy(i,j,k) = zero
-        dfz(i,j,k) = zero
-      end if
+      end do
     end do
   end do
 end do
-
diff --git a/src/include/upwind_second.h b/src/include/upwind_second.h
index f6004e9..fd899c5 100644
--- a/src/include/upwind_second.h
+++ b/src/include/upwind_second.h
@@ -15,162 +15,164 @@ idz = half / cctk_delta_space(3)
 ! Calculate appropriate one sided derivatives for all active cells.
 ! Cells are active if eh_mask is greater or equal to zero. They are interior
 ! cells if eh_mask == 0 and boundary cells if eh_mask > 0.
-do k = kzl, kzr
-  do j = jyl, jyr
-    do i = ixl, ixr
-
-      ! If this is an active cell...
-      if ( eh_mask(i,j,k) .ge. 0 ) then
-
-        ! If this is not a boundary cell...
-        if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),1) ) ) then
-
-          ! If the cell to the left is not a boundary cell...
-          if ( ( .not. btest(eh_mask(i-1,j,k),0) ) .or. &
-                     ( eh_mask(i-1,j,k) .eq. -2 ) ) then
-
-            ! Calculate left handed one sided derivative
-            al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
-          else
+do l = 1, eh_number_level_sets
+  do k = kzl, kzr
+    do j = jyl, jyr
+      do i = ixl, ixr
+  
+        ! If this is an active cell...
+        if ( eh_mask(i,j,k,l) .ge. 0 ) then
+  
+          ! If this is not a boundary cell...
+          if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
+  
+            ! If the cell to the left is not a boundary cell...
+            if ( ( .not. btest(eh_mask(i-1,j,k,l),0) ) .or. &
+                       ( eh_mask(i-1,j,k,l) .eq. -2 ) ) then
+  
+              ! Calculate left handed one sided derivative
+              al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
+            else
+              al = zero
+            end if
+  
+            ! If the cell to the right is not a boundary cell...
+            if ( ( .not. btest(eh_mask(i+1,j,k,l),1) ) .or. &
+                       ( eh_mask(i+1,j,k,l) .eq. -2 ) ) then 
+  
+              ! Calculate right handed one sided derivative
+              ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
+            else
+              ar = zero
+            end if
+  
+          ! Else if the cell is a left boundary cell.
+          else if ( btest(eh_mask(i,j,k,l),0) ) then
+  
+            ! calculate only right handed one sided difference.
             al = zero
-          end if
-
-          ! If the cell to the right is not a boundary cell...
-          if ( ( .not. btest(eh_mask(i+1,j,k),1) ) .or. &
-                     ( eh_mask(i+1,j,k) .eq. -2 ) ) then 
-
-            ! Calculate right handed one sided derivative
-            ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+            ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
+  
+          ! Else it must be a right boundary cell.
           else
+  
+            ! So calculate only the left handed one sided difference.
+            al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
             ar = zero
           end if
-
-        ! Else if the cell is a left boundary cell.
-        else if ( btest(eh_mask(i,j,k),0) ) then
-
-          ! calculate only right handed one sided difference.
-          al = zero
-          ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
-
-        ! Else it must be a right boundary cell.
-        else
-
-          ! So calculate only the left handed one sided difference.
-          al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
-          ar = zero
-        end if
-
-        ! Do the same for the y-derivative.
-        if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),3) ) ) then
-          if ( ( .not. btest(eh_mask(i,j-1,k),2) ) .or. &
-                     ( eh_mask(i,j-1,k) .eq. -2 ) ) then
-            bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
-          else
+  
+          ! Do the same for the y-derivative.
+          if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+            if ( ( .not. btest(eh_mask(i,j-1,k,l),2) ) .or. &
+                       ( eh_mask(i,j-1,k,l) .eq. -2 ) ) then
+              bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
+            else
+              bl = zero
+            end if
+            if ( ( .not. btest(eh_mask(i,j+1,k,l),3) ) .or. &
+                       ( eh_mask(i,j+1,k,l) .eq. -2 ) ) then 
+              br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
+            else
+              br = zero
+            end if
+          else if ( btest(eh_mask(i,j,k,l),2) ) then
             bl = zero
-          end if
-          if ( ( .not. btest(eh_mask(i,j+1,k),3) ) .or. &
-                     ( eh_mask(i,j+1,k) .eq. -2 ) ) then 
-            br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
+            br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
           else
+            bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
             br = zero
           end if
-        else if ( btest(eh_mask(i,j,k),2) ) then
-          bl = zero
-          br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
-        else
-          bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
-          br = zero
-        end if
-
-        ! And the z-derivative
-        if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-             ( .not. btest(eh_mask(i,j,k),5) ) ) then
-          if ( ( .not. btest(eh_mask(i,j,k-1),4) ) .or. &
-                     ( eh_mask(i,j,k-1) .eq. -2 ) ) then
-            cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
-          else
+  
+          ! And the z-derivative
+          if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+               ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+            if ( ( .not. btest(eh_mask(i,j,k-1,l),4) ) .or. &
+                       ( eh_mask(i,j,k-1,l) .eq. -2 ) ) then
+              cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
+            else
+              cl = zero
+            end if
+            if ( ( .not. btest(eh_mask(i,j,k+1,l),5) ) .or. &
+                       ( eh_mask(i,j,k+1,l) .eq. -2 ) ) then 
+              cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
+            else
+              cr = zero
+            end if
+          else if ( btest(eh_mask(i,j,k,l),4) ) then
             cl = zero
-          end if
-          if ( ( .not. btest(eh_mask(i,j,k+1),5) ) .or. &
-                     ( eh_mask(i,j,k+1) .eq. -2 ) ) then 
-            cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
+            cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
           else
+            cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
             cr = zero
           end if
-        else if ( btest(eh_mask(i,j,k),4) ) then
-          cl = zero
-          cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
-        else
-          cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
-          cr = zero
-        end if
-
-        ! Get the negative and positive parts of the one sided derivatives in
-        ! the x-direction.
-        alminus = half*(al-abs(al))
-        alplus = half*(al+abs(al))
-        arminus = half*(ar-abs(ar))
-        arplus = half*(ar+abs(ar))
-
-        ! Ditto for the y-direction.
-        blminus = half*(bl-abs(bl))
-        blplus = half*(bl+abs(bl))
-        brminus = half*(br-abs(br))
-        brplus = half*(br+abs(br))
-
-        ! Ditto for the z-direction.
-        clminus = half*(cl-abs(cl))
-        clplus = half*(cl+abs(cl))
-        crminus = half*(cr-abs(cr))
-        crplus = half*(cr+abs(cr))
-
-        ! If f>0 choose the correct one sided derivative depending on the
-        ! magnitude of the negative and positive parts
-        if ( f(i,j,k) .gt. 0 ) then
-          if ( abs(alplus) .gt. abs(arminus) ) then
-            dfx(i,j,k) = alplus
-          else
-            dfx(i,j,k) = arminus
-          end if
-          if ( abs(blplus) .gt. abs(brminus) ) then
-            dfy(i,j,k) = blplus
-          else
-            dfy(i,j,k) = brminus
-          end if
-          if ( abs(clplus) .gt. abs(crminus) ) then
-            dfz(i,j,k) = clplus
-          else
-            dfz(i,j,k) = crminus
-          end if
-        endif
-
-        ! Ditto if f<0.
-        if ( f(i,j,k) .lt. 0 ) then
-          if ( abs(alminus) .gt. abs(arplus) ) then
-            dfx(i,j,k) = alminus
-          else
-            dfx(i,j,k) = arplus
-          end if
-          if ( abs(blminus) .gt. abs(brplus) ) then
-            dfy(i,j,k) = blminus
-          else
-            dfy(i,j,k) = brplus
-          end if
-          if ( abs(clminus) .gt. abs(crplus) ) then
-            dfz(i,j,k) = clminus
-          else
-            dfz(i,j,k) = crplus
+  
+          ! Get the negative and positive parts of the one sided derivatives in
+          ! the x-direction.
+          alminus = half*(al-abs(al))
+          alplus = half*(al+abs(al))
+          arminus = half*(ar-abs(ar))
+          arplus = half*(ar+abs(ar))
+  
+          ! Ditto for the y-direction.
+          blminus = half*(bl-abs(bl))
+          blplus = half*(bl+abs(bl))
+          brminus = half*(br-abs(br))
+          brplus = half*(br+abs(br))
+  
+          ! Ditto for the z-direction.
+          clminus = half*(cl-abs(cl))
+          clplus = half*(cl+abs(cl))
+          crminus = half*(cr-abs(cr))
+          crplus = half*(cr+abs(cr))
+  
+          ! If f>0 choose the correct one sided derivative depending on the
+          ! magnitude of the negative and positive parts
+          if ( f(i,j,k,l) .gt. 0 ) then
+            if ( abs(alplus) .gt. abs(arminus) ) then
+              dfx(i,j,k,l) = alplus
+            else
+              dfx(i,j,k,l) = arminus
+            end if
+            if ( abs(blplus) .gt. abs(brminus) ) then
+              dfy(i,j,k,l) = blplus
+            else
+              dfy(i,j,k,l) = brminus
+            end if
+            if ( abs(clplus) .gt. abs(crminus) ) then
+              dfz(i,j,k,l) = clplus
+            else
+              dfz(i,j,k,l) = crminus
+            end if
+          endif
+  
+          ! Ditto if f<0.
+          if ( f(i,j,k,l) .lt. 0 ) then
+            if ( abs(alminus) .gt. abs(arplus) ) then
+              dfx(i,j,k,l) = alminus
+            else
+              dfx(i,j,k,l) = arplus
+            end if
+            if ( abs(blminus) .gt. abs(brplus) ) then
+              dfy(i,j,k,l) = blminus
+            else
+              dfy(i,j,k,l) = brplus
+            end if
+            if ( abs(clminus) .gt. abs(crplus) ) then
+              dfz(i,j,k,l) = clminus
+            else
+              dfz(i,j,k,l) = crplus
+            end if
           end if
+  
+        ! If the cell is not active set the derivatives to zero.
+        else
+          dfx(i,j,k,l) = zero
+          dfy(i,j,k,l) = zero
+          dfz(i,j,k,l) = zero
         end if
-
-      ! If the cell is not active set the derivatives to zero.
-      else
-        dfx(i,j,k) = zero
-        dfy(i,j,k) = zero
-        dfz(i,j,k) = zero
-      end if
+      end do
     end do
   end do
 end do
diff --git a/src/include/upwind_second2.h b/src/include/upwind_second2.h
index d2e3c17..2f7fdb0 100644
--- a/src/include/upwind_second2.h
+++ b/src/include/upwind_second2.h
@@ -4,90 +4,90 @@
 ! $Header$
 
 ! If this is an active cell...
-if ( eh_mask(i,j,k) .ge. 0 ) then
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
 
   ! If this is not a boundary cell...
-  if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),1) ) ) then
+  if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),1) ) ) then
 
     ! If the cell to the left is not a boundary cell...
-    if ( ( .not. btest(eh_mask(i-1,j,k),0) ) .or. &
-               ( eh_mask(i-1,j,k) .eq. -2 ) ) then
+    if ( ( .not. btest(eh_mask(i-1,j,k,l),0) ) .or. &
+               ( eh_mask(i-1,j,k,l) .eq. -2 ) ) then
 
       ! Calculate left handed one sided derivative
-      al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
+      al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
     else
       al = zero
     end if
 
     ! If the cell to the right is not a boundary cell...
-    if ( ( .not. btest(eh_mask(i+1,j,k),1) ) .or. &
-               ( eh_mask(i+1,j,k) .eq. -2 ) ) then 
+    if ( ( .not. btest(eh_mask(i+1,j,k,l),1) ) .or. &
+               ( eh_mask(i+1,j,k,l) .eq. -2 ) ) then 
 
       ! Calculate right handed one sided derivative
-      ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+      ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
     else
       ar = zero
     end if
 
   ! Else if the cell is a left boundary cell.
-  else if ( btest(eh_mask(i,j,k),0) ) then
+  else if ( btest(eh_mask(i,j,k,l),0) ) then
 
     ! calculate only right handed one sided difference.
     al = zero
-    ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
+    ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
 
   ! Else it must be a right boundary cell.
   else
 
     ! So calculate only the left handed one sided difference.
-    al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
+    al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
     ar = zero
   end if
 
   ! Do the same for the y-derivative.
-  if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),3) ) ) then
-    if ( ( .not. btest(eh_mask(i,j-1,k),2) ) .or. &
-               ( eh_mask(i,j-1,k) .eq. -2 ) ) then
-      bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),3) ) ) then
+    if ( ( .not. btest(eh_mask(i,j-1,k,l),2) ) .or. &
+               ( eh_mask(i,j-1,k,l) .eq. -2 ) ) then
+      bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
     else
       bl = zero
     end if
-    if ( ( .not. btest(eh_mask(i,j+1,k),3) ) .or. &
-               ( eh_mask(i,j+1,k) .eq. -2 ) ) then 
-      br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
+    if ( ( .not. btest(eh_mask(i,j+1,k,l),3) ) .or. &
+               ( eh_mask(i,j+1,k,l) .eq. -2 ) ) then 
+      br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
     else
       br = zero
     end if
-  else if ( btest(eh_mask(i,j,k),2) ) then
+  else if ( btest(eh_mask(i,j,k,l),2) ) then
     bl = zero
-    br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
+    br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
   else
-    bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
+    bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
     br = zero
   end if
 
   ! And the z-derivative
-  if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
-       ( .not. btest(eh_mask(i,j,k),5) ) ) then
-    if ( ( .not. btest(eh_mask(i,j,k-1),4) ) .or. &
-               ( eh_mask(i,j,k-1) .eq. -2 ) ) then
-      cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
+  if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
+       ( .not. btest(eh_mask(i,j,k,l),5) ) ) then
+    if ( ( .not. btest(eh_mask(i,j,k-1,l),4) ) .or. &
+               ( eh_mask(i,j,k-1,l) .eq. -2 ) ) then
+      cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
     else
       cl = zero
     end if
-    if ( ( .not. btest(eh_mask(i,j,k+1),5) ) .or. &
-               ( eh_mask(i,j,k+1) .eq. -2 ) ) then 
-      cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
+    if ( ( .not. btest(eh_mask(i,j,k+1,l),5) ) .or. &
+               ( eh_mask(i,j,k+1,l) .eq. -2 ) ) then 
+      cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
     else
       cr = zero
     end if
-  else if ( btest(eh_mask(i,j,k),4) ) then
+  else if ( btest(eh_mask(i,j,k,l),4) ) then
     cl = zero
-    cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
+    cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
   else
-    cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
+    cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
     cr = zero
   end if
 
@@ -112,46 +112,46 @@ if ( eh_mask(i,j,k) .ge. 0 ) then
 
   ! If f>0 choose the correct one sided derivative depending on the
   ! magnitude of the negative and positive parts
-  if ( f(i,j,k) .gt. 0 ) then
+  if ( f(i,j,k,l) .gt. 0 ) then
     if ( abs(alplus) .gt. abs(arminus) ) then
-      dfx(i,j,k) = alplus
+      dfx(i,j,k,l) = alplus
     else
-      dfx(i,j,k) = arminus
+      dfx(i,j,k,l) = arminus
     end if
     if ( abs(blplus) .gt. abs(brminus) ) then
-      dfy(i,j,k) = blplus
+      dfy(i,j,k,l) = blplus
     else
-      dfy(i,j,k) = brminus
+      dfy(i,j,k,l) = brminus
     end if
     if ( abs(clplus) .gt. abs(crminus) ) then
-      dfz(i,j,k) = clplus
+      dfz(i,j,k,l) = clplus
     else
-      dfz(i,j,k) = crminus
+      dfz(i,j,k,l) = crminus
     end if
   endif
 
   ! Ditto if f<0.
-  if ( f(i,j,k) .lt. 0 ) then
+  if ( f(i,j,k,l) .lt. 0 ) then
     if ( abs(alminus) .gt. abs(arplus) ) then
-      dfx(i,j,k) = alminus
+      dfx(i,j,k,l) = alminus
     else
-      dfx(i,j,k) = arplus
+      dfx(i,j,k,l) = arplus
     end if
     if ( abs(blminus) .gt. abs(brplus) ) then
-      dfy(i,j,k) = blminus
+      dfy(i,j,k,l) = blminus
     else
-      dfy(i,j,k) = brplus
+      dfy(i,j,k,l) = brplus
     end if
     if ( abs(clminus) .gt. abs(crplus) ) then
-      dfz(i,j,k) = clminus
+      dfz(i,j,k,l) = clminus
     else
-      dfz(i,j,k) = crplus
+      dfz(i,j,k,l) = crplus
     end if
   end if
 
 ! If the cell is not active set the derivatives to zero.
 else
-  dfx(i,j,k) = zero
-  dfy(i,j,k) = zero
-  dfz(i,j,k) = zero
+  dfx(i,j,k,l) = zero
+  dfy(i,j,k,l) = zero
+  dfz(i,j,k,l) = zero
 end if
diff --git a/src/include/upwind_shift_second2.h b/src/include/upwind_shift_second2.h
index f5fff46..57b7930 100644
--- a/src/include/upwind_shift_second2.h
+++ b/src/include/upwind_shift_second2.h
@@ -2,106 +2,106 @@
 ! $Header$
 
 ! If this is an active cell...
-if ( eh_mask(i,j,k) .ge. 0 ) then
+if ( eh_mask(i,j,k,l) .ge. 0 ) then
 
   ! If the shift in the x-direction is <0.
   if ( betax(i,j,k) .lt. zero ) then
 
     ! If both of the two nearest cells to the right are active cells...
-    if ( ( eh_mask(i+1,j,k) .ge. 0 ) .and. ( eh_mask(i+2,j,k) .ge. 0 ) ) then
+    if ( ( eh_mask(i+1,j,k,l) .ge. 0 ) .and. ( eh_mask(i+2,j,k,l) .ge. 0 ) ) then
 
       ! Calculate the right handed one sided derivative.
-      dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                           four * f(i+1,j,k) - f(i+2,j,k) )
+      dfx(i,j,k,l) = idx * ( -three * f(i,j,k,l) + &
+                           four * f(i+1,j,k,l) - f(i+2,j,k,l) )
 
     ! Else if only the nearest neighbour to the right is active...
-    else if ( eh_mask(i+1,j,k) .ge. 0 ) then
+    else if ( eh_mask(i+1,j,k,l) .ge. 0 ) then
 
       ! Calculate the centered derivative.
-      dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+      dfx(i,j,k,l) = idx * ( f(i+1,j,k,l) - f(i-1,j,k,l) )
 
     ! Else it must be a boundary cell...
     else
 
       ! So calculate the left handed one sided derivative.
-      dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                           four * f(i-1,j,k) + f(i-2,j,k) )
+      dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - &
+                           four * f(i-1,j,k,l) + f(i-2,j,k,l) )
     end if
 
   ! Else if the shift is >= 0.
   else
 
     ! If both of the two nearest cells to the left are active cells...
-    if ( ( eh_mask(i-1,j,k) .ge. 0 ) .and. ( eh_mask(i-2,j,k) .ge. 0 ) ) then
+    if ( ( eh_mask(i-1,j,k,l) .ge. 0 ) .and. ( eh_mask(i-2,j,k,l) .ge. 0 ) ) then
 
       ! Calculate the left handed one sided derivative.
-      dfx(i,j,k) = idx * ( three * f(i,j,k) - &
-                           four * f(i-1,j,k) + f(i-2,j,k) )
+      dfx(i,j,k,l) = idx * ( three * f(i,j,k,l) - &
+                           four * f(i-1,j,k,l) + f(i-2,j,k,l) )
 
     ! Else if only the nearest neighbour to the left is active...
-    else if ( eh_mask(i-1,j,k) .ge. 0 ) then
+    else if ( eh_mask(i-1,j,k,l) .ge. 0 ) then
 
       ! Calculate the centered derivative.
-      dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+      dfx(i,j,k,l) = idx * ( f(i+1,j,k,l) - f(i-1,j,k,l) )
 
     ! Else it must be a boundary cell...
     else
 
       ! So calculate the left handed one sided derivative.
-      dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
-                           four * f(i+1,j,k) - f(i+2,j,k) )
+      dfx(i,j,k,l) = idx * ( -three * f(i,j,k,l) + &
+                           four * f(i+1,j,k,l) - f(i+2,j,k,l) )
     end if
   end if
 
   ! Ditto for the y-derivative.
   if ( betay(i,j,k) .lt. zero ) then
-    if ( ( eh_mask(i,j+1,k) .ge. 0 ) .and. ( eh_mask(i,j+2,k) .ge. 0 ) ) then
-      dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                           four * f(i,j+1,k) - f(i,j+2,k) )
-    else if ( eh_mask(i,j+1,k) .ge. 0 ) then
-      dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
+    if ( ( eh_mask(i,j+1,k,l) .ge. 0 ) .and. ( eh_mask(i,j+2,k,l) .ge. 0 ) ) then
+      dfy(i,j,k,l) = idy * ( -three * f(i,j,k,l) + &
+                           four * f(i,j+1,k,l) - f(i,j+2,k,l) )
+    else if ( eh_mask(i,j+1,k,l) .ge. 0 ) then
+      dfy(i,j,k,l) = idy * ( f(i,j+1,k,l) - f(i,j-1,k,l) )
     else
-      dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                           four * f(i,j-1,k) + f(i,j-2,k) )
+      dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - &
+                           four * f(i,j-1,k,l) + f(i,j-2,k,l) )
     end if
   else
-    if ( ( eh_mask(i,j-1,k) .ge. 0 ) .and. ( eh_mask(i,j-2,k) .ge. 0 ) ) then
-      dfy(i,j,k) = idy * ( three * f(i,j,k) - &
-                           four * f(i,j-1,k) + f(i,j-2,k) )
-    else if ( eh_mask(i-1,j,k) .ge. 0 ) then
-      dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
+    if ( ( eh_mask(i,j-1,k,l) .ge. 0 ) .and. ( eh_mask(i,j-2,k,l) .ge. 0 ) ) then
+      dfy(i,j,k,l) = idy * ( three * f(i,j,k,l) - &
+                           four * f(i,j-1,k,l) + f(i,j-2,k,l) )
+    else if ( eh_mask(i-1,j,k,l) .ge. 0 ) then
+      dfy(i,j,k,l) = idy * ( f(i,j+1,k,l) - f(i,j-1,k,l) )
     else
-      dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
-                           four * f(i,j+1,k) - f(i,j+2,k) )
+      dfy(i,j,k,l) = idy * ( -three * f(i,j,k,l) + &
+                           four * f(i,j+1,k,l) - f(i,j+2,k,l) )
     end if
   end if
 
   ! Ditto for the z-derivative.
   if ( betaz(i,j,k) .lt. zero ) then
-    if ( ( eh_mask(i,j,k+1) .ge. 0 ) .and. ( eh_mask(i,j,k+2) .ge. 0 ) ) then
-      dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                           four * f(i,j,k+1) - f(i,j,k+2) )
-    else if ( eh_mask(i,j,k+1) .ge. 0 ) then
-      dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
+    if ( ( eh_mask(i,j,k+1,l) .ge. 0 ) .and. ( eh_mask(i,j,k+2,l) .ge. 0 ) ) then
+      dfz(i,j,k,l) = idz * ( -three * f(i,j,k,l) + &
+                           four * f(i,j,k+1,l) - f(i,j,k+2,l) )
+    else if ( eh_mask(i,j,k+1,l) .ge. 0 ) then
+      dfz(i,j,k,l) = idz * ( f(i,j,k+1,l) - f(i,j,k-1,l) )
     else
-      dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                           four * f(i,j,k-1) + f(i,j,k-2) )
+      dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - &
+                           four * f(i,j,k-1,l) + f(i,j,k-2,l) )
     end if
   else
-	    if ( ( eh_mask(i,j,k-1) .ge. 0 ) .and. ( eh_mask(i,j,k-2) .ge. 0 ) ) then
-      dfz(i,j,k) = idz * ( three * f(i,j,k) - &
-                           four * f(i,j,k-1) + f(i,j,k-2) )
-    else if ( eh_mask(i-1,j,k) .ge. 0 ) then
-      dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
+	    if ( ( eh_mask(i,j,k-1,l) .ge. 0 ) .and. ( eh_mask(i,j,k-2,l) .ge. 0 ) ) then
+      dfz(i,j,k,l) = idz * ( three * f(i,j,k,l) - &
+                           four * f(i,j,k-1,l) + f(i,j,k-2,l) )
+    else if ( eh_mask(i-1,j,k,l) .ge. 0 ) then
+      dfz(i,j,k,l) = idz * ( f(i,j,k+1,l) - f(i,j,k-1,l) )
     else
-      dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
-                           four * f(i,j,k+1) - f(i,j,k+2) )
+      dfz(i,j,k,l) = idz * ( -three * f(i,j,k,l) + &
+                           four * f(i,j,k+1,l) - f(i,j,k+2,l) )
     end if
   end if
 
 ! If the cell is not active set the derivatives to zero.
 else
-  dfx(i,j,k) = zero
-  dfy(i,j,k) = zero
-  dfz(i,j,k) = zero
+  dfx(i,j,k,l) = zero
+  dfy(i,j,k,l) = zero
+  dfz(i,j,k,l) = zero
 end if