
      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


program t2

!***********************************************************************
!     This code evolves binary black hole spacetimes using the ADM 3+1
!  formalism.  A initial spacelike hypersurface is formed by using
!  an exact form for a known spacetime, or by solving the time-independent
!  Einstein equations, also known as the constraint equations.  This
!  initial hypersurface is evolved using the time-dependent Einstein
!  equations.  In this manner, a foliation of the full spacetime can
!  be computed.
!     The hypersurfaces are freely evolved, i.e.--the constraint
!  equations are not used in the solution of hypersurfaces subsequent
!  to the initial one.  The evolution is a three-level update scheme,
!  using time-centered and space-centered derivatives.  At the boundaries
!  of the computational domain, the outer edges of the grid and the
!  boundary of the excised singularity, we use one-sided differences
!  to compute spatial derivatives.
!     Two important subroutines used are the apparent
!  horizon solver and a generic interpolator/extrapolator.  The program
!  is designed in a modular fasion to allow a variety of finite differencing
!  schemes and conforms to a DAGH compatible data structure.
!
!     Randall R. Correll
!     Mijan F. Huq
!     Scott A. Klasky
!     Richard A. Matzner
!
!     COPYRIGHT 1995
!
! Translation to HPF: Tomasz Haupt (1996)             haupt@npac.syr.edu
!***********************************************************************


use types
use hpfprocessors

integer(I4B), parameter       :: nx=33,ny=33,nz=33
real(SP)                      :: xmin,ymin,zmin
real(SP)                      :: h,amp,w,dt
real(SP), dimension(nx,ny,nz) :: g11,g12,g13,g22,g23,g33
real(SP), dimension(nx,ny,nz) :: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
real(SP), dimension(nx,ny,nz) :: k11,k12,k13,k21,k22,k23,k31,k32,k33
real(SP), dimension(nx,ny,nz) :: k11m1,k12m1,k13m1,k21m1,k22m1,k23m1, &
                                 k31m1,k32m1,k33m1
real(SP), dimension(nx,ny,nz) :: 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) :: nt = 20
real(SP)     :: time

!system clock
integer count,count_rate,count_max
real clock

!.......................................................................
INTERFACE
subroutine  ivp(g11,g12,g13,g22,g23,g33, &
            k11,k12,k13,k21,k22,k23,k31,k32,k33, &
            nx,ny,nz,xmin,ymin,zmin,h,amp,w)   
use types
use hpfprocessors
real(SP), dimension(:,:,:) :: g11,g12,g13,g22,g23,g33
real(SP), dimension(:,:,:) :: k11,k12,k13,k21,k22,k23,k31,k32,k33
integer(I4B) :: nx,ny,nz 
real(SP) :: xmin,ymin,zmin,h
real(SP) :: amp, w

!HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
!HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
!HPF$ ALIGN WITH *g11 :: k11,k12,k13,k21,k22,k23,k31,k32,k33

end subroutine ivp

subroutine setgauge(alpha,beta1,beta2,beta3)
use hpfprocessors
use types
real(SP), dimension(:,:,:) :: alpha,beta1,beta2,beta3
!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: alpha
!HPF$ ALIGN WITH *alpha                  :: beta1, beta2, beta3
end subroutine setgauge


subroutine initialize(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,dt,time)

use types
use hpfprocessors

real(SP), dimension(:,:,:),intent(INOUT):: g11,g12,g13,g22,g23,g33
real(SP), dimension(:,:,:),intent(INOUT):: g11m1,g12m1,g13m1,g22m1,g23m1,g33m1
real(SP), dimension(:,:,:),intent(INOUT):: 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)   :: 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
real(SP) h,dt,time
end subroutine initialize

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

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

write(6,*) ' --- start execution --- '

call system_clock(count,count_rate,count_max)
if(count_rate == 0) then
   write(6,*) 'Oops, no system clock'
   write(6,*) 'STOP'
   stop
else
   write(6,*) 'system clock'
   write(6,*) 'count,rate,max:',count,count_rate,count_max
   clock=real(count)/count_rate
   write(6,*) 'start time:',clock
endif

h=0.25
xmin=-3.9375
ymin=-3.9375
zmin=-3.9375
amp=1.0d-05
w=1.0d0
time=0.0d0
dt=h*0.25

! initialize t=-1 data

call ivp(g11,g12,g13,g22,g23,g33, &
            k11,k12,k13,k21,k22,k23,k31,k32,k33, &
            nx,ny,nz,xmin,ymin,zmin,h,amp,w)

call setgauge(alpha,beta1,beta2,beta3)

call initialize(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,dt,time)

!reinitialize t=0 data

call ivp(g11,g12,g13,g22,g23,g33, &
            k11,k12,k13,k21,k22,k23,k31,k32,k33, &
            nx,ny,nz,xmin,ymin,zmin,h,amp,w)


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

write(6,*) '===== END OF JOB ===='
call system_clock(count,count_rate,count_max)
write(6,*) 'system clock'
write(6,*) 'count,rate,max:',count,count_rate,count_max
clock= real(count)/count_rate - clock
write(6,*) 'end time:',clock, ' seconds'

!at hoc verification
write(6,*) 'g11,g12,g13; g22,g23,g33'
write(6,*) g11(2,3,4),g12(3,4,5),g13(4,5,6)
write(6,*) g22(4,2,3),g23(5,4,3),g33(6,5,4)

stop 
end







