      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_k3b(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,t&
     &23,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)):: te&
     &mp1,temp2,temp3 
!HPF$      align with k11 :: temp1,temp2,temp3
      call diff_x (nx,ny,nz,h,k21,temp1)
      call diff_y (nx,ny,nz,h,k21,temp2)
      call diff_z (nx,ny,nz,h,k21,temp3)
      t21 = t21+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta2,temp1)
      call diff_y (nx,ny,nz,h,beta2,temp2)
      call diff_z (nx,ny,nz,h,beta2,temp3)
      t21 = ((t21-(k11*temp1))-(k21*temp2))-(k31*temp3)
      call diff_x (nx,ny,nz,h,beta1,temp1)
      call diff_x (nx,ny,nz,h,beta2,temp2)
      call diff_x (nx,ny,nz,h,beta3,temp3)
      t21 = t21+k21*temp1+k22*temp2+k23*temp3
      call diff_x (nx,ny,nz,h,k22,temp1)
      call diff_y (nx,ny,nz,h,k22,temp2)
      call diff_z (nx,ny,nz,h,k22,temp3)
      t22 = t22+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta2,temp1)
      call diff_y (nx,ny,nz,h,beta2,temp2)
      call diff_z (nx,ny,nz,h,beta2,temp3)
      t22 = ((t22-(k12*temp1))-(k22*temp2))-(k32*temp3)
      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)
      t22 = t22+k21*temp1+k22*temp2+k23*temp3
      call diff_x (nx,ny,nz,h,k23,temp1)
      call diff_y (nx,ny,nz,h,k23,temp2)
      call diff_z (nx,ny,nz,h,k23,temp3)
      t23 = t23+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta2,temp1)
      call diff_y (nx,ny,nz,h,beta2,temp2)
      call diff_z (nx,ny,nz,h,beta2,temp3)
      t23 = ((t23-(k13*temp1))-(k23*temp2))-(k33*temp3)
      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)
      t23 = t23+k21*temp1+k22*temp2+k23*temp3
      call diff_x (nx,ny,nz,h,k31,temp1)
      call diff_y (nx,ny,nz,h,k31,temp2)
      call diff_z (nx,ny,nz,h,k31,temp3)
      t31 = t31+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta3,temp1)
      call diff_y (nx,ny,nz,h,beta3,temp2)
      call diff_z (nx,ny,nz,h,beta3,temp3)
      t31 = ((t31-(k11*temp1))-(k21*temp2))-(k31*temp3)
      call diff_x (nx,ny,nz,h,beta1,temp1)
      call diff_x (nx,ny,nz,h,beta2,temp2)
      call diff_x (nx,ny,nz,h,beta3,temp3)
      t31 = t31+k31*temp1+k32*temp2+k33*temp3
      call diff_x (nx,ny,nz,h,k32,temp1)
      call diff_y (nx,ny,nz,h,k32,temp2)
      call diff_z (nx,ny,nz,h,k32,temp3)
      t32 = t32+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta3,temp1)
      call diff_y (nx,ny,nz,h,beta3,temp2)
      call diff_z (nx,ny,nz,h,beta3,temp3)
      t32 = ((t32-(k12*temp1))-(k22*temp2))-(k32*temp3)
      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)
      t32 = t32+k31*temp1+k32*temp2+k33*temp3
      call diff_x (nx,ny,nz,h,k33,temp1)
      call diff_y (nx,ny,nz,h,k33,temp2)
      call diff_z (nx,ny,nz,h,k33,temp3)
      t33 = t33+beta1*temp1+beta2*temp2+beta3*temp3
      call diff_x (nx,ny,nz,h,beta3,temp1)
      call diff_y (nx,ny,nz,h,beta3,temp2)
      call diff_z (nx,ny,nz,h,beta3,temp3)
      t33 = ((t33-(k13*temp1))-(k23*temp2))-(k33*temp3)
      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)
      t33 = t33+k31*temp1+k32*temp2+k33*temp3
      return
      end subroutine

