
      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_k3a(k11,k12,k13,k21,k22,k23,k31,k32,k33,          &
                     beta1,beta2,beta3,                             &
                     t11,t12,t13,t21,t22,t23,t31,t32,t33,           &
                     nx,ny,nz,h)

use types
use hpfprocessors
use diff_interface

real(SP), dimension(:,:,:), intent(IN) :: k11,k12,k13,k21,k22,k23,k31,k32,k33
real(SP), dimension(:,:,:), intent(IN) :: beta1,beta2,beta3
real(SP), dimension(:,:,:), intent(INOUT):: t11,t12,t13,t21,t22,t23,t31,t32,t33

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: k11 
!HPF$ ALIGN WITH *k11 :: k12,k13,k21,k22,k23,k31,k32,k33
!HPF$ ALIGN WITH *k11 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
!HPF$ ALIGN WITH *k11 :: beta1,beta2,beta3

integer(I4B) :: nx,ny,nz
real(SP)     :: h

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


!     Compute Lie_beta(K^a_b) terms

!     Lie_beta(K^1_1)
!k21
      call diff_x(nx,ny,nz,h,k11,temp1)
      call diff_y(nx,ny,nz,h,k11,temp2)
      call diff_z(nx,ny,nz,h,k11,temp3)

      t11 = t11 + beta1 * temp1  + beta2 * temp2 + beta3 * temp3

!k22
      call diff_x(nx,ny,nz,h,beta1,temp1)
      call diff_y(nx,ny,nz,h,beta1,temp2)
      call diff_z(nx,ny,nz,h,beta1,temp3)

      t11 = t11 - k11 * temp1 - k21 * temp2 - k31 * temp3

!k23
      call diff_y(nx,ny,nz,h,beta1,temp1)
      call diff_y(nx,ny,nz,h,beta2,temp2)
      call diff_y(nx,ny,nz,h,beta3,temp3)

      t11 = t11 + k11 * temp1 + k12 * temp2 + k13 * temp3

!     Lie_beta(K^1_2)

!k24
      call diff_x(nx,ny,nz,h,k12,temp1)
      call diff_y(nx,ny,nz,h,k12,temp2)
      call diff_z(nx,ny,nz,h,k12,temp3)

      t12 = t12 + beta1 * temp1  + beta2 * temp2 + beta3 * temp3

!k25
      call diff_x(nx,ny,nz,h,beta1,temp1)
      call diff_y(nx,ny,nz,h,beta1,temp2)
      call diff_z(nx,ny,nz,h,beta1,temp3)

      t12 = t12 - k12 * temp1 - k22 * temp2 - k32 * temp3

!k26
      call diff_y(nx,ny,nz,h,beta1,temp1)
      call diff_y(nx,ny,nz,h,beta2,temp2)
      call diff_y(nx,ny,nz,h,beta3,temp3)

      t12 = t12 + k11 * temp1 + k12 * temp2 + k13 * temp3

!     Lie_beta(K^1_3)

!k27
      call diff_x(nx,ny,nz,h,k13,temp1)
      call diff_y(nx,ny,nz,h,k13,temp2)
      call diff_z(nx,ny,nz,h,k13,temp3)

      t13 = t13 + beta1 * temp1  + beta2 * temp2 + beta3 * temp3

!k28
      call diff_x(nx,ny,nz,h,beta1,temp1)
      call diff_y(nx,ny,nz,h,beta1,temp2)
      call diff_z(nx,ny,nz,h,beta1,temp3)

      t13 = t13 - k13 * temp1 - k23 * temp2 - k33 * temp3

!k29
      call diff_z(nx,ny,nz,h,beta1,temp1)
      call diff_z(nx,ny,nz,h,beta2,temp2)
      call diff_z(nx,ny,nz,h,beta3,temp3)

      t13 = t13 + k11 * temp1 + k12 * temp2 + k13 * temp3


return
end




