      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_k(gup11, gup12, gup13, gup22, gup23, gup33, k11,&
     & k12, k13, k21, k22, k23, k31, k32, k33, k11m1, k12m1, k13m1, k21m&
     &1, k22m1, k23m1, k31m1, k32m1, k33m1, c111, c112, c113, c122, c123&
     &, c133, c211, c212, c213, c222, c223, c233, c311, c312, c313, c322&
     &, c323, c333, r11, r12, r13, r22, r23, r33, alpha, beta1, beta2, b&
     &eta3, nx, ny, nz, h, dt)
      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(inout):: k11m1,k12m1,k13m1,k21&
     &m1,k22m1,k23m1,k31m1,k32m1,k33m1 
      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(in):: r11,r12,r13,r22,r23,r33 

      real (SP),  dimension(:,:,:),intent(in):: alpha,beta1,beta2,beta3 

!HPF$      distribute * (*,*,block) onto * PROC :: k11
!HPF$      align with * k11 :: k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * k11 :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1&
!HPF$&,k32m1,k33m1
!HPF$      align with * k11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * k11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * k11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with * k11 :: r11,r12,r13,r22,r23,r33
!HPF$      align with * k11 :: gup11,gup12,gup13,gup22,gup23,gup33
!HPF$      align with * k11 :: alpha,beta1,beta2,beta3
      integer (I4B):: nx,ny,nz
      real (SP):: h,dt
      real (SP),  dimension(size (k11,1),size (k11,2),size (k11,3)):: t1&
     &1,t12,t13,t21,t22,t23,t31,t32,t33 
!HPF$      align with k11 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
      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,c&
     &123,c133 
         real (SP),  dimension(:,:,:),intent(in):: c211,c212,c213,c222,c&
     &223,c233 
         real (SP),  dimension(:,:,:),intent(in):: c311,c312,c313,c322,c&
     &323,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
         end subroutine

         subroutine evolve_k2(gup11, gup12, gup13, gup22, gup23, gup33, &
     &t11, t12, t13, t21, t22, t23, t31, t32, t33)
         use types
         use hpfprocessors
         real (SP),  dimension(:,:,:),intent(in):: gup11,gup12,gup13,gup&
     &22,gup23,gup33 
         real (SP),  dimension(:,:,:),intent(inout):: t11,t12,t13,t21,t2&
     &2,t23,t31,t32,t33 
!HPF$      distribute * (*,*,block) onto * PROC :: gup11
!HPF$      align with * gup11 :: gup12,gup13,gup22,gup23,gup33
!HPF$      align with * gup11 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
         end subroutine

         subroutine evolve_k3a(k11, k12, k13, k21, k22, k23, k31, k32, k&
     &33, beta1, beta2, beta3, t11, t12, t13, t21, t22, t23, t31, t32, t&
     &33, nx, ny, nz, h)
         use types
         use hpfprocessors
         use diff_interface
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(in):: beta1,beta2,beta3 
         real (SP),  dimension(:,:,:),intent(inout):: t11,t12,t13,t21,t2&
     &2,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
         end subroutine

         subroutine evolve_k3b(k11, k12, k13, k21, k22, k23, k31, k32, k&
     &33, beta1, beta2, beta3, t11, t12, t13, t21, t22, t23, t31, t32, t&
     &33, nx, ny, nz, h)
         use types
         use hpfprocessors
         use diff_interface
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(in):: beta1,beta2,beta3 
         real (SP),  dimension(:,:,:),intent(inout):: t11,t12,t13,t21,t2&
     &2,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
         end subroutine

         subroutine evolve_k4(gup11, gup12, gup13, gup22, gup23, gup33, &
     &k11, k12, k13, k21, k22, k23, k31, k32, k33, r11, r12, r13, r22, r&
     &23, r33, t11, t12, t13, t21, t22, t23, t31, t32, t33, alpha)
         use types
         use hpfprocessors
         real (SP),  dimension(:,:,:),intent(in):: gup11,gup12,gup13,gup&
     &22,gup23,gup33 
         real (SP),  dimension(:,:,:),intent(in):: k11,k12,k13,k21,k22,k&
     &23,k31,k32,k33 
         real (SP),  dimension(:,:,:),intent(in):: r11,r12,r13,r22,r23,r&
     &33 
         real (SP),  dimension(:,:,:),intent(inout):: t11,t12,t13,t21,t2&
     &2,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
         end subroutine

         subroutine evolve_k5(k11m1, k12m1, k13m1, k21m1, k22m1, k23m1, &
     &k31m1, k32m1, k33m1, t11, t12, t13, t21, t22, t23, t31, t32, t33, &
     &dt)
         use types
         use hpfprocessors
         real (SP),  dimension(:,:,:),intent(inout):: k11m1,k12m1,k13m1,&
     &k21m1,k22m1,k23m1,k31m1,k32m1,k33m1 
         real (SP),  dimension(:,:,:),intent(inout):: t11,t12,t13,t21,t2&
     &2,t23,t31,t32,t33 
!HPF$      distribute * (*,*,block) onto * PROC :: k11m1
!HPF$      align with * k11m1 :: k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32&
!HPF$&m1,k33m1
!HPF$      align with * k11m1 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
         real (SP):: dt
         end subroutine

      end interface
      call 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)
      call evolve_k2 (gup11,gup12,gup13,gup22,gup23,gup33,t11,t12,t13,t2&
     &1,t22,t23,t31,t32,t33)
      call evolve_k3a (k11,k12,k13,k21,k22,k23,k31,k32,k33,beta1,beta2,b&
     &eta3,t11,t12,t13,t12,t22,t23,t31,t32,t33,nx,ny,nz,h)
      call evolve_k3b (k11,k12,k13,k21,k22,k23,k31,k32,k33,beta1,beta2,b&
     &eta3,t11,t12,t13,t12,t22,t23,t31,t32,t33,nx,ny,nz,h)
      call evolve_k4 (gup11,gup12,gup13,gup22,gup23,gup33,k11,k12,k13,k2&
     &1,k22,k23,k31,k32,k33,r11,r12,r13,r22,r23,r33,t11,t12,t13,t21,t22,&
     &t23,t31,t32,t33,alpha)
      call evolve_k5 (k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32m1,k3&
     &3m1,t11,t12,t13,t21,t22,t23,t31,t32,t33,dt)
      return
      end subroutine

