
      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  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)   
!----------------------------------------------------------------------------
! Mijan F. Huq
! Center for Relativity,
! University of Texas at Austin
!
! Program to generate linearized gravitional wave initial data.
! References: 
!   Saul A. Teukolsky, Phys Rev D, vol 26, no 4, 1982
!   Kenneth Eppley, "Pure Gravitational Waves", in Sources...
!
! Comments:
!  Did not fix with L'Hopital's rule terms with x^2 + y^2 in the denominator
!  instead chose to off-center grid from origin of coordinate system...
!
! HPF Implementation by T. Haupt  NPAC at Syracuse University 
!                 (haupt@npac.syr,edu)
!----------------------------------------------------------------------------

use hpfprocessors
use types

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

real(SP),dimension(size(g11,1),size(g11,2),size(g11,3)) :: xx,yy,zz,r,A_,B_,C_

!HPF$ ALIGN WITH g11 :: xx,yy,zz,r,A_,B_,C_

!.......................................................................
INTERFACE
subroutine setradial(r,A_,B_,C_,amp, w)
use hpfprocessors
use types
real(SP), dimension(:,:,:) ::  r, A_, B_, C_
real(SP)  amp,w
!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: r
!HPF$ ALIGN WITH *r                      :: A_, B_, C_
end subroutine setradial

subroutine metcomp(g11,g12,g13,g22,g23,g33,A_,B_,C_,r,xx,yy,zz) 
use hpfprocessors
use types
real(SP),dimension(:,:,:)  :: g11,g12,g13,g22,g23,g33
real(SP),dimension(:,:,:)  :: A_, B_, C_, r, xx, yy, zz
!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: g11
!HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
!HPF$ ALIGN WITH *g11 :: A_, B_, C_, r, xx, yy, zz
end subroutine metcomp
END INTERFACE
!.......................................................................


!do i=1,nx
!  xx(i,:,:)=xmin+(i-1)*h
!enddo

!do i=1,ny
!  yy(:,i,:)=ymin+(i-1)*h
!enddo

!do i=1,nz
!  zz(:,:,i)=zmin+(i-1)*h
!enddo

forall(i=1:nx,j=1:ny,k=1:nz) xx(i,j,k)=xmin+(i-1)*h
forall(i=1:nx,j=1:ny,k=1:nz) yy(i,j,k)=ymin+(j-1)*h
forall(i=1:nx,j=1:ny,k=1:nz) zz(i,j,k)=zmin+(k-1)*h

r = sqrt( xx*xx + yy*yy + zz*zz )

call setradial(r,A_,B_,C_,amp, w)
call metcomp(g11,g12,g13,g22,g23,g33,A_,B_,C_,r,xx,yy,zz)

write(6,*) g11(2,3,4)
write(6,*) g12(2,3,4)
write(6,*) g13(2,3,4)
write(6,*) g22(2,3,4)
write(6,*) g23(2,3,4)
write(6,*) g33(2,3,4)

k11 = 0.0d0
k12 = 0.0d0
k13 = 0.0d0
k21 = 0.0d0
k22 = 0.0d0
k23 = 0.0d0
k31 = 0.0d0
k32 = 0.0d0
k33 = 0.0d0

return
end

!---------------------------------------------------------------------------
!    Routine to set up Radial functions
!
!   Form time symmetric solution by superposing f(t+r) - f(t-r)
!----------------------------------------------------------------------------

subroutine setradial(r,A_,B_,C_,amp, w)
!---------------------------------------------------------------------------
use types
use hpfprocessors

real(SP), dimension(:,:,:) ::  r, A_, B_, C_
real(SP)  amp,w

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: r
!HPF$ ALIGN WITH *r                      :: A_, B_, C_

real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: fo, fi, fo1, fi1, fo2
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: fi2, fo3, fi3, fo4, fi4 

!HPF$ ALIGN WITH r :: fo, fi, fo1, fi1, fo2, fi2, fo3, fi3, fo4, fi4

