
      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 evolve_k4(gup11,gup12,gup13,gup22,gup23,gup33,         &
                    k11,k12,k13,k21,k22,k23,k31,k32,k33,          &
                    r11,r12,r13,r22,r23,r33,                      & 
                    t11,t12,t13,t21,t22,t23,t31,t32,t33,          &
                    alpha)

use types
use hpfprocessors

real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
real(SP), dimension(:,:,:), intent(IN) :: k11,k12,k13,k21,k22,k23,k31,k32,k33
real(SP), dimension(:,:,:), intent(IN) :: r11,r12,r13,r22,r23,r33
real(SP), dimension(:,:,:), intent(INOUT):: t11,t12,t13,t21,t22,t23,t31,t32,t33
real(SP), dimension(:,:,:), intent(IN) :: alpha

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: k11 
!HPF$ ALIGN WITH *k11 :: k12,k13,k21,k22,k23,k31,k32,k33
!HPF$ ALIGN WITH *k11 :: r11,r12,r13,r22,r23,r33
!HPF$ ALIGN WITH *k11 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
!HPF$ ALIGN WITH *k11 :: gup11,gup12,gup13,gup22,gup23,gup33 
!HPF$ ALIGN WITH *k11 :: alpha

real(SP),dimension(size(k11,1),size(k11,2),size(k11,3)) :: temp
!HPF$ ALIGN WITH k11 :: temp


!     Compute alpha*(R^a_b + K*K^a_b) terms

      temp = k11 + k22 + k33

      t11 = t11 + alpha*( (gup11*r11 +gup12*r12 +gup13*r13) +temp*k11 )
      t12 = t12 + alpha*( (gup11*r12 +gup12*r22 +gup13*r23) +temp*k12 )
      t13 = t13 + alpha*( (gup11*r13 +gup12*r23 +gup13*r33) +temp*k13 )
      t21 = t21 + alpha*( (gup12*r11 +gup22*r12 +gup23*r13) +temp*k21 )
      t22 = t22 + alpha*( (gup12*r12 +gup22*r22 +gup23*r23) +temp*k22 )
      t23 = t23 + alpha*( (gup12*r13 +gup22*r23 +gup23*r33) +temp*k23 )
      t31 = t31 + alpha*( (gup13*r11 +gup23*r12 +gup33*r13) +temp*k31 )
      t32 = t32 + alpha*( (gup13*r12 +gup23*r22 +gup33*r23) +temp*k32 )
      t33 = t33 + alpha*( (gup13*r13 +gup23*r23 +gup33*r33) +temp*k33 )

return
end




