      module types
      integer ,  parameter:: I4B=SELECTED_INT_KIND (9) 
      integer ,  parameter:: I2B=SELECTED_INT_KIND (4) 
      integer ,  parameter:: I1B=SELECTED_INT_KIND (2) 
      integer ,  parameter:: SP=KIND (1.0D0) 
      end module


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


      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

         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

         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

      end interface
      end module


      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,t&
     &23,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)):: te&
     &mp 
!HPF$      align with k11 :: temp
      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 subroutine