! Set up generating function and its derivatives for ingoing and outgoing

  fi=-2*amp*r/w**2 * exp(-r**2/w**2)
  fo= 2*amp*r/w**2 * exp(-r**2/w**2)

  fi1=2*amp*exp(-(r**2/w**2))*(-w**2+2*r**2)/w**4  
  fo1=2*amp*exp(-(r**2/w**2))*(-w**2+2*r**2)/w**4  

  fi2=-4*amp*r*exp(-r**2)*(-3*w**2+2*r**2)/w**6  
  fo2= 4*amp*r*exp(-r**2)*(-3*w**2+2*r**2)/w**6  

  fi3=4*amp*exp(-r**2/w**2)*(3*w**4-12*r**2*w**2+4*r**4)/w**8
  fo3=4*amp*exp(-r**2/w**2)*(3*w**4-12*r**2*w**2+4*r**4)/w**8

  fi4=-8*amp*r*exp(-r**2/w**2)*(15*w**4-20*r**2*w**2+4*r**4)/w**10
  fo4= 8*amp*r*exp(-r**2/w**2)*(15*w**4-20*r**2*w**2+4*r**4)/w**10

      A_ = 3.0d0 * ( ( fo2 - fi2 )/r**3 +     & 
              ( fo1 + fi1 ) * 3.0d0 / r**4 +  &
              ( fo  - fi  ) * 3.0d0 / r**5)

      B_ = - 1.0d0 /r**2 *(                         &
                 (fo3 + fi3) + 3.0d0/r*(fo2-fi2) +  &
                 6.0d0/r**2 *(fo1 + fi1) +          &
                 6.0d0/r**3 *(fo - fi) )

      C_ = 0.25/r *(                                       &
                 ( fo4 - fi4) + 2.0d0/r *(fo3 + fi3)  +    &
                 9.0d0/r**2 *(fo2 - fi2) +                 &
                 21.0d0/r**3 *(fo1 + fi1) +                &
                 21.0d0/r**4 *(fo - fi) )                  
      return
      end
!-----------------------------------------------------------------------
subroutine metcomp(g11,g12,g13,g22,g23,g33,A_,B_,C_,r,xx,yy,zz) 
!-----------------------------------------------------------------------
use hpfprocessors
use types
implicit none
real(SP),dimension(:,:,:)  :: g11,g12,g13,g22,g23,g33
real(SP),dimension(:,:,:)  :: A_, B_, C_, r, xx, yy, zz

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: g11
!HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
!HPF$ ALIGN WITH *g11 :: A_, B_, C_, r, xx, yy, zz

real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: frr, frth, frph, fthph 
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: d1fthth,d1fphph,d2fthth
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: d2fphph,gthph,gphph
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: grr,grth,grph,gthth
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: Lrx,Lry,Lrz,Ltx,Lty
real(SP),dimension(size(r,1),size(r,2),size(r,3)) :: Ltz,Lpx,Lpy,Lpz 

!HPF$ ALIGN WITH g11 :: frr, frth, frph, fthph, d1fthth,d1fphph,d2fthth
!HPF$ ALIGN WITH g11 :: d2fphph,gthph,gphph,grr,grth,grph,gthth
!HPF$ ALIGN WITH g11 :: Lrx,Lry,Lrz,Ltx,Lty,Ltz,Lpx,Lpy,Lpz 

!........................................................................
INTERFACE
subroutine angfunc(m,x, y, z, r,                 &
               frr, frth, frph, fthph, d1fthth,  &
               d1fphph,d2fthth,d2fphph)
use hpfprocessors
use types
integer(I4B)   m

real(SP), dimension(:,:,:) :: x, y, z, r
real(SP), dimension(:,:,:) :: frr, frth, frph, fthph
real(SP), dimension(:,:,:) :: d1fthth,d1fphph,d2fthth,d2fphph

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: r
!HPF$ ALIGN WITH *r :: x, y, z, frr, frth, frph, fthph
!HPF$ ALIGN WITH *r ::  d1fthth,d1fphph,d2fthth,d2fphph
end subroutine angfunc
END INTERFACE
!........................................................................

      integer(I4B)   mm

      
!      want m=0 waves

      mm = 0

!     Get values of angular functions

      call angfunc(0,xx, yy, zz, r, frr, frth, frph, fthph,     &
                  d1fthth, d1fphph,d2fthth,d2fphph)                 

