
      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







! evolve_g.f 
!-----------------------------------------------------------------------

subroutine evolve_g(g11,g12,g13,g22,g23,g33, &
                    g11m1,g12m1,g13m1,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)) :: & 
     t11,t12,t13,t21,t22,t23,t31,t32,t33

!HPF$ ALIGN WITH  g11 ::  t11,t12,t13,t21,t22,t23,t31,t32,t33


!--------- evolve_g1 ----------------------------------------------------
!   Compute the derivative of beta's 

      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)

!--------- evolve_g2 ---------------------------------------------------


     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
      
         

! Lower second index of Tij

!    tij = tik * gkj => store tik = temp1k for every i,j

     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


! -------- evolve_g34 --------------------------------------------------


     t11 = t11 + t11      
     t12 = t12 + t21      
     t13 = t13 + t31      
     t22 = t22 + t22      
     t23 = t23 + t32      
     t33 = t33 + t33      
     t21 = t12
     t31 = t13
     t32 = t23

! -------- evolve_g5 ----------------------------------------------------
!    tij = tij - 2*alph(1)%d*(g1i*K1j+g2i*K2j+g3i*K3j)
! but allow for g21 = g12 , g31 = g13, g32 = g23

      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 )

! ------- evolve_g6 ----------------------------------------------------

      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

! here t11 is g11p1 , etc.

      g11m1 = t11
      g12m1 = t12
      g13m1 = t13
      g22m1 = t22
      g23m1 = t23
      g33m1 = t33


return 
end subroutine evolve_g
























