
      MODULE types
      
!    symbolic names for 4, 2, and 1 byte integers
	
	INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
	INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
	INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)

! and for single and double precision reals
!	INTEGER, PARAMETER :: SP = KIND(1.0)
!	INTEGER, PARAMETER :: DP = KIND(1.0D0)
	INTEGER, PARAMETER :: SP = KIND(1.0D0)

      END MODULE types


	
module hpfprocessors
!HPF$ PROCESSORS PROC(4)
end module hpfprocessors


module diff_interface

INTERFACE
subroutine diff_x(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_x

subroutine diff_y(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_y

subroutine diff_z(nx,ny,nz,h,f,df)
use types
use hpfprocessors
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df
end subroutine diff_z

END INTERFACE
end module diff_interface







subroutine diff_x(nx,ny,nz,h,f,df)
use types
use hpfprocessors

implicit none
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df

real(SP)                    :: h2,hi
integer(I4B)                :: i,j,k

h2 = 0.5d0/h
hi = 1.0d0/h
	
! x derivative


!        interior
          FORALL (k=2:nz-1,j=2:ny-1,i=2:nx-1) &
          df(i,j,k) = h2 * ( f(i+1,j,k) - f(i-1,j,k) )

!        outer bounaries
	  FORALL (k=1:nz,j=1:ny)         
          df(1,j,k) = h2 * ( -3.0d0*f(1,j,k)+ 4.0d0*f(1+1,j,k) - f(1+2,j,k))
          df(nx,j,k) = h2 * (3.0d0*f(nx,j,k)- 4.0d0*f(nx-1,j,k) + f(nx-2,j,k))
          END FORALL

          FORALL (k=1:nz,i=2:nx-1)
          df(i,1,k) = h2 *(f(i+1,1,k)-f(i-1,1,k) )
          df(i,ny,k) = h2 *(f(i+1,ny,k)-f(i-1,ny,k) )
          END FORALL

          FORALL (j=1:ny,i=2:nx-1)
          df(i,j,1) = h2 *(f(i+1,j,1)-f(i-1,j,1) )
          df(i,j,nz) = h2 *(f(i+1,j,nz)-f(i-1,j,nz) )
          END FORALL


      return
      end

subroutine diff_y(nx,ny,nz,h,f,df)
use types
use hpfprocessors

implicit none
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df

real(SP)                    :: h2,hi
integer(I4B)                :: i,j,k

h2 = 0.5d0/h
hi = 1.0d0/h

! Y derivative

!       interior

	  FORALL ( k=2:nz-1,j=2:ny-1,i=2:nx-1) &
          df(i,j,k) = h2 * ( f(i,j+1,k) - f(i,j-1,k) )

!       outer boundaries

	  FORALL (k=1:nz,i=1:nx)
          df(i,1,k) = h2 * ( -3.0d0*f(i,1,k) + 4.0d0*f(i,1+1,k) - f(i,1+2,k))
          df(i,ny,k) = h2 * (3.0d0*f(i,ny,k)- 4.0d0*f(i,ny-1,k) + f(i,ny-2,k))
	  END FORALL

	  FORALL (k=1:nz,j=2:ny-1)
          df(1,j,k) = h2 *(f(1,j+1,k)-f(1,j-1,k) )
          df(nx,j,k) = h2 *(f(nx,j+1,k)-f(nx,j-1,k) )
	  END FORALL

	  FORALL (i=1:nx,j=2:ny-1)
          df(i,j,1) = h2 *(f(i,j+1,1)-f(i,j-1,1) )
          df(i,j,nz) = h2 *(f(i,j+1,nz)-f(i,j-1,nz) )
	  END FORALL


      return
      end

subroutine diff_z(nx,ny,nz,h,f,df)
use types
use hpfprocessors

implicit none
integer(I4B)                :: nx,ny,nz
real(SP) h
real(SP), dimension(:,:,:), intent(IN)  :: f
real(SP), dimension(:,:,:), intent(OUT) :: df

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
!HPF$ ALIGN WITH *f                      :: df

real(SP)                    :: h2,hi
integer(I4B)                :: i,j,k

real(SP), dimension(nx,ny,2)       :: tmp1
real(SP), dimension(nx,ny,nz-1:nz) :: tmp2
!HPF$ ALIGN (i,j,k) WITH f(i,j,k)  :: tmp1,tmp2


h2 = 0.5d0/h
hi = 1.0d0/h
! Z derivative

!       interior

          FORALL ( k=2:nz-1,j=2:ny-1,i=2:nx-1)   & 
          df(i,j,k) = h2 * ( f(i,j,k+1) - f(i,j,k-1) )

!       outer boundaries

!	  FORALL (j=1:ny,i=1:nx,k=1:1)           & 
!          df(i,j,k) = h2 * ( -3.0d0*f(i,j,k)+ 4.0d0*f(i,j,k+1) - f(i,j,k+2))

          tmp1(:,:,1:2) =f(:,:,2:3)-f(:,:,1:2)
          df(:,:,1:1)=h2*(3.0d0*tmp1(:,:,1:1)-tmp1(:,:,2:2))
          
!	  FORALL (j=1:ny,i=1:nx,k=nz:nz)           & 
!          df(i,j,k) = h2 * (3.0d0*f(i,j,k)- 4.0d0*f(i,j,k-1) + f(i,j,k-2))

          tmp2(:,:,nz-1:nz) =f(:,:,nz-1:nz)-f(:,:,nz-2:nz-1)
          df(:,:,nz:nz)=h2*(3.0d0*tmp2(:,:,nz:nz)-tmp2(:,:,nz-1:nz-1))
          
	  FORALL (i=1:nx,k=2:nz-1)
          df(i,1,k) = h2 *(f(i,1,k+1)-f(i,1,k-1) )
          df(i,ny,k) = h2 *(f(i,ny,k+1)-f(i,ny,k-1) )
	  END FORALL

	  FORALL (j=1:ny,k=2:nz-1)
          df(1,j,k) = h2 *(f(1,j,k+1)-f(1,j,k-1) )
          df(nx,j,k) = h2 *(f(nx,j,k+1)-f(nx,j,k-1) )
	  END FORALL


      return
      end










