
      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







module evol_interface

!.........................................................................
INTERFACE
subroutine tensors(g11,g12,g13,g22,g23,g33, &
                   gup11,gup12,gup13,gup22,gup23,gup33, &
                   c111,c112,c113,c122,c123,c133, & 
                   c211,c212,c213,c222,c223,c233, & 
                   c311,c312,c313,c322,c323,c333, &
                   r11,r12,r13,r22,r23,r33, &
                   nx,ny,nz,h)
use hpfprocessors
use types

real(SP), dimension(:,:,:),intent(IN)  :: g11,g12,g13,g22,g23,g33
real(SP), dimension(:,:,:),intent(OUT) :: gup11,gup12,gup13,gup22,gup23,gup33
real(SP), dimension(:,:,:),intent(OUT) :: c111,c112,c113,c122,c123,c133
real(SP), dimension(:,:,:),intent(OUT) :: c211,c212,c213,c222,c223,c233
real(SP), dimension(:,:,:),intent(OUT) :: c311,c312,c313,c322,c323,c333
real(SP), dimension(:,:,:),intent(OUT) :: r11,r12,r13,r22,r23,r33

!HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
!HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
!HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
!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 :: r11,r12,r13,r22,r23,r33

integer(I4B)  nx, ny, nz
real(SP) h
end subroutine tensors


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
end subroutine evolve_g

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
end subroutine evolve_k

subroutine swap(a,b)
use types
use hpfprocessors

real(SP), dimension(:,:,:), intent(INOUT) :: a,b
!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC  :: a
!HPF$ ALIGN WITH *a :: b
end subroutine swap

      subroutine constraint(nx,ny,nz,                  &
                k11,k12,k13,k21,k22,k23,k31,k32,k33,   &
                c111,c112,c113,c122,c123,c133,         &
                c211,c212,c213,c222,c223,c233,         &
                c311,c312,c313,c322,c323,c333,         &
                r11,r12,r13,r22,r23,r33,               &
                gup11,gup12,gup13,gup22,gup23,gup33,   &
                hamnorm,mom1norm,mom2norm,mom3norm,h)

use types
use hpfprocessors
use diff_interface

integer(I4B)  nx, ny, nz

real(SP), dimension(:,:,:), intent(IN)  ::  k11,k12,k13
real(SP), dimension(:,:,:), intent(IN)  ::  k21,k22,k23
real(SP), dimension(:,:,:), intent(IN)  ::  k31,k32,k33

real(SP), dimension(:,:,:), intent(IN)  ::  r11, r12, r13
real(SP), dimension(:,:,:), intent(IN)  ::       r22, r23
real(SP), dimension(:,:,:), intent(IN)  ::            r33

real(SP), dimension(:,:,:), intent(IN)  ::  gup11,gup12,gup13
real(SP), dimension(:,:,:), intent(IN)  ::        gup22,gup23
real(SP), dimension(:,:,:), intent(IN)  ::              gup33

!     connection coefficients:  C^i_jk
real(SP), dimension(:,:,:), intent(IN)  ::  c111, c112, c113
real(SP), dimension(:,:,:), intent(IN)  ::        c122, c123
real(SP), dimension(:,:,:), intent(IN)  ::              c133
real(SP), dimension(:,:,:), intent(IN)  ::  c211, c212, c213
real(SP), dimension(:,:,:), intent(IN)  ::        c222, c223
real(SP), dimension(:,:,:), intent(IN)  ::              c233
real(SP), dimension(:,:,:), intent(IN)  ::  c311, c312, c313
real(SP), dimension(:,:,:), intent(IN)  ::        c322, c323
real(SP), dimension(:,:,:), intent(IN)  ::              c333

!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 :: gup11,gup12,gup13,gup22,gup23,gup33
!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  

real(SP) ::  hamnorm,mom1norm,mom2norm,mom3norm,h 
end subroutine constraint

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

end module evol_interface




subroutine evolve(g11,g12,g13,g22,g23,g33,                         &
            g11m1,g12m1,g13m1,g22m1,g23m1,g33m1,                   &
            k11,k12,k13,k21,k22,k23,k31,k32,k33,                   &
            k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32m1,k33m1, &
            alpha,beta1,beta2,beta3,                               &
            nx,ny,nz,h,time,dt,nt)

use types
use hpfprocessors
use evol_interface 

