
      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_k1(alpha,                                       &
                    c111,c112,c113,c122,c123,c133,                &
                    c211,c212,c213,c222,c223,c233,                &
                    c311,c312,c313,c322,c323,c333,                &
                    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) :: alpha
real(SP), dimension(:,:,:), intent(IN) :: c111,c112,c113,c122,c123,c133
real(SP), dimension(:,:,:), intent(IN) :: c211,c212,c213,c222,c223,c233
real(SP), dimension(:,:,:), intent(IN) :: c311,c312,c313,c322,c323,c333
real(SP), dimension(:,:,:), intent(OUT):: t11,t12,t13,t21,t22,t23,t31,t32,t33

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: alpha 
!HPF$ ALIGN WITH *alpha :: c111,c112,c113,c122,c123,c133
!HPF$ ALIGN WITH *alpha :: c211,c212,c213,c222,c223,c233
!HPF$ ALIGN WITH *alpha :: c311,c312,c313,c322,c323,c333
!HPF$ ALIGN WITH *alpha :: t11,t12,t13,t21,t22,t23,t31,t32,t33

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

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

!     Compute D^a(D_b(alpha)) terms

      call diff_x(nx,ny,nz,h,alpha,temp1)
      call diff_y(nx,ny,nz,h,alpha,temp2)
      call diff_z(nx,ny,nz,h,alpha,temp3)

!   Calculate second derivatives of alpha

      call diff_x(nx,ny,nz,h,temp1,t11)
      call diff_y(nx,ny,nz,h,temp1,t12)
      call diff_z(nx,ny,nz,h,temp1,t13)
      call diff_y(nx,ny,nz,h,temp2,t22)
      call diff_z(nx,ny,nz,h,temp2,t23)
      call diff_z(nx,ny,nz,h,temp3,t33)
      t21 = t12
      t31 = t13
      t32 = t23

! form covariant derivative D_a D_b alpha
! D_a D_b alpha = t{ab} - c1{ab}*temp1 - c2{ab}*temp2 - c3{ab}*temp3

      t11 = t11 - c111 * temp1 - c211 * temp2 - c311 * temp3
      t12 = t12 - c112 * temp1 - c212 * temp2 - c312 * temp3
      t13 = t13 - c113 * temp1 - c213 * temp2 - c313 * temp3
      t22 = t22 - c122 * temp1 - c222 * temp2 - c322 * temp3
      t23 = t23 - c123 * temp1 - c223 * temp2 - c323 * temp3
      t33 = t33 - c133 * temp1 - c233 * temp2 - c333 * temp3
      t21 = t12
      t31 = t13
      t32 = t23

return
end
