
      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







subroutine evolve_k(gup11,gup12,gup13,gup22,gup23,gup33,          &
                    k11,k12,k13,k21,k22,k23,k31,k32,k33,          &
                    k11m1,k12m1,k13m1,k21m1,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,beta3,                     &
                    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,k21m1,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,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)) :: &
                                         t11,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,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(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 evolve_k1

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,gup22,gup23,gup33
real(SP), dimension(:,:,:), intent(INOUT):: t11,t12,t13,t21,t22,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 evolve_k2

subroutine evolve_k3a(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,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 evolve_k3a

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,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 evolve_k3b

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,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 evolve_k4

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,t22,t23,t31,t32,t33

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: k11m1 
!HPF$ ALIGN WITH *k11m1 :: k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32m1,k33m1
!HPF$ ALIGN WITH *k11m1 :: t11,t12,t13,t21,t22,t23,t31,t32,t33
real(SP) dt
end subroutine evolve_k5

END INTERFACE
!...........................................................................

!compute D^a(D_b(alpha)) terms

     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)

!raise the first index of t{ab}
 
     call evolve_k2(gup11,gup12,gup13,gup22,gup23,gup33,          &
                    t11,t12,t13,t21,t22,t23,t31,t32,t33) 

!compute Lie_beta(K^a_b) terms

     call evolve_k3a(k11,k12,k13,k21,k22,k23,k31,k32,k33,          &
                    beta1,beta2,beta3,                            &
                    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,beta3,                            &
                    t11,t12,t13,t12,t22,t23,t31,t32,t33,          &
                    nx,ny,nz,h)

!compute alpha*(R^a_b + K*K^a_b) terms

     call 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)

! evolve one timestep via leapfrog method ( i.e. time-centered)

     call evolve_k5(k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32m1,k33m1, &
                    t11,t12,t13,t21,t22,t23,t31,t32,t33,dt)


return
end






