      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_g(g11, g12, g13, g22, g23, g33, g11m1, g12m1, g1&
     &3m1, g22m1, g23m1, g33m1, c111, c112, c113, c122, c123, c133, c211&
     &, c212, c213, c222, c223, c233, c311, c312, c313, c322, c323, c333&
     &, k11, k12, k13, k21, k22, k23, k31, k32, k33, alpha, beta1, beta2&
     &, beta3, h, dt, nx, ny, nz)
      use types
      use hpfprocessors
      use diff_interface
      real (SP),  dimension(:,:,:),intent(in):: g11,g12,g13,g22,g23,g33 

      real (SP),  dimension(:,:,:),intent(out):: g11m1,g12m1,g13m1,g22m1&
     &,g23m1,g33m1 
      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):: k11,k12,k13,k21,k22,k23,&
     &k31,k32,k33 
      real (SP),  dimension(:,:,:),intent(in):: alpha,beta1,beta2,beta3 

      integer (I4B):: nx,ny,nz
      real (SP):: h,dt
!HPF$      distribute g11 * (*,*,block) onto * PROC
!HPF$      align with * g11 :: g12,g13,g22,g23,g33
!HPF$      align with * g11 :: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
!HPF$      align with * g11 :: c111,c112,c113,c122,c123,c133
!HPF$      align with * g11 :: c211,c212,c213,c222,c223,c233
!HPF$      align with * g11 :: c311,c312,c313,c322,c323,c333
!HPF$      align with * g11 :: k11,k12,k13,k21,k22,k23,k31,k32,k33
!HPF$      align with * g11 :: alpha,beta1,beta2,beta3
      real (SP),  dimension(size (g11,1),size (g11,2),size (g11,3)):: t1&
     &1,t12,t13,t21,t22,t23,t31,t32,t33 
!HPF$      align with g11 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
      call diff_x (nx,ny,nz,h,beta1,t11)
      call diff_y (nx,ny,nz,h,beta1,t12)
      call diff_z (nx,ny,nz,h,beta1,t13)
      call diff_x (nx,ny,nz,h,beta2,t21)
      call diff_y (nx,ny,nz,h,beta2,t22)
      call diff_z (nx,ny,nz,h,beta2,t23)
      call diff_x (nx,ny,nz,h,beta3,t31)
      call diff_y (nx,ny,nz,h,beta3,t32)
      call diff_z (nx,ny,nz,h,beta3,t33)
      t11 = t11+c111*beta1+c112*beta2+c113*beta3
      t12 = t12+c112*beta1+c122*beta2+c123*beta3
      t13 = t13+c113*beta1+c123*beta2+c133*beta3
      t21 = t21+c211*beta1+c212*beta2+c213*beta3
      t22 = t22+c212*beta1+c222*beta2+c223*beta3
      t23 = t23+c213*beta1+c223*beta2+c233*beta3
      t31 = t31+c311*beta1+c312*beta2+c313*beta3
      t32 = t32+c312*beta1+c322*beta2+c323*beta3
      t33 = t33+c313*beta1+c323*beta2+c333*beta3
      t11 = t11*g11+t12*g12+t13*g13
      t12 = t11*g12+t12*g22+t13*g23
      t13 = t11*g13+t12*g23+t13*g33
      t21 = t21*g11+t22*g12+t23*g13
      t22 = t21*g12+t22*g22+t23*g23
      t23 = t21*g13+t22*g23+t23*g33
      t31 = t31*g11+t32*g12+t33*g13
      t32 = t31*g12+t32*g22+t33*g23
      t33 = t31*g13+t32*g23+t33*g33
      t11 = t11+t11
      t12 = t12+t21
      t13 = t13+t31
      t22 = t22+t22
      t23 = t23+t32
      t33 = t33+t33
      t21 = t12
      t31 = t13
      t32 = t23
      t11 = t11-(2.0d0*alpha*(g11*k11+g12*k21+g13*k31))
      t12 = t12-(2.0d0*alpha*(g11*k12+g12*k22+g13*k32))
      t13 = t13-(2.0d0*alpha*(g11*k13+g12*k23+g13*k33))
      t21 = t21-(2.0d0*alpha*(g12*k11+g22*k21+g23*k31))
      t22 = t22-(2.0d0*alpha*(g12*k12+g22*k22+g23*k32))
      t23 = t23-(2.0d0*alpha*(g12*k13+g22*k23+g23*k33))
      t31 = t31-(2.0d0*alpha*(g13*k11+g23*k21+g33*k31))
      t32 = t32-(2.0d0*alpha*(g13*k12+g23*k22+g33*k32))
      t33 = t33-(2.0d0*alpha*(g13*k13+g23*k23+g33*k33))
      t11 = 2.0d0*dt*t11+g11m1
      t12 = 2.0d0*dt*t12+g12m1
      t13 = 2.0d0*dt*t13+g13m1
      t22 = 2.0d0*dt*t22+g22m1
      t23 = 2.0d0*dt*t23+g23m1
      t33 = 2.0d0*dt*t33+g33m1
      g11m1 = t11
      g12m1 = t12
      g13m1 = t13
      g22m1 = t22
      g23m1 = t23
      g33m1 = t33
      return
      end subroutine