real(SP), dimension(:,:,:) :: g11,g12,g13,g22,g23,g33
real(SP), dimension(:,:,:) :: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
real(SP), dimension(:,:,:) :: k11,k12,k13,k21,k22,k23,k31,k32,k33
real(SP), dimension(:,:,:) :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1, &
                                 k31m1,k32m1,k33m1
real(SP), dimension(:,:,:) :: alpha,beta1,beta2,beta3

!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 :: k11,k12,k13,k21,k22,k23,k31,k32,k33
!HPF$ ALIGN WITH *g11 :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1,k31m1,k32m1,k33m1
!HPF$ ALIGN WITH *g11 :: alpha,beta1,beta2,beta3

integer(I4B) nx,ny,nz,nt
real(SP) h,time,dt
real(SP) :: tol=100.0

real(SP), dimension(size(g11,1),size(g11,2),size(g11,3)) ::  &
                        gup11,gup12,gup13,gup22,gup23,gup33, &
                        c111,c112,c113,c122,c123,c133,       & 
                        c211,c212,c213,c222,c223,c233,       & 
                        c311,c312,c313,c322,c323,c333,       &
                        r11,r12,r13,r22,r23,r33,             &
                        t7,t8,t9,temp1

!HPF$ ALIGN WITH g11 :: gup11,gup12,gup13,gup22,gup23,gup33  
!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 :: r11,r12,r13,r22,r23,r33
!HPF$ ALIGN WITH g11 :: t7,t8,t9,temp1

real(SP) ::  hamnorm,mom1norm,mom2norm,mom3norm 

! the evolution sequence is:
! evolve * swap * tensors
! while tensors = inverse * connect * ricci is input for evolve
! therefore tensors must be called once before the evolution loop

call tensors(g11,g12,g13,g22,g23,g33,          &
             gup11,gup12,gup13,gup22,gup23,gup33,                   &
             c111,c112,c113,c122,c123,c133,                         & 
             c211,c212,c213,c222,c223,c233,                         & 
             c311,c312,c313,c322,c323,c333,                         &
             r11,r12,r13,r22,r23,r33,                               &
             nx,ny,nz,h)

! constraints monitors

call  constraint(nx,ny,nz,                             &
                k11,k12,k13,k21,k22,k23,k31,k32,k33,   &
                c111,c112,c113,c122,c123,c133,         &
                c211,c212,c213,c222,c223,c233,         &
                c311,c312,c313,c322,c323,c333,         &
                r11,r12,r13,r22,r23,r33,               &
                gup11,gup12,gup13,gup22,gup23,gup33,   &
                hamnorm,mom1norm,mom2norm,mom3norm,h) 

it=0
write(6,*) 'Iter   Time    HamErr    Mom1Err    Mom2Err    Mom3Err'
write(6,'(i5,f9.5,1P,4e12.3)') it,time, hamnorm,mom1norm,mom2norm,mom3norm
 
!---------------------- EVOLUTION LOOP ------------------------------------

do it=1,nt
 time=time+dt

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

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


!swap time levels

	call swap(g11,g11m1)
	call swap(g12,g12m1)
	call swap(g13,g13m1)
	call swap(g22,g22m1)
	call swap(g23,g23m1)
	call swap(g33,g33m1)

	call swap(k11,k11m1)
	call swap(k12,k12m1)
	call swap(k13,k13m1)
	call swap(k21,k21m1)
	call swap(k22,k22m1)
	call swap(k23,k23m1)
	call swap(k31,k31m1)
	call swap(k32,k32m1)
	call swap(k33,k33m1)

 call tensors(g11,g12,g13,g22,g23,g33,                              &
             gup11,gup12,gup13,gup22,gup23,gup33,                   &
             c111,c112,c113,c122,c123,c133,                         & 
             c211,c212,c213,c222,c223,c233,                         & 
             c311,c312,c313,c322,c323,c333,                         &
             r11,r12,r13,r22,r23,r33,                               &
             nx,ny,nz,h)

 call  constraint(nx,ny,nz,                            &
                k11,k12,k13,k21,k22,k23,k31,k32,k33,   &
                c111,c112,c113,c122,c123,c133,         &
                c211,c212,c213,c222,c223,c233,         &
                c311,c312,c313,c322,c323,c333,         &
                r11,r12,r13,r22,r23,r33,               &
                gup11,gup12,gup13,gup22,gup23,gup33,   &
                hamnorm,mom1norm,mom2norm,mom3norm,h) 


 write(6,'(i5,f9.5,1P,4e12.3)') it,time, hamnorm,mom1norm,mom2norm,mom3norm

 if(hamnorm.gt.tol) exit
enddo

return
end