!         Form spherical coordinate metric components
      grr = 1.0d0 + A_ * frr

      grth = B_ *frth * r
      grph = B_ *frph * sqrt(xx**2 + yy**2)
      gthth=( 1.0d0 + C_*d1fthth + A_ * d2fthth )*r**2
      gthph=((A_-2.0d0*C_)*fthph)*r* sqrt(xx**2 + yy**2)

      gphph=(1.0d0 + C_*d1fphph + A_*d2fphph)*(xx**2 + yy**2) 

!         Form transformation matrix

      Lrx = xx / r
      Lry = yy / r
      Lrz = zz / r

      Ltx = xx * zz / r**2 /sqrt(xx**2 + yy**2)
      Lty = yy * zz / r**2 /sqrt(xx**2 + yy**2)
      Ltz = -sqrt(xx**2 + yy**2) / r**2

      Lpx = -yy / (xx**2 + yy**2)
      Lpy = xx / (xx**2 + yy**2)
      Lpz = 0.0d0



!         Form Cartesian metric components
      g11 = grr * Lrx**2 + 2.0d0*grth*Lrx*Ltx +         &
           2.0d0*grph * Lrx * Lpx + gthth * Ltx**2 +    &
           2.0d0 * gthph * Ltx * Lpx + gphph * Lpx**2   

      g12 = grr * Lrx*Lry + 2.0d0*grth*Lrx*Lty +        &
           2.0d0*grph * Lrx * Lpy + gthth * Ltx*Lty +   &
           2.0d0 * gthph * Ltx * Lpy + gphph * Lpx*Lpy  

      g13 = grr * Lrx*Lrz + 2.0d0*grth*Lrx*Ltz +        &
           2.0d0*grph * Lrx * Lpz + gthth * Ltx*Ltz +   &
           2.0d0 * gthph * Ltx * Lpz + gphph * Lpx*Lpz  


      g22 = grr * Lry**2 + 2.0d0*grth*Lry*Lty +         &
           2.0d0*grph * Lry * Lpy + gthth * Lty**2 +    &
           2.0d0 * gthph * Lty * Lpy + gphph * Lpy**2   

      g23 = grr * Lry*Lrz + 2.0d0*grth*Lry*Ltz +        &
           2.0d0*grph * Lry * Lpz + gthth * Lty*Ltz +   &
           2.0d0 * gthph * Lty * Lpz + gphph * Lpy*Lpz  

      g33 = grr * Lrz**2 + 2.0d0*grth*Lrz*Ltz +         &
           2.0d0*grph * Lrz * Lpz + gthth * Ltz**2 +    &
           2.0d0 * gthph * Ltz * Lpz + gphph * Lpz**2   
		
      return
      end
!-----------------------------------------------------------------------
! Routine to calculate angular functions

subroutine angfunc(m,x, y, z, r,                        &
                   frr, frth, frph, fthph, d1fthth,     &
                   d1fphph,d2fthth,d2fphph)
!-----------------------------------------------------------------------
use hpfprocessors
use types
implicit none

integer(I4B)   m

real(SP), dimension(:,:,:) :: x, y, z, r
real(SP), dimension(:,:,:) :: frr, frth, frph, fthph
real(SP), dimension(:,:,:) :: d1fthth,d1fphph,d2fthth,d2fphph

!HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: r
!HPF$ ALIGN WITH *r :: x, y, z, frr, frth, frph, fthph
!HPF$ ALIGN WITH *r ::d1fthth,d1fphph,d2fthth,d2fphph  

real(SP), dimension(size(r,1),size(r,2),size(r,3)) :: sint, cost, sinp, cosp

!HPF$ ALIGN WITH r :: sint, cost, sinp, cosp

      sint = sqrt((x**2 + y**2)/(x**2 + y**2 + z**2))
      cost = z / sqrt(x**2 + y**2 + z**2)
         
      sinp = y/sqrt(x**2 + y**2)
      cosp = x/sqrt(x**2 + y**2)


      frr = 2.0d0 - 3.0d0*sint**2
      frth= -3.0d0 * sint * cost
      frph = 0.0d0
      fthph= 0.0d0

      d1fthth= 3.0d0* sint**2
      d1fphph= -d1fthth

      d2fthth=-1.0d0
      d2fphph= 3.0d0*sint**2 - 1.0d0

      return
      end
!-----------------------------------------------------------------------













