Merge pull request #400 from matt-frey/array-index

Use named indices for n-rank arrays

src/3d/Makefile.am

@@ -6,6 +6,7 @@ AM_FCFLAGS =				\
 
 bin_PROGRAMS = epic3d
 epic3d_SOURCES = 				\
+	utils/dimensions.f90			\
 	utils/deriv1d.f90			\
 	utils/options.f90			\
 	utils/parameters.f90			\

src/3d/fields/fields.f90

@@ -3,6 +3,7 @@
 !     and functions.
 ! =============================================================================
 module fields
+    use dimensions, only : n_dim, I_X, I_Y, I_Z
     use parameters, only : dx, dxi, extent, lower, nx, ny, nz
     use constants, only : zero
     implicit none
@@ -16,7 +17,7 @@ module fields
     ! and from 0 to ny-1 in y
     double precision, allocatable, dimension(:, :, :, :) :: &
         velog,     &   ! velocity vector field (u, v, w)
-        vortg,     &   ! vorticity vector field (\omegax, \omegay, \omegaz)
+        vortg,     &   ! vorticity vector field (\xi, \eta, \zeta)
         vtend,     &   ! vorticity tendency
         velgradg       ! velocity gradient tensor
                        ! ordering: du/dx, du/dy,
@@ -25,9 +26,9 @@ module fields
                        ! the derivatives dv/dx, du/dz, dv/dz and dw/dz
                        ! are calculated on the fly with vorticity
                        ! or the assumption of incompressibility (du/dx + dv/dy + dw/dz = 0):
-                       !    dv/dx = \omegaz + du/dy
-                       !    du/dz = \omegay + dw/dx
-                       !    dv/dz = dw/dy - \omegax
+                       !    dv/dx = \zeta + du/dy
+                       !    du/dz = \eta + dw/dx
+                       !    dv/dz = dw/dy - \xi
                        !    dw/dz = - (du/dx + dv/dy)
 
     double precision, allocatable, dimension(:, :, :) :: &
@@ -44,6 +45,14 @@ module fields
         nparg,     &   ! number of parcels per grid box
         nsparg         ! number of small parcels per grid box
 
+    ! velocity strain indices
+    integer, parameter ::
+                        , I_DUDX = 1 & ! index for du/dx strain component
+                        , I_DUDY = 2 & ! index for du/dy strain component
+                        , I_DVDY = 3 & ! index for dv/dy strain component
+                        , I_DWDX = 4 & ! index for dw/dx strain component
+                        , I_DWDY = 5   ! index for dw/dy strain component
+
     contains
 
         ! Allocate all fields
@@ -52,7 +61,7 @@ subroutine field_alloc
                 return
             endif
 
-            allocate(velog(-1:nz+1, 0:ny-1, 0:nx-1, 3))
+            allocate(velog(-1:nz+1, 0:ny-1, 0:nx-1, n_dim))
             allocate(velgradg(-1:nz+1, 0:ny-1, 0:nx-1, 5))
 
             allocate(volg(-1:nz+1, 0:ny-1, 0:nx-1))
@@ -61,9 +70,9 @@ subroutine field_alloc
             allocate(sym_volg(-1:nz+1, 0:ny-1, 0:nx-1))
 #endif
 
-            allocate(vortg(-1:nz+1, 0:ny-1, 0:nx-1, 3))
+            allocate(vortg(-1:nz+1, 0:ny-1, 0:nx-1, n_dim))
 
-            allocate(vtend(-1:nz+1, 0:ny-1, 0:nx-1, 3))
+            allocate(vtend(-1:nz+1, 0:ny-1, 0:nx-1, n_dim))
 
             allocate(tbuoyg(-1:nz+1, 0:ny-1, 0:nx-1))
 
@@ -99,12 +108,12 @@ subroutine field_default
         ! @param[out] i lower, zonal cell index
         ! @param[out] j lower, vertical cell index
         pure subroutine get_index(pos, i, j, k)
-            double precision, intent(in)  :: pos(3)
+            double precision, intent(in)  :: pos(n_dim)
             integer,          intent(out) :: i, j, k
 
-            i = floor((pos(1) - lower(1)) * dxi(1))
-            j = floor((pos(2) - lower(2)) * dxi(2))
-            k = floor((pos(3) - lower(3)) * dxi(3))
+            i = floor((pos(I_X) - lower(I_X)) * dxi(I_X))
+            j = floor((pos(I_Y) - lower(I_Y)) * dxi(I_Y))
+            k = floor((pos(I_Z) - lower(I_Z)) * dxi(I_Z))
         end subroutine get_index
 
 
@@ -133,11 +142,11 @@ end subroutine periodic_index_shift
         ! @param[out] pos position of (i, j, k) in the domain
         pure subroutine get_position(i, j, k, pos)
             integer,          intent(in)  :: i, j, k
-            double precision, intent(out) :: pos(3)
+            double precision, intent(out) :: pos(n_dim)
 
-            pos(1) = lower(1) + i * dx(1)
-            pos(2) = lower(2) + j * dx(2)
-            pos(3) = lower(3) + k * dx(3)
+            pos(I_X) = lower(I_X) + i * dx(I_X)
+            pos(I_Y) = lower(I_Y) + j * dx(I_Y)
+            pos(I_Z) = lower(I_Z) + k * dx(I_Z)
 
         end subroutine get_position