      program onebh
c*********************************************************************
c NCSA Copyright Statement
c Copyright (c) 1993 and 1994 by University of Illinois Board of Trustees 
c
c This code was written by Peter Anninos, David Bernstein, David Hobill,
c Edward Seidel, Larry Smarr, and John Towns at NCSA, with support from
c the University of Illinois, the State of Illinois, the National Science
c Foundation, and other federal agencies.
c
c We encourage the use of this code for studying black hole spacetimes,
c as long as credit and notification are provided.  Please
c retain the copyright notice and include this statement
c with any copies that you make, and send notification of the use of 
c this code to mosaic@jean-luc.ncsa.uiuc.edu.  We would appreciate receiving
c any copies of work done with this code.  Also, if you make any useful
c modifications or extensions to this code, please let us know and send us
c a copy.  Thanks very much, and enjoy!
c**********************************************************************
c
c  Main driver to integrate the metric and momentum variables of the
c  Einstein equations in the 3+1 form and to control I/O operations.
c
c  This code solves the single, axisymmetric black hole spacetime
c  which is distorted by the presence of an initially time symmetric
c  gravitational wave (Brill wave).
c
c  Also computed are the Zerilli function, Newman-Penrose scalars,
c  Riemann invariants, Bel-Robinson vectors, apparent horizons
c  (location, mass and embedding) and radiative losses.
c
c  Leap frog with 1/2 step extrapolation is used for the evolution,
c  and centered finite difference formulas are used for spatial
c  derivatives for 2nd order accuracy in space and time. The evolved
c  variables are the metric components (g11,g22,g33) and the extrinsic 
c  curvature components (h11,h22,h33,h12). The lapse (g00) and 
c  nonvanishing shift components (g01,g02) are computed every time
c  slice. The conformal factor (psi, and phi=4*ln(psi)) and all its
c  derivatives (psix,psiy,psixx,psiyy,psixy) are computed on the
c  initial time slice and do not evolve.
c
c  Descriptions of other key variables and documentation of output
c  files are included in the radhor subroutine which computes and
c  prints out the radiation indicators and apparent horizon (location,
c  embedding and gaussian curvature).
c
c  There is a scale parameter 'm' which is the mass of the unperturbed
c  wormhole. It is taken to be 2 in this code.
c
c**********************************************************************
c
      include 'param.h'
      common /admmablk/ adm,igdet

      real coty(im,jm)
      real psi(im,jm),psix(im,jm),psixx(im,jm),phi(im,jm)
      real psixy(im,jm),psiy(im,jm),psiyy(im,jm)

      real g11(im,jm),g22(im,jm),g33(im,jm)
      real g11p(im,jm),g22p(im,jm),g33p(im,jm)
      real g11o(jm),g22o(jm),g33o(jm)
      real g11ex(im,jm),g22ex(im,jm),g33ex(im,jm)

      real h11(im,jm),h22(im,jm),h12(im,jm),h33(im,jm)
      real h11p(im,jm),h22p(im,jm),h12p(im,jm),h33p(im,jm)
      real h11o(jm),h22o(jm),h12o(jm),h33o(jm)
      real h11ex(im,jm),h22ex(im,jm),h12ex(im,jm),h33ex(im,jm)

      real g00(im,jm),g00p(im,jm),g00ppp(im,jm)
      real g00pp(im,jm),g00ex(im,jm)

      real g01(im,jm),g02(im,jm),chi(im,jm)
      real g01p(im,jm),g02p(im,jm)
      real g01ex(im,jm),g02ex(im,jm)

      real ric11(im,jm),ric22(im,jm),ric12(im,jm),ric33(im,jm)
      real scurv(im,jm)
      character*10 var,file
      integer pt

c  Setup the output data files
        open(unit=37,file='info.dat')

        open(unit=46,file='psi4.xdat')
        open(unit=47,file='psi0.xdat')
        open(unit=48,file='psi4d.xdat')
        open(unit=49,file='psi0d.xdat')
        open(unit=50,file='bre.tdat')
        open(unit=51,file='brt.tdat')
        open(unit=52,file='brn.tdat')
        open(unit=54,file='crvi.tdat')
        open(unit=55,file='crvj.tdat')
        open(unit=57,file='y13.tdat')
        open(unit=58,file='y23.tdat')
        open(unit=59,file='psi0.tdat')
        open(unit=61,file='psi2.tdat')
        open(unit=63,file='psi4.tdat')

        open(unit=64,file='ahm.dat')
        open(unit=70,file='embd.dat')
        open(unit=71,file='cpoe.dat')
        open(unit=72,file='hor.dat')

        open(unit=90,file='zpsi2.dat')
        open(unit=91,file='zpsi4.dat')

        open(unit=82,file='tezpsi2.dat')
        open(unit=83,file='tezpsi4.dat')
        open(unit=84,file='tezpsit.dat')
        open(unit=86,file='tebre.dat')
        open(unit=88,file='tepsi4.dat')

        open(unit=92,file='dedtzpsi2.dat')
        open(unit=93,file='dedtzpsi4.dat')
        open(unit=94,file='dedtzpsit.dat')
        open(unit=96,file='dedtbre.dat')
        open(unit=98,file='dedtpsi4.dat')

        write(37,*)' Total time for the run is: ',tmax,' M'
        write(37,*)' Grid resolution: ',imax,' x ',jmax
        write(37,*)'   delta eta    : ',delx
        write(37,*)'   delta theta  : ',dely
        write(37,*)'   delta t      : ',delt
        write(6,*)' Total time for the run is: ',tmax,' M'
        write(6,*)' Grid resolution: ',imax,' x ',jmax
        write(6,*)'   delta eta    : ',delx
        write(6,*)'   delta theta  : ',dely
        write(6,*)'   delta t      : ',delt
        write(6,*)' 3D data output every ',prnt3d,' M'
 
c  Set up the initial data and compute the conformal factor
      call initdat (g11,g22,g33,g01,g02,
     &g11p,g22p,g33p,g01p,g02p,g11o,g22o,g33o,h11,h22,h12,h33,
     &h11p,h22p,h12p,h33p,h11o,h12o,h22o,h33o,
     &psi,psix,psiy,psixx,psiyy,psixy,coty,ric11,ric22,ric12,ric33,phi)

      
        it=1
        pt=1
        nt=0
        time=0.

c  determine the lapse depending on the chosen slicing condition
c     ich = 2 is for the analytic Schwartzchild lapse for static black
c             holes
c     otherwise choose alpha=1. everywhere initially (maximal slice
c             solution symmetric accross boundary)
        if (ich.eq.2)then
          do 286 j=1,jmax+1
          do 285 i=1,imax+3 
            xx=(float(i)-2.)*delx
            htan=tanh(0.5*xx)
            g00(i,j)=htan
            g00p(i,j)=htan
            g00pp(i,j)=htan
            g00ppp(i,j)=htan
285       continue
286       continue
        else
          do 291 j=1,jmax+1
          do 290 i=1,imax+3
            g00(i,j)=1.
            g00p(i,j)=1.
            g00pp(i,j)=1.
            g00ppp(i,j)=1.
290       continue
291       continue
        end if

c  Write out the HDF arrays on the initial time slice.
      i3dprint=1
      if(i3dprint0.eq.0)i3dprint=0
      if(i3dprint.eq.1)then
	write(6,*)' printing metric data at t = ',time
        var=' g00 '
        file= 'g00hdf'
        call writehdf(var,file,g00,time,im,jm)
        var=' g11 '
        file= 'g11hdf'
        call writehdf(var,file,g11,time,im,jm)
        var=' g22 '
        file= 'g22hdf'
        call writehdf(var,file,g22,time,im,jm)
        var=' g33 '
        file= 'g33hdf'
        call writehdf(var,file,g33,time,im,jm)
        var=' h11 '
        file= 'h11hdf'
        call writehdf(var,file,h11,time,im,jm)
        var=' h12 '
        file= 'h12hdf'
        call writehdf(var,file,h12,time,im,jm)
        var=' h22 '
        file= 'h22hdf'
        call writehdf(var,file,h22,time,im,jm)
        var=' h33 '
        file= 'h33hdf'
        call writehdf(var,file,h33,time,im,jm)
        var=' g01 '
        file= 'g01hdf'
        call writehdf(var,file,g01,time,im,jm)
        var=' g02 '
        file= 'g02hdf'
        call writehdf(var,file,g02,time,im,jm)
      endif

c  compute apparent horizon and radiation variables on first slice
      call radhor(nt,g11,g22,g33,h11,h22,h33,h12,psi,phi,g00,
     $            g01,g02,ric11,ric22,ric12,ric33,i3dprint)
      call zerilli(psi,g11,g22,g33,time,psix)
      i3dprint=0
      nmax=adm*tmax/delt+1


c
c  The main loop to evolve the Einstein equations
c=======================================================================
c
      do 1000 nt=1,nmax

c  the metric component evolution
c-----------------------------------------------------------------------
        do 301 j=2,jmax
        do 300 i=2,imax+2
          psin=1./psi(i,j)
          pxpin=psix(i,j)*psin
          pypin=psiy(i,j)*psin
          a=g11ex(i,j)
          b=g22ex(i,j)
          d=g33ex(i,j)
          v1=g01(i,j)
          v2=g02(i,j)
 
          v1x=.5*(g01(i+1,j)-g01(i-1,j))/delx
          v1y=.5*(g01(i,j+1)-g01(i,j-1))/dely
          v2x=.5*(g02(i+1,j)-g02(i-1,j))/delx
          v2y=.5*(g02(i,j+1)-g02(i,j-1))/dely
          ax=.5*(g11ex(i+1,j)-g11ex(i-1,j))/delx
          bx=.5*(g22ex(i+1,j)-g22ex(i-1,j))/delx
          dx=.5*(g33ex(i+1,j)-g33ex(i-1,j))/delx
          ay=.5*(g11ex(i,j+1)-g11ex(i,j-1))/dely
          by=.5*(g22ex(i,j+1)-g22ex(i,j-1))/dely
          dy=.5*(g33ex(i,j+1)-g33ex(i,j-1))/dely
 
          g11p(i,j)=g11(i,j)
          g22p(i,j)=g22(i,j)
          g33p(i,j)=g33(i,j)

          s1 = 4*a*pypin*v2
          s2 = ay*v2
          s3 = 2*a*v1x
          s4 = 4*a*pxpin*v1
          s5 = ax*v1
          s6 = -2*h11(i,j)*g00(i,j)
          g11(i,j)=g11(i,j)+delt*(s1+s2+s3+s4+s5+s6)

          s1 = 2*b*v2y
          s2 = 4*b*pypin*v2
          s3 = by*v2
          s4 = 4*b*pxpin*v1
          s5 = bx*v1
          s6 = -2*h22(i,j)*g00(i,j)
          g22(i,j)=g22(i,j)+delt*(s1+s2+s3+s4+s5+s6)

          s1 = 4*d*pypin*v2
          s2 = dy*v2
          s3 = 2*coty(i,j)*d*v2
          s4 = 4*d*pxpin*v1
          s5 = dx*v1
          s6 = -2*h33(i,j)*g00(i,j)
          g33(i,j)=g33(i,j)+delt*(s1+s2+s3+s4+s5+s6)
300     continue
301     continue

c  boundary conditions at z axis (y=0)
        do 310 i=2,imax+2
          g11p(i,1)=g11(i,1)
          g22p(i,1)=g22(i,1)
          g33p(i,1)=g33(i,1)
          g11(i,1)=g11(i,2)
          g22(i,1)=g22(i,2)
          g33(i,1)=g33(i,2)
c  boundary conditions at equator
          g11p(i,jmax+1)=g11(i,jmax+1)
          g22p(i,jmax+1)=g22(i,jmax+1)
          g33p(i,jmax+1)=g33(i,jmax+1)
          g11(i,jmax+1)=g11(i,jmax)
          g22(i,jmax+1)=g22(i,jmax)
          g33(i,jmax+1)=g33(i,jmax)
310     continue

c  boundary conditions at the throat (i=2)
        do 315 j=1,jmax+1
          g11p(1,j)=g11(1,j)
          g22p(1,j)=g22(1,j)
          g33p(1,j)=g33(1,j)
          g11(1,j)=g11(3,j)
          g22(1,j)=g22(3,j)
          g33(1,j)=g33(3,j)
c  boundary conditions at the outer edge of the grid (i=imax+3)
          g11p(imax+3,j)=g11(imax+3,j)
          g22p(imax+3,j)=g22(imax+3,j)
          g33p(imax+3,j)=g33(imax+3,j)
          g11(imax+3,j)=g11o(j)
          g22(imax+3,j)=g22o(j)
          g33(imax+3,j)=g33o(j)
315     continue


c  the Ricci tensor components
c----------------------------------------------------------------------
        do 326 j=2,jmax
        do 325 i=2,imax+2
          d=g33(i,j)
          a=g11(i,j)
          b=g22(i,j)

          ax=.5*(g11(i+1,j)-g11(i-1,j))/delx
          bx=.5*(g22(i+1,j)-g22(i-1,j))/delx
          dx=.5*(g33(i+1,j)-g33(i-1,j))/delx
          ay=.5*(g11(i,j+1)-g11(i,j-1))/dely
          by=.5*(g22(i,j+1)-g22(i,j-1))/dely
          dy=.5*(g33(i,j+1)-g33(i,j-1))/dely
          ayy=(g11(i,j+1)-2.*g11(i,j)+g11(i,j-1))/dydy
          bxx=(g22(i+1,j)-2.*g22(i,j)+g22(i-1,j))/dxdx
          dxx=(g33(i+1,j)-2.*g33(i,j)+g33(i-1,j))/dxdx
          dyy=(g33(i,j+1)-2.*g33(i,j)+g33(i,j-1))/dydy
          deltd2=g33(i+1,j+1)-g33(i-1,j+1)-g33(i+1,j-1)+g33(i-1,j-1)
          dxy=.25*deltd2/dxdy
        
          p=psi(i,j)
          pxx=psixx(i,j)
          pxy=psixy(i,j)
          pyy=psiyy(i,j)
          px=psix(i,j)
          py=psiy(i,j)
          pxpin=px/p
          pypin=py/p
          pxxpin=pxx/p
          pyypin=pyy/p
          pxypin=pxy/p

      s1 = -2*a*pyypin/b
      s2 = -2*a*pypin**2/b
      s3 = -a*dy*pypin/(b*d)
      s4 = -2*a*coty(i,j)*pypin/b
      s5 = a*by*pypin/b**2
      s6 = -2*ay*pypin/b
      s7 = -4*pxxpin
      s8 = 4*pxpin**2
      s9 = -dx*pxpin/d
      s10 = -bx*pxpin/b
      s11 = 2*ax*pxpin/a
      s12 = -ay*dy/(b*d)/4.0
      s13 = -dxx/d/2.0
      s14 = dx**2/d**2/4.0
      s15 = ax*dx/(a*d)/4.0
      s16 = -ay*coty(i,j)/b/2.0
      s17 = ay*by/b**2/4.0
      s18 = -bxx/b/2.0
      s19 = bx**2/b**2/4.0
      s20 = ax*bx/(a*b)/4.0
      s21 = -ayy/b/2.0
      s22 = ay**2/(a*b)/4.0

          r11a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
          r11b=s15+s16+s17+s18+s19+s20+s21+s22
          ric11(i,j)=r11a+r11b

      s1 = 6*pxpin*pypin
      s2 = bx*pypin/b
      s3 = -2*pxypin
      s4 = ay*pxpin/a
      s5 = dx*dy/d**2/4.0
      s6 = bx*dy/(b*d)/4.0
      s7 = -dxy/d/2.0
      s8 = -coty(i,j)*dx/d/2.0
      s9 = ay*dx/(a*d)/4.0
      s10 = bx*coty(i,j)/b/2.0

          r12a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10
          ric12(i,j)=r12a

      s1 = -4*pyypin
      s2 = 4*pypin**2
      s3 = -dy*pypin/d
      s4 = -2*coty(i,j)*pypin
      s5 = 2*by*pypin/b
      s6 = -ay*pypin/a
      s7 = -2*b*pxxpin/a
      s8 = -2*b*pxpin**2/a
      s9 = -b*dx*pxpin/(a*d)
      s10 = -2*bx*pxpin/a
      s11 = ax*b*pxpin/a**2
      s12 = -dyy/d/2.0
      s13 = dy**2/d**2/4.0
      s14 = -coty(i,j)*dy/d
      s15 = by*dy/(b*d)/4.0
      s16 = -bx*dx/(a*d)/4.0
      s17 = by*coty(i,j)/b/2.0
      s18 = ay*by/(a*b)/4.0
      s19 = -bxx/a/2.0
      s20 = bx**2/(a*b)/4.0
      s21 = ax*bx/a**2/4.0
      s22 = -ayy/a/2.0
      s23 = ay**2/a**2/4.0
      s24 = 1

          r22a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
          r22b=s15+s16+s17+s18+s19+s20+s21+s22+s23+s24
          ric22(i,j)=r22a+r22b

      s1 = -2*d*pyypin/b
      s2 = -2*d*pypin**2/b
      s3 = -2*dy*pypin/b
      s4 = -4*coty(i,j)*d*pypin/b
      s5 = by*d*pypin/b**2
      s6 = -ay*d*pypin/(a*b)
      s7 = -2*d*pxxpin/a
      s8 = -2*d*pxpin**2/a
      s9 = -2*dx*pxpin/a
      s10 = -bx*d*pxpin/(a*b)
      s11 = ax*d*pxpin/a**2
      s12 = -dyy/b/2.0
      s13 = dy**2/(b*d)/4.0
      s14 = -coty(i,j)*dy/b
      s15 = by*dy/b**2/4.0
      s16 = -ay*dy/(a*b)/4.0
      s17 = -dxx/a/2.0
      s18 = dx**2/(a*d)/4.0
      s19 = -bx*dx/(a*b)/4.0
      s20 = ax*dx/a**2/4.0
      s21 = by*coty(i,j)*d/b**2/2.0
      s22 = -ay*coty(i,j)*d/(a*b)/2.0
      s23 = d/b

          r33a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
          r33b=s15+s16+s17+s18+s19+s20+s21+s22+s23
          ric33(i,j)=r33a+r33b
325     continue
326     continue

c  boundary conditions at z axis (y=0)
        do 390 i=2,imax+2
          ric11(i,1)=ric11(i,2)
          ric22(i,1)=ric22(i,2)
          ric12(i,1)=-ric12(i,2)
          ric33(i,1)=ric33(i,2)
c  boundary conditions at equator
          ric11(i,jmax+1)=ric11(i,jmax)
          ric22(i,jmax+1)=ric22(i,jmax)
          ric12(i,jmax+1)=-ric12(i,jmax)
          ric33(i,jmax+1)=ric33(i,jmax)
390     continue

c  boundary conditions at the throat (i=2)
        do 395 j=1,jmax+1
          ric11(1,j)=ric11(3,j)
          ric22(1,j)=ric22(3,j)
          ric12(1,j)=-ric12(3,j)
          ric33(1,j)=ric33(3,j)
c  boundary conditions at the outer edge of the grid (i=imax+3)
          ric11(imax+3,j)=ric11(imax+2,j)
          ric22(imax+3,j)=ric22(imax+2,j)
          ric12(imax+3,j)=ric12(imax+2,j)
          ric33(imax+3,j)=ric33(imax+2,j)
395     continue


c  extrapolate to the next half-timestep
c-------------------------------------------------------------------
        do 331 j=1,jmax+1
                do 330 i=1,imax+3
                        h11ex(i,j)=1.5*h11(i,j)-.5*h11p(i,j)
                        h12ex(i,j)=1.5*h12(i,j)-.5*h12p(i,j)
                        h22ex(i,j)=1.5*h22(i,j)-.5*h22p(i,j)
                        h33ex(i,j)=1.5*h33(i,j)-.5*h33p(i,j)
                        g00ex(i,j)=1.5*g00(i,j)-.5*g00p(i,j)
                        g01ex(i,j)=1.5*g01(i,j)-.5*g01p(i,j)
                        g02ex(i,j)=1.5*g02(i,j)-.5*g02p(i,j)
330             continue
331     continue


c  extrinsic curvature evolution
c-------------------------------------------------------------------
        do 351 j=2,jmax
        do 350 i=2,imax+2
          a=g11(i,j)
          b=g22(i,j)
          d=g33(i,j)
          ha=h11ex(i,j)
          hb=h22ex(i,j)
          hc=h12ex(i,j)
          hd=h33ex(i,j)
          v1=g01ex(i,j)
          v2=g02ex(i,j)
          p=psi(i,j)
          p4=p**4
          psin=1./psi(i,j)
          pxpin=psix(i,j)*psin
          pypin=psiy(i,j)*psin

          ax=.5*(g11(i+1,j)-g11(i-1,j))/delx
          bx=.5*(g22(i+1,j)-g22(i-1,j))/delx
          dx=.5*(g33(i+1,j)-g33(i-1,j))/delx
          ay=.5*(g11(i,j+1)-g11(i,j-1))/dely
          by=.5*(g22(i,j+1)-g22(i,j-1))/dely
          dy=.5*(g33(i,j+1)-g33(i,j-1))/dely
          ayy=(g11(i,j+1)-2.*g11(i,j)+g11(i,j-1))/dydy
          bxx=(g22(i+1,j)-2.*g22(i,j)+g22(i-1,j))/dxdx
          dxx=(g33(i+1,j)-2.*g33(i,j)+g33(i-1,j))/dxdx
          dyy=(g33(i,j+1)-2.*g33(i,j)+g33(i,j-1))/dydy
          v1x=.5*(g01ex(i+1,j)-g01ex(i-1,j))/delx
          v1y=.5*(g01ex(i,j+1)-g01ex(i,j-1))/dely
          v2x=.5*(g02ex(i+1,j)-g02ex(i-1,j))/delx
          v2y=.5*(g02ex(i,j+1)-g02ex(i,j-1))/dely
          hax=.5*(h11ex(i+1,j)-h11ex(i-1,j))/delx
          hbx=.5*(h22ex(i+1,j)-h22ex(i-1,j))/delx
          hcx=.5*(h12ex(i+1,j)-h12ex(i-1,j))/delx
          hdx=.5*(h33ex(i+1,j)-h33ex(i-1,j))/delx
          hay=.5*(h11ex(i,j+1)-h11ex(i,j-1))/dely
          hby=.5*(h22ex(i,j+1)-h22ex(i,j-1))/dely
          hcy=.5*(h12ex(i,j+1)-h12ex(i,j-1))/dely
          hdy=.5*(h33ex(i,j+1)-h33ex(i,j-1))/dely
          g00x=.5*(g00(i+1,j)-g00(i-1,j))/delx
          g00y=.5*(g00(i,j+1)-g00(i,j-1))/dely
          g00xx=(g00(i+1,j)-2.*g00(i,j)+g00(i-1,j))/dxdx
          g00yy=(g00(i,j+1)-2.*g00(i,j)+g00(i,j-1))/dydy
          dgxy=g00(i+1,j+1)-g00(i-1,j+1)-g00(i+1,j-1)+g00(i-1,j-1)
          g00xy=.25*dgxy/dxdy

          h11p(i,j)=h11(i,j)
          h22p(i,j)=h22(i,j)
          h12p(i,j)=h12(i,j)
          h33p(i,j)=h33(i,j)

      s1 = 2*hc*v2x
      s2 = 4*ha*pypin*v2
      s3 = hay*v2
      s4 = 2*ha*v1x
      s5 = 4*ha*pxpin*v1
      s6 = hax*v1
      s7 = -2*a*g00y*pypin/(b*p4)
      s8 = 2*g00x*pxpin/p4
      s9 = -ay*g00y/(b*p4)/2.0
      s10 = -g00xx/p4
      s11 = ax*g00x/(a*p4)/2.0
      s12 = ric11(i,j)*g00(i,j)/p4
      s13 = ha*hd*g00(i,j)/d
      s14 = -2*hc**2*g00(i,j)/b
      s15 = ha*hb*g00(i,j)/b
      s16 = -ha**2*g00(i,j)/a
      sc11=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14+s15+s16

          h11(i,j)=h11(i,j)+delt*sc11

      s1 = hc*v2y
      s2 = hb*v2x
      s3 = 4*hc*pypin*v2
      s4 = hcy*v2
      s5 = ha*v1y
      s6 = hc*v1x
      s7 = 4*hc*pxpin*v1
      s8 = hcx*v1
      s9 = 2*g00x*pypin/p4
      s10 = 2*g00y*pxpin/p4
      s11 = bx*g00y/(b*p4)/2.0
      s12 = -g00xy/p4
      s13 = ay*g00x/(a*p4)/2.0
      s14 = ric12(i,j)*g00(i,j)/p4
      s15 = hc*hd*g00(i,j)/d
      s16 = -hb*hc*g00(i,j)/b
      s17 = -ha*hc*g00(i,j)/a
      sc12=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14+s15+s16+s17

          h12(i,j)=h12(i,j)+delt*sc12

      s1 = 2*hb*v2y
      s2 = 4*hb*pypin*v2
      s3 = hby*v2
      s4 = 2*hc*v1y
      s5 = 4*hb*pxpin*v1
      s6 = hbx*v1
      s7 = 2*g00y*pypin/p4
      s8 = -2*b*g00x*pxpin/(a*p4)
      s9 = -g00yy/p4
      s10 = by*g00y/(b*p4)/2.0
      s11 = -bx*g00x/(a*p4)/2.0
      s12 = ric22(i,j)*g00(i,j)/p4
      s13 = hb*hd*g00(i,j)/d
      s14 = -2*hc**2*g00(i,j)/a
      s15 = -hb**2*g00(i,j)/b
      s16 = ha*hb*g00(i,j)/a
      sc22=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14+s15+s16

          h22(i,j)=h22(i,j)+delt*sc22

      s1 = ric33(i,j)*g00(i,j)/p4
      s2 = 4*hd*pypin*v2
      s3 = hdy*v2
      s4 = 2*coty(i,j)*hd*v2
      s5 = 4*hd*pxpin*v1
      s6 = hdx*v1
      s7 = -2*d*g00y*pypin/(b*p4)
      s8 = -2*d*g00x*pxpin/(a*p4)
      s9 = -dy*g00y/(b*p4)/2.0
      s10 = -coty(i,j)*d*g00y/(b*p4)
      s11 = -dx*g00x/(a*p4)/2.0
      s12 = -hd**2*g00(i,j)/d
      s13 = hb*hd*g00(i,j)/b
      s14 = ha*hd*g00(i,j)/a
      sc33=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14

          h33(i,j)=h33(i,j)+delt*sc33
350     continue
351     continue

c  boundary conditions at z axis (y=0)
        do 360 i=2,imax+2
          h11p(i,1)=h11(i,1)
          h22p(i,1)=h22(i,1)
          h12p(i,1)=h12(i,1)
          h33p(i,1)=h33(i,1)
          h11(i,1)=h11(i,2)
          h22(i,1)=h22(i,2)
          h12(i,1)=-h12(i,2)
          h33(i,1)=h33(i,2)
c  boundary conditions at equator
          h11p(i,jmax+1)=h11(i,jmax+1)
          h22p(i,jmax+1)=h22(i,jmax+1)
          h12p(i,jmax+1)=h12(i,jmax+1)
          h33p(i,jmax+1)=h33(i,jmax+1)
          h11(i,jmax+1)=h11(i,jmax)
          h22(i,jmax+1)=h22(i,jmax)
          h12(i,jmax+1)=-h12(i,jmax)
          h33(i,jmax+1)=h33(i,jmax)
360     continue

c  boundary conditions at the throat (i=2)
        do 365 j=1,jmax+1
          h11p(1,j)=h11(1,j)
          h22p(1,j)=h22(1,j)
          h12p(1,j)=h12(1,j)
          h33p(1,j)=h33(1,j)
          h11(1,j)=h11(3,j)
          h22(1,j)=h22(3,j)
          h12(1,j)=-h12(3,j)
          h33(1,j)=h33(3,j)
c  boundary conditions at the outer edge of the grid (i=imax+3)
          h11p(imax+3,j)=h11(imax+3,j)
          h22p(imax+3,j)=h22(imax+3,j)
          h12p(imax+3,j)=h12(imax+3,j)
          h33p(imax+3,j)=h33(imax+3,j)
          h11(imax+3,j)=h11o(j)
          h22(imax+3,j)=h22o(j)
          h12(imax+3,j)=h12o(j)
          h33(imax+3,j)=h33o(j)
365     continue


c  extrapolate to the next half-timestep
c-------------------------------------------------------------------
        do 376 j=1,jmax+1
                do 375 i=1,imax+3
                        g01p(i,j)=g01(i,j)
                        g02p(i,j)=g02(i,j)
                        g11ex(i,j)=1.5*g11(i,j)-.5*g11p(i,j)
                        g22ex(i,j)=1.5*g22(i,j)-.5*g22p(i,j)
                        g33ex(i,j)=1.5*g33(i,j)-.5*g33p(i,j)
375             continue
376     continue


c  the scalar curvature 
c-------------------------------------------------------------------
        do 501 j=2,jmax
        do 500 i=2,imax+2
          psin=1./psi(i,j)
          psin2=psin*psin
          psin4=psin2*psin2
          psi4=1./psin4
          pxpin=psix(i,j)*psin
          pypin=psiy(i,j)*psin
          ax=.5*(g11ex(i+1,j)-g11ex(i-1,j))/delx
          bx=.5*(g22ex(i+1,j)-g22ex(i-1,j))/delx
          dx=.5*(g33ex(i+1,j)-g33ex(i-1,j))/delx
          ay=.5*(g11ex(i,j+1)-g11ex(i,j-1))/dely
          by=.5*(g22ex(i,j+1)-g22ex(i,j-1))/dely
          dy=.5*(g33ex(i,j+1)-g33ex(i,j-1))/dely
          d=g33ex(i,j)
          a=g11ex(i,j)
          b=g22ex(i,j)
          q=a*b
          ha=h11(i,j)
          hb=h22(i,j)
          hc=h12(i,j)
          hd=h33(i,j)
          sc1=b*b*ha*ha
          sc2=a*a*hb*hb
          sc3=2.*a*b*hc*hc
          scurv(i,j)=(sc1+sc2+sc3)/(q*q)+hd*hd/(d*d)
500     continue
501     continue


c  solve for the lapse
c-------------------------------------------------------------------
c  if the maximal slicing is chosen, solve the maximal slicing eq'n
        if (ich.eq.3)then
          do 532 j=2,jmax
          do 531 i=2,imax+2
            g00ppp(i,j)=g00pp(i,j)
            g00pp(i,j)=g00p(i,j)
            g00p(i,j)=g00(i,j)
            g00(i,j)=1.5*g00p(i,j)-0.5*g00ppp(i,j)
531       continue
532       continue

          call maximal (coty,g00,g11,g22,g33,psi,psix,psiy,scurv,nt)
        end if

c  if the lapse is determined by the determinant of the three metric
        if(ich.eq.4) then
          do 571 j=1,jmax+1
          do 570 i=1,imax+3
            detg=g33(i,2)*g11(i,2)*g22(i,2)
            detg2=detg*detg
            detg4=detg2*detg2
            g00(i,j)=detg
570       continue
571       continue
        end if



c  solve the shift equation
c------------------------------------------------------------------
        call shiftv (chi,g00,g01,g02,g11,g22,h12,nt)



c  setup to compute radiation variables and to output data
c------------------------------------------------------------------
        time=float(nt)*delt/adm

        if(float(pt)*nprnt3d*adm-(float(nt)).le.0.5) then
        i3dprint=1
	write(6,*)' printing data at t = ',time

        var=' g00 '
        file= 'g00hdf'
        call writehdf(var,file,g00,time,im,jm)
        var=' g11 '
        file= 'g11hdf'
        call writehdf(var,file,g11,time,im,jm)
        var=' g22 '
        file= 'g22hdf'
        call writehdf(var,file,g22,time,im,jm)
        var=' g33 '
        file= 'g33hdf'
        call writehdf(var,file,g33,time,im,jm)
        var=' h11 '
        file= 'h11hdf'
        call writehdf(var,file,h11,time,im,jm)
        var=' h12 '
        file= 'h12hdf'
        call writehdf(var,file,h12,time,im,jm)
        var=' h22 '
        file= 'h22hdf'
        call writehdf(var,file,h22,time,im,jm)
        var=' h33 '
        file= 'h33hdf'
        call writehdf(var,file,h33,time,im,jm)
        var=' g01 '
        file= 'g01hdf'
        call writehdf(var,file,g01,time,im,jm)
        var=' g02 '
        file= 'g02hdf'
        call writehdf(var,file,g02,time,im,jm)
        pt=pt+1
        end if

c  compute apparent horizon and radiation variables
c--------------------------------------------------------------------
        if((float(it)*nprnt2d-(float(nt)).le.0.5).or.
     $     (nt.eq.nmax))then
          call radhor(nt,g11,g22,g33,h11,h22,h33,h12,psi,phi,g00,
     $                g01,g02,ric11,ric22,ric12,ric33,i3dprint)
          call zerilli(psi,g11,g22,g33,time,psix)
          i3dprint=0
          it=it+1
        endif

c  move forward to the next time step
1000  continue

c
c Done main evolution loop
c====================================================================
c

c
c  compute de/dt (radiative losses) from Zerilli function
c----------------------------------------------------------------------
c
        rewind(90)
        rewind(91)

        read(90,*,end=8555) t2,z22
        read(91,*,end=8555) t2,z42
        do 843 i=1,nmax
        t1=t2
        z21=z22
        z41=z42
        read(90,*,end=8555) t2,z22
        read(91,*,end=8555) t2,z42

        dedtz2=((z22-z21)/(t2-t1))**2/(32.*pi)
        dedtz4=((z42-z41)/(t2-t1))**2/(32.*pi)

        write(92,101) .5*(t2+t1),dedtz2
        write(93,101) .5*(t2+t1),dedtz4
        write(94,101) .5*(t2+t1),dedtz2+dedtz4
 843    continue
 8555   continue


c
c  integrate de/dt expressions to find total radiated energy
c----------------------------------------------------------------------
c
        rewind(92)
        rewind(93)
        rewind(94)
        rewind(96)
        rewind(98)

        do 833 j=1,7
          if(j.eq.4.or.j.eq.6)goto 833
          iunit=91+j
          ene=0.
          read(iunit,*,end=8455) t2,dedt2
          ene=.5*(dedt2)*t2+ene
          write(iunit-10,101) t2,ene

          do 844 i=1,nmax
            t1=t2
            dedt1=dedt2
            read(iunit,*,end=8455) t2,dedt2
            ene=.5*(dedt1+dedt2)*(t2-t1)+ene
            write(iunit-10,101) t2,ene
 844      continue
 8455     continue
          if(j.eq.1)write(6,*) 'total z2 energy radiated = ',ene
          if(j.eq.2)write(6,*) 'total z4 energy radiated = ',ene
          if(j.eq.3)write(6,*) 'total zt energy radiated = ',ene
          if(j.eq.4)write(6,*) 'total bra energy radiated = ',ene
          if(j.eq.5)write(6,*) 'total bri energy radiated = ',ene
          if(j.eq.6)write(6,*) 'total p4a energy radiated = ',ene
          if(j.eq.7)write(6,*) 'total p4i energy radiated = ',ene
 833    continue
 101    format(1x,e14.7,1x,e14.7)
        stop
        end







      subroutine initdat (g11,g22,g33,g01,g02,
     &g11p,g22p,g33p,g01p,g02p,g11o,g22o,g33o,h11,h22,h12,h33,
     &h11p,h22p,h12p,h33p,h11o,h12o,h22o,h33o,
     &psi,psix,psiy,psixx,psiyy,psixy,coty,ric11,ric22,ric12,ric33,phi)
c**********************************************************************
c
c  This subroutine sets up the initial data and solves the initial
c  value problem for the conformal factor.
c
c**********************************************************************
c
      include 'param.h'
      real coty(im,jm)
      real psi(im,jm),psix(im,jm),psixx(im,jm),phi(im,jm)
      real psixy(im,jm),psiy(im,jm),psiyy(im,jm)
      real g11(im,jm),g22(im,jm),g33(im,jm)
      real g11p(im,jm),g22p(im,jm),g33p(im,jm)
      real g11o(jm),g22o(jm),g33o(jm)
      real g01(im,jm),g02(im,jm)
      real g01p(im,jm),g02p(im,jm)
      real h11(im,jm),h22(im,jm),h12(im,jm),h33(im,jm)
      real h11p(im,jm),h22p(im,jm),h12p(im,jm),h33p(im,jm)
      real h11o(jm),h22o(jm),h12o(jm),h33o(jm)
      real ric11(im,jm),ric22(im,jm),ric12(im,jm),ric33(im,jm)
      real etah(im),xih(jm)
      real s(im,jm)
      character*10 var,file
  
c  Set up Brill wave and coordinate grid
      write(6,*)'Brill wave coefficients: x0=',x0,'  s0=',s0,' a=',
     1bwamp
      write(6,*)'Using sin**',pow,' q function'
      write(37,*)'Brill wave coefficients: x0=',x0,'  s0=',s0,' a=',
     1bwamp
      write(37,*)'Using sin**',pow,' q function'
      do 120 j=1,jmax+1
        y=(float(j)-1.5)*dely
        g11o(j)=1.
        g22o(j)=1.
        g33o(j)=1.
        h11o(j)=0.
        h12o(j)=0.
        h22o(j)=0.
        h33o(j)=0.
        if(i3bh.eq.0)then
          ang=1.0
        else
          ang=exp(-((y-0.0)/(6.*dely))**2)
        endif
        do 110 i=1,imax+3
          coty(i,j)=cos(y)/sin(y)
          x=(float(i)-2.)*delx
          arg1=(x-x0)*(x-x0)/s0**2
          arg2=(x+x0)*(x+x0)/s0**2
          st=exp(-arg1)+exp(-arg2)
          s(i,j)=bwamp*sin(y)**pow*ang*st

          g11p(i,j)=exp(2.*s(i,j))
          g22p(i,j)=exp(2.*s(i,j))
          g33p(i,j)=1.
          g01p(i,j)=0.
          g02p(i,j)=0.
          h11p(i,j)=0.
          h22p(i,j)=0.
          h12p(i,j)=0.
          h33p(i,j)=0.

          g11(i,j)=exp(2.*s(i,j))
          g22(i,j)=exp(2.*s(i,j))
          g33(i,j)=1.
          g01(i,j)=0.
          g02(i,j)=0.
          h11(i,j)=0.
          h22(i,j)=0.
          h12(i,j)=0.
          h33(i,j)=0.
 110    continue
 120  continue
 
c  Compute the conformal factor and its derivatives 
      call conform(s,psi,psix,psiy,psixx,psiyy,psixy,coty,phi)

c  Compute the Ricci tensor components
        do 326 j=2,jmax
        do 325 i=2,imax+2
        d=g33(i,j)
        a=g11(i,j)
        b=g22(i,j)

        ax=.5*(g11(i+1,j)-g11(i-1,j))/delx
        bx=.5*(g22(i+1,j)-g22(i-1,j))/delx
        dx=.5*(g33(i+1,j)-g33(i-1,j))/delx
        ay=.5*(g11(i,j+1)-g11(i,j-1))/dely
        by=.5*(g22(i,j+1)-g22(i,j-1))/dely
        dy=.5*(g33(i,j+1)-g33(i,j-1))/dely
        ayy=(g11(i,j+1)-2.*g11(i,j)+g11(i,j-1))/dydy
        bxx=(g22(i+1,j)-2.*g22(i,j)+g22(i-1,j))/dxdx
        dxx=(g33(i+1,j)-2.*g33(i,j)+g33(i-1,j))/dxdx
        dyy=(g33(i,j+1)-2.*g33(i,j)+g33(i,j-1))/dydy
        deltd2=g33(i+1,j+1)-g33(i-1,j+1)-g33(i+1,j-1)+g33(i-1,j-1)
        dxy=.25*deltd2/dxdy
        
        p=psi(i,j)
        pxx=psixx(i,j)
        pxy=psixy(i,j)
        pyy=psiyy(i,j)
        px=psix(i,j)
        py=psiy(i,j)
        pxpin=px/p
        pypin=py/p
        pxxpin=pxx/p
        pyypin=pyy/p
        pxypin=pxy/p

      s1 = -2*a*pyypin/b
      s2 = -2*a*pypin**2/b
      s3 = -a*dy*pypin/(b*d)
      s4 = -2*a*coty(i,j)*pypin/b
      s5 = a*by*pypin/b**2
      s6 = -2*ay*pypin/b
      s7 = -4*pxxpin
      s8 = 4*pxpin**2
      s9 = -dx*pxpin/d
      s10 = -bx*pxpin/b
      s11 = 2*ax*pxpin/a
      s12 = -ay*dy/(b*d)/4.0
      s13 = -dxx/d/2.0
      s14 = dx**2/d**2/4.0
      s15 = ax*dx/(a*d)/4.0
      s16 = -ay*coty(i,j)/b/2.0
      s17 = ay*by/b**2/4.0
      s18 = -bxx/b/2.0
      s19 = bx**2/b**2/4.0
      s20 = ax*bx/(a*b)/4.0
      s21 = -ayy/b/2.0
      s22 = ay**2/(a*b)/4.0

        r11a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
        r11b=s15+s16+s17+s18+s19+s20+s21+s22
        ric11(i,j)=r11a+r11b

      s1 = 6*pxpin*pypin
      s2 = bx*pypin/b
      s3 = -2*pxypin
      s4 = ay*pxpin/a
      s5 = dx*dy/d**2/4.0
      s6 = bx*dy/(b*d)/4.0
      s7 = -dxy/d/2.0
      s8 = -coty(i,j)*dx/d/2.0
      s9 = ay*dx/(a*d)/4.0
      s10 = bx*coty(i,j)/b/2.0

        r12a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10
        ric12(i,j)=r12a

      s1 = -4*pyypin
      s2 = 4*pypin**2
      s3 = -dy*pypin/d
      s4 = -2*coty(i,j)*pypin
      s5 = 2*by*pypin/b
      s6 = -ay*pypin/a
      s7 = -2*b*pxxpin/a
      s8 = -2*b*pxpin**2/a
      s9 = -b*dx*pxpin/(a*d)
      s10 = -2*bx*pxpin/a
      s11 = ax*b*pxpin/a**2
      s12 = -dyy/d/2.0
      s13 = dy**2/d**2/4.0
      s14 = -coty(i,j)*dy/d
      s15 = by*dy/(b*d)/4.0
      s16 = -bx*dx/(a*d)/4.0
      s17 = by*coty(i,j)/b/2.0
      s18 = ay*by/(a*b)/4.0
      s19 = -bxx/a/2.0
      s20 = bx**2/(a*b)/4.0
      s21 = ax*bx/a**2/4.0
      s22 = -ayy/a/2.0
      s23 = ay**2/a**2/4.0
      s24 = 1

        r22a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
        r22b=s15+s16+s17+s18+s19+s20+s21+s22+s23+s24
        ric22(i,j)=r22a+r22b

      s1 = -2*d*pyypin/b
      s2 = -2*d*pypin**2/b
      s3 = -2*dy*pypin/b
      s4 = -4*coty(i,j)*d*pypin/b
      s5 = by*d*pypin/b**2
      s6 = -ay*d*pypin/(a*b)
      s7 = -2*d*pxxpin/a
      s8 = -2*d*pxpin**2/a
      s9 = -2*dx*pxpin/a
      s10 = -bx*d*pxpin/(a*b)
      s11 = ax*d*pxpin/a**2
      s12 = -dyy/b/2.0
      s13 = dy**2/(b*d)/4.0
      s14 = -coty(i,j)*dy/b
      s15 = by*dy/b**2/4.0
      s16 = -ay*dy/(a*b)/4.0
      s17 = -dxx/a/2.0
      s18 = dx**2/(a*d)/4.0
      s19 = -bx*dx/(a*b)/4.0
      s20 = ax*dx/a**2/4.0
      s21 = by*coty(i,j)*d/b**2/2.0
      s22 = -ay*coty(i,j)*d/(a*b)/2.0
      s23 = d/b

        r33a=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11+s12+s13+s14
        r33b=s15+s16+s17+s18+s19+s20+s21+s22+s23
        ric33(i,j)=r33a+r33b

325     continue
326     continue

c  boundary conditions at z axis (y=0)
        do 390 i=2,imax+2
        ric11(i,1)=ric11(i,2)
        ric22(i,1)=ric22(i,2)
        ric12(i,1)=-ric12(i,2)
        ric33(i,1)=ric33(i,2)
c  boundary conditions at equator
        ric11(i,jmax+1)=ric11(i,jmax)
        ric22(i,jmax+1)=ric22(i,jmax)
        ric12(i,jmax+1)=-ric12(i,jmax)
        ric33(i,jmax+1)=ric33(i,jmax)
390     continue
c  boundary conditions at the throat (i=2)
        do 395 j=1,jmax+1
        ric11(1,j)=ric11(3,j)
        ric22(1,j)=ric22(3,j)
        ric12(1,j)=-ric12(3,j)
        ric33(1,j)=ric33(3,j)
c  boundary conditions at the outer edge of the grid (i=imax+3)
        ric11(imax+3,j)=ric11(imax+2,j)
        ric22(imax+3,j)=ric22(imax+2,j)
        ric12(imax+3,j)=ric12(imax+2,j)
        ric33(imax+3,j)=ric33(imax+2,j)
395     continue

  
c  Now write out the conformal factor to an hdf file
      do 140 j=1,jm
140     xih(j)=(float(j)-1.5)*dely
      do 130 i=1,im
130     etah(i)=(float(i)-2.)*delx
      call inithdf (etah,xih,im,jm)
      var=' psi '
      file='psihdf'
      call writehdf (var,file,psi,0.,im,jm)
      return
      end
 




 
 
      subroutine conform(s,psi,psix,psiy,psixx,psiyy,psixy,coty,phi)
c********************************************************************
c
c This subroutine calculates the conformal factor in conformal
c coordinates and returns it in the array psi(i,j).  The array s(i,j)
c contains the q function.  It also calculates the derivatives of psi
c and returns them in the appropriate arrays.
c
c********************************************************************
c
      include 'param.h'
      real coty(im,jm)
      real psi(im,jm),psix(im,jm),psixx(im,jm),phi(im,jm)
      real psixy(im,jm),psiy(im,jm),psiyy(im,jm)
      real s(im,jm),a(im,jm),b(im,jm),c(im,jm),rd(im,jm),ld(im,jm),
     &     rhs(im,jm)
 
      do 120 j=2,jmax
      do 110 i=2,imax+2
      sxx=(s(i+1,j)-2.*s(i,j)+s(i-1,j))/dxdx
      syy=(s(i,j+1)-2.*s(i,j)+s(i,j-1))/dydy
      cotthet=coty(i,j)
      a(i,j)=(0.25*(sxx+syy-1.)*dxdx-2.-2.*dxdx/dydy)
      b(i,j)=(1.+.5*dely*cotthet)*dxdx/dydy
      c(i,j)=(1.-.5*dely*cotthet)*dxdx/dydy
      ld(i,j)=1.
      rd(i,j)=1.
      rhs(i,j)=0.
 110  continue
 120  continue
c
c  Apply boundary conditions to coefficients and right hand side
c

c  Along z-axis use symmetric boundary condition
      j=2
      do 150 i=3,imax+1
      a(i,j)=a(i,j)+c(i,j)
      c(i,j)=0.
 150  continue

c  At corner of z-axis and outer boundary
      j=2
      i=imax+2
      a(i,j)=a(i,j)+c(i,j)
      a(i,j)=a(i,j)-delx*rd(i,j)
      ld(i,j)=ld(i,j)+rd(i,j)
      rhs(i,j)=rhs(i,j)-rd(i,j)*2.*delx*exp(xmax/2.)
      rd(i,j)=0.
      c(i,j)=0.

c  Along outer boundary (Robin condition)
      i=imax+2
      do 160 j=3,jmax-1
      a(i,j)=a(i,j)-delx*rd(i,j)
      ld(i,j)=ld(i,j)+rd(i,j)
      rhs(i,j)=rhs(i,j)-rd(i,j)*2.*delx*exp(xmax/2.)
      rd(i,j)=0.
      
 160  continue

c  At corner of outer boundary and equator
      i=imax+2
      j=jmax
      a(i,j)=a(i,j)-delx*rd(i,j)
      ld(i,j)=ld(i,j)+rd(i,j)
      rhs(i,j)=rhs(i,j)-rd(i,j)*2.*delx*exp(xmax/2.)
      a(i,j)=a(i,j)+b(i,j)
      rd(i,j)=0.
      b(i,j)=0.

c  Along equator (symmetric)
      j=jmax
      do 170 i=3,imax+1
      a(i,j)=a(i,j)+b(i,j)
      b(i,j)=0.
 170  continue

c  At corner of equator and throat 
      i=2
      j=jmax
      a(i,j)=a(i,j)+b(i,j)
      rd(i,j)=rd(i,j)+ld(i,j)
      b(i,j)=0.
      ld(i,j)=0.

c  Along throat (symmetric)
      i=2
      do 180 j=3,jmax-1
      rd(i,j)=rd(i,j)+ld(i,j)
      ld(i,j)=0.
 180  continue

c  At coner of throat and z-axis
      i=2
      j=2
      rd(i,j)=rd(i,j)+ld(i,j)
      a(i,j)=a(i,j)+c(i,j)
      c(i,j)=0.
      ld(i,j)=0.
      
c
c  Use spherical solution as initial guess
c
      do 187 j=1,jmax+1
      do 185 i=1,imax+3
      psi(i,j)=2.*cosh(float(i-2)*delx/2.)
185   continue
187   continue
      
c
c  Call elliptic equation solver
c
      eps=1.e-11
      
      call bicgst(2,imax+2,2,jmax,a,b,c,rd,ld,psi,rhs,eps,rmax,ier)

  
c  Use boundary conditions to obtain c.f.  on z-axis and equator
      do 210 i=2,imax+2
      psi(i,1)=psi(i,2)
      psi(i,jmax+1)=psi(i,jmax)
 210  continue
  
c  Boundary conditions on throat and outer boundary (Robin condition)
      do 220 j=1,jmax+1
      psi(1,j)=psi(3,j)
      psi(imax+3,j)=psi(imax+1,j)-psi(imax+2,j)*delx
     &               +2.*delx*exp(xmax/2.)
 220  continue
 
c  Now compute the first and second order derivatives of psi.
      do 161 j=2,jmax
      do 159 i=2, imax+2
        psix(i,j)=.5*(psi(i+1,j)-psi(i-1,j))/delx
        psiy(i,j)=.5*(psi(i,j+1)-psi(i,j-1))/dely
        psixx(i,j)=(psi(i+1,j)-2.*psi(i,j)+psi(i-1,j))/dxdx
        psiyy(i,j)=(psi(i,j+1)-2.*psi(i,j)+psi(i,j-1))/dydy
        delta2=psi(i+1,j+1)-psi(i-1,j+1)-psi(i+1,j-1)+psi(i-1,j-1)
        psixy(i,j)=.25*delta2/dxdy
159     continue
161   continue

c  Use boundary conditions to obtain c.f.  on z-axis and equator
      do 215 j=1,jmax+1
      do 214 i=1,imax+3
        phi(i,j)=4.*log(psi(i,j))
 214  continue
 215  continue
      return
      end







      subroutine maximal(coty,g00,g11,g22,g33,psi,psix,psiy,scurv,nt)
c**********************************************************************
c
c  This subroutine solves the maximal slicing equation in conformal
c  coordinates and returns the lapse in the array g00(i,j).
c
c**********************************************************************
c
      include 'param.h'
      real coty(im,jm)
      real psi(im,jm),psix(im,jm),psiy(im,jm)
      real g11(im,jm),g22(im,jm),g33(im,jm)
      real g00(im,jm)
      real scurv(im,jm)
      integer nt,i,j
      real a(im,jm),b(im,jm),c(im,jm),ld(im,jm),rd(im,jm),
     &     rhs(im,jm)
 
c  Determine the coefficients for the maximal slicing equation.
      do 120 j=2,jmax
      do 110 i=2,imax+2
      cotthet=coty(i,j)
      psi4in=1./psi(i,j)**4
      p4g11in=psi4in/g11(i,j)
      p4g22in=psi4in/g22(i,j)
      g11x=0.5*(g11(i+1,j)-g11(i-1,j))/delx
      g11y=0.5*(g11(i,j+1)-g11(i,j-1))/dely
      g22x=0.5*(g22(i+1,j)-g22(i-1,j))/delx
      g22y=0.5*(g22(i,j+1)-g22(i,j-1))/dely
      g33x=0.5*(g33(i+1,j)-g33(i-1,j))/delx
      g33y=0.5*(g33(i,j+1)-g33(i,j-1))/dely
      c1=(cotthet+2.*psiy(i,j)/psi(i,j)+g33y/(2.*g33(i,j))
     &     -g22y/(2.*g22(i,j))+g11y/(2.*g11(i,j)))*p4g22in
      c2=(2.*psix(i,j)/psi(i,j)+g33x/(2.*g33(i,j))
     &     +g22x/(2.*g22(i,j))-g11x/(2.*g11(i,j)))*p4g11in
      ld(i,j)=(-0.5*c2*delx + p4g11in)*dydy
      c(i,j)=(-0.5*c1*dely + p4g22in)*dxdx
      a(i,j)=-(dx2dy2*scurv(i,j) + 2.*dxdx*p4g22in
     &     +2.*dydy*p4g11in)
      b(i,j)=(0.5*c1*dely + p4g22in)*dxdx
      rd(i,j)=(0.5*c2*delx + p4g11in)*dydy
      rhs(i,j)=0.
 110  continue
 120  continue
 
c  Apply boundary conditions to coefficients and right hand side
c  along z-axis
      j=2
      do 150 i=3,imax+1
      a(i,j)=a(i,j)+c(i,j)
      c(i,j)=0.
 150  continue
      
c  corner of z-axis and outer boundary
      j=2
      i=imax+2
      a(i,j)=a(i,j)+c(i,j)
      rhs(i,j)=-rd(i,j)*tanh((xmax+delx)/2.)
      rd(i,j)=0.
      c(i,j)=0.

c  along outer boundary
      i=imax+2
      do 155 j=3,jmax-1
      rhs(i,j)=-rd(i,j)*tanh((xmax+delx)/2.)
      rd(i,j)=0.
 155  continue

c  corner of outer boundary and equator
      i=imax+2
      j=jmax
      rhs(i,j)=-rd(i,j)*tanh((xmax+delx)/2.)
      rd(i,j)=0.
      a(i,j)=a(i,j)+b(i,j)
      b(i,j)=0.

c  along equator
      j=jmax
      do 157 i=3,imax+1
      a(i,j)=a(i,j)+b(i,j)
      b(i,j)=0.
 157  continue

c  corner of equator and throat
      i=2
      j=jmax
      a(i,j)=a(i,j)+b(i,j)
      rd(i,j)=rd(i,j)+ld(i,j)
      ld(i,j)=0.
      b(i,j)=0.

c  along throat
      i=2
      do 160 j=3,jmax-1
      rd(i,j)=rd(i,j)+ld(i,j)
      ld(i,j)=0.
 160  continue

c  corner of throat and z-axis
      i=2
      j=2
      rd(i,j)=rd(i,j)+ld(i,j)
      a(i,j)=a(i,j)+c(i,j)
      c(i,j)=0.
      ld(i,j)=0.
 
c  Use previous solution as initial guess
c  Call elliptic equation solver

      eps=1.e-9
      
      call bicgst(2,imax+2,2,jmax,a,b,c,rd,ld,g00,rhs,eps,rmax,ier)


c  Use boundary conditions to obtain lapse on z-axis and equator
      do 210 i=2,imax+2
      g00(i,1)=g00(i,2)
      g00(i,jmax+1)=g00(i,jmax)
 210  continue
  
c  Boundary conditions on throat and outer boundary
      do 220 j=1,jmax+1
      g00(1,j)=g00(3,j)
      g00(imax+3,j)=tanh((xmax+delx)/2)
 220  continue
      return
      end








      subroutine shiftv (chi,g00,g01,g02,g11,g22,h12,nt)
c**********************************************************************
c
c  This subroutine calculates the shift vector in conformal coordinates
c  and returns it in the arrays g01(i,j) and g02(i,j).
c
c**********************************************************************
c
      include 'param.h'
      real chi(im,jm),chix,chiy
      real g11(im,jm),g22(im,jm)
      real g00(im,jm)
      real g01(im,jm),g02(im,jm)
      real h12(im,jm)
      real a(im,jm),b(im,jm),c(im,jm),ld(im,jm),rd(im,jm),
     &     rhs(im,jm)

c  Determine the coefficients for the shift vector potential problem.
      do 120 j=2,jmax
      do 110 i=2,imax+2
      ld(i,j)=g22(i,j)*dydy
      c(i,j)=g11(i,j)*dxdx
      a(i,j)=-2.*(g11(i,j)*dxdx+g22(i,j)*dydy)
      b(i,j)=c(i,j)
      rd(i,j)=ld(i,j)
 110  continue
 120  continue

c  Apply boundary conditions to coefficients and right hand side
c  along z-axis (anitsymmetric)
      j=2
      do 150 i=3,imax+1
      a(i,j)=a(i,j)-c(i,j)
      c(i,j)=0.
 150  continue

c  At corner of z-axis and outer boundary
      j=2
      i=imax+2
      a(i,j)=a(i,j)-c(i,j)
      rd(i,j)=0.
      c(i,j)=0.

c  along outer boundary
      i=imax+2
      do 183 j=3,jmax-1
      rd(i,j)=0.
 183  continue

c  at corner of outer boundary and equator
      i=imax+2
      j=jmax
      a(i,j)=a(i,j)-b(i,j)
      b(i,j)=0.
      rd(i,j)=0.

c  along equator (antisymmetric)
      j=jmax
      do 155 i=3,imax+1
      a(i,j)=a(i,j)-b(i,j)
      b(i,j)=0.
 155  continue

c  at corner of equator and throat
      i=2
      j=jmax
      a(i,j)=a(i,j)-b(i,j)
      rd(i,j)=rd(i,j)-ld(i,j)
      ld(i,j)=0.
      b(i,j)=0.

c  along throat (antisymmetric)
      i=2
      do 180 j=1,jmax-1
      rd(i,j)=rd(i,j)-ld(i,j)
      ld(i,j)=0.
 180  continue

c  at corner of throat and z-axis
      i=2
      j=2
      rd(i,j)=rd(i,j)-ld(i,j)
      a(i,j)=a(i,j)-c(i,j)
      c(i,j)=0.
      ld(i,j)=0.

      do 187 j=2,jmax
      do 185 i=2,imax+2
      rhs(i,j)=2.*g00(i,j)*h12(i,j)*dxdx*dydy
185   continue
187   continue

c  Use previous solution as initial guess
c  Call elliptic equation solver

      eps=1.e-10
      
      call bicgst(2,imax+2,2,jmax,a,b,c,rd,ld,chi,rhs,eps,rmax,ier)


c  Use boundary conditions to obtain potential on z-axis and equator
        do 645 i=2, imax+2
        chi(i,1)=-chi(i,2)
        chi(i,jmax+1)=-chi(i,jmax)
645     continue
        do 650 j=1, jmax+1
c  boundary conditions on throat and outer boundary
        chi(1,j)=-chi(3,j)
        chi(imax+3,j)=0.
650     continue
 
c  Now compute the shift vector from its potential.
        do 161 j=2,jmax
        do 159 i=2, imax+2
        chix=.5*(chi(i+1,j)-chi(i-1,j))/delx
        chiy=.5*(chi(i,j+1)-chi(i,j-1))/dely
        g01(i,j)=chiy
        g02(i,j)=chix
159     continue
161     continue

c  boundary conditions at z axis (y=0)
        do 670 i=2,imax+2
        g01(i,1)=g01(i,2)
        g02(i,1)=-g02(i,2)
c  boundary conditions at equator
        g01(i,jmax+1)=g01(i,jmax)
        g02(i,jmax+1)=-g02(i,jmax)
670     continue
c  boundary conditions at the throat (i=2)
        do 675 j=1,jmax+1
        g01(1,j)=-g01(3,j)
        g02(1,j)=g02(3,j)
c  boundary conditions at the outer edge of the grid (i=imax+3)
        g01(imax+3,j)=g01(imax+2,j)
        g02(imax+3,j)=g02(imax+2,j)
675     continue
      return
      end







      subroutine zerilli(psi,g11,g22,g33,time,psix)
c*********************************************************************
c
c  This routine calculates the projected wave amplitudes. Similarly
c  with the radial derivatives, determine the coefficients
c  and contributions to the integrand.
c
c  NOTE:  The normalization is just asin eqns 15-19 of PRD 45, p3544
c  (1992), except that here we divide by the result of the AMD mass of
c  the entire spacetime. Also, the time has been rescaled ADM mass,
c  whereas in the PRD article, the time was rescaled by the 'bare' mass
c  of the wormhole (m=2).
c     
c**********************************************************************
c
        include 'param.h'
        common /admmablk/ adm,igdet
	dimension sa(2*jmax,12),amp(12)
	real h1,h2,k,k1,k2,h1o,h2o,ko,k1o,k2o
	dimension rsch(im)
c       dimension rtort(im)
        dimension zpsi2(im),zpsi4(im)
        dimension g11(im,jm),g22(im,jm),g33(im,jm)
        dimension psi(im,jm),psix(im,jm)

c Compute zpsi2 and zpsi4 (\ell=2 and \ell=4 Zerilli functions) at all
c points in the spacetime (useful for making mivies if desired)
 	do 830 i=2,imax+1
          j=1
          y=(float(j)-1.5)*dely

	  projct0=sin(y)
	  projct1=sin(y)*(4.-3.*sin(y)**2)
 	  projct2=sin(y)*(3.*cos(y)**2-1.)
	  projct3=sin(y)*(4.-9.*sin(y)**2)
	  projct4=sin(y)*(35.*cos(y)**4-30.*cos(y)**2+3.)
	  projct5=sin(y)**3*(7.*cos(y)**2-1.)
	  projct6=sin(y)*(35.*sin(y)**4-60.*sin(y)**2+24.)
	  projct7=sin(y)*(175.*sin(y)**4-180.*sin(y)**2+24.)
	  projct8=sin(y)**3

	  psi2=psi(i,2)*psi(i,2)
	  psi4=psi2*psi2
	  psin=1./psi(i,2)
	  b=g22(i,2)
	  d=g33(i,2)
	  g1=psi4*g11(i,2)
	  g2=psi4*g22(i,2)
	  g3=psi4*g33(i,2)
	  bx=.5*(g22(i+1,2)-g22(i-1,2))/delx
	  dx=.5*(g33(i+1,2)-g33(i-1,2))/delx
	  dg2=(4.*psix(i,2)*psin*b+bx-2.*b)*psi4
	  dg3=(4.*psix(i,2)*psin*d+dx-2.*d)*psi4

	  sa(j,1)=projct8*(g2-g3)
	  sa(j,2)=g2*projct1+g3*projct3
	  sa(j,3)=projct2*g1
	  sa(j,4)=projct8*(dg2-dg3)
	  sa(j,5)=dg2*projct1+dg3*projct3
	  sa(j,6)=projct5*(g2-g3)
	  sa(j,7)=g2*projct6+g3*projct7
	  sa(j,8)=projct4*g1
	  sa(j,9)=projct5*(dg2-dg3)
	  sa(j,10)=dg2*projct6+dg3*projct7
	  sa(j,11)=g2*projct0
	  sa(j,12)=g1*projct0

	  do 810 j=2,jmax
	    y=(float(j)-1.5)*dely

	    projct0=sin(y)
	    projct1=sin(y)*(4.-3.*sin(y)**2)
	    projct2=sin(y)*(3.*cos(y)**2-1.)
	    projct3=sin(y)*(4.-9.*sin(y)**2)
	    projct4=sin(y)*(35.*cos(y)**4-30.*cos(y)**2+3.)
	    projct5=sin(y)**3*(7.*cos(y)**2-1.)
	    projct6=sin(y)*(35.*sin(y)**4-60.*sin(y)**2+24.)
	    projct7=sin(y)*(175.*sin(y)**4-180.*sin(y)**2+24.)
	    projct8=sin(y)**3

            psi2=psi(i,j)*psi(i,j)   
            psi4=psi2*psi2
	    psin=1./psi(i,j)
	    b=g22(i,j)
	    d=g33(i,j)
            g1=psi4*g11(i,j)
            g2=psi4*g22(i,j)
            g3=psi4*g33(i,j)
	    bx=.5*(g22(i+1,j)-g22(i-1,j))/delx
	    dx=.5*(g33(i+1,j)-g33(i-1,j))/delx
	    dg2=(4.*psix(i,j)*psin*b+bx-2.*b)*psi4
	    dg3=(4.*psix(i,j)*psin*d+dx-2.*d)*psi4

	    sa(j,1)=projct8*(g2-g3)
	    sa(j,2)=g2*projct1+g3*projct3
	    sa(j,3)=projct2*g1
	    sa(j,4)=projct8*(dg2-dg3)
	    sa(j,5)=dg2*projct1+dg3*projct3
	    sa(j,6)=projct5*(g2-g3)
	    sa(j,7)=g2*projct6+g3*projct7
	    sa(j,8)=projct4*g1
	    sa(j,9)=projct5*(dg2-dg3)
	    sa(j,10)=dg2*projct6+dg3*projct7
	    sa(j,11)=g2*projct0
	    sa(j,12)=g1*projct0
810	  continue
          do 812 j=1,jmax
            do 811 ni=1,12
               sa(jmax+j,ni)=sa(jmax-j+1,ni)
               y=(float(jmax+j)-1.5)*dely
811         continue
812       continue

	  j=2*jmax
   	  y=pi

	  projct0=sin(y)
	  projct1=sin(y)*(4.-3.*sin(y)**2)
	  projct2=sin(y)*(3.*cos(y)**2-1.)
	  projct3=sin(y)*(4.-9.*sin(y)**2)
	  projct4=sin(y)*(35.*cos(y)**4-30.*cos(y)**2+3.)
	  projct5=sin(y)**3*(7.*cos(y)**2-1.)
	  projct6=sin(y)*(35.*sin(y)**4-60.*sin(y)**2+24.)
	  projct7=sin(y)*(175.*sin(y)**4-180.*sin(y)**2+24.)
	  projct8=sin(y)**3
 
	  psi2=psi(i,jmax)*psi(i,jmax)
	  psi4=psi2*psi2
	  psin=1./psi(i,jmax)
	  b=g22(i,jmax)
	  d=g33(i,jmax)
	  g1=psi4*g11(i,jmax)
	  g2=psi4*g22(i,jmax)
	  g3=psi4*g33(i,jmax)
	  bx=.5*(g22(i+1,jmax)-g22(i-1,jmax))/delx
	  dx=.5*(g33(i+1,jmax)-g33(i-1,jmax))/delx
	  dg2=(4.*psix(i,jmax)*psin*b+bx-2.*b)*psi4
	  dg3=(4.*psix(i,jmax)*psin*d+dx-2.*d)*psi4

	  sa(j,1)=projct8*(g2-g3)
	  sa(j,2)=g2*projct1+g3*projct3
	  sa(j,3)=projct2*g1
	  sa(j,4)=projct8*(dg2-dg3)
	  sa(j,5)=dg2*projct1+dg3*projct3
	  sa(j,6)=projct5*(g2-g3)
	  sa(j,7)=g2*projct6+g3*projct7
	  sa(j,8)=projct4*g1
	  sa(j,9)=projct5*(dg2-dg3)
	  sa(j,10)=dg2*projct6+dg3*projct7
	  sa(j,11)=g2*projct0
	  sa(j,12)=g1*projct0

c  do the integrations over the sphere using Simpson's rule
          do 821 ni=1,12
            amp(ni)=0.
            amp(ni)=amp(ni)+(sa(1,ni)+sa(2*jmax-1,ni))*dely/3.
            do 820 j=2,2*jmax-4,2
              amp(ni)=amp(ni)+(4.*sa(j,ni)+2.*sa(j+1,ni))*dely/3.
820         continue
            amp(ni)=amp(ni)+4.*sa(2*jmax-2,ni)*dely/3.
c  add a trapazoidal rule part
            amp(ni)=amp(ni)+.5*sa(2*jmax-1,ni)*dely
821       continue

c  the normalization factor is mult by 1/2 since integ is from 0 to 2pi
          rsch2=.5*amp(11)
	  if(rsch2.le.0.) then
	    rsch2 = 100.
	    iscrewup=1
          else
	    rsch(i)=sqrt(rsch2)
 	  endif
	  r=rsch(i)
 	  r2=r*r
c Assume we are perturbing a blak hole of mass = 2.
	  smass=2.
	  sfact=1.-2.*smass/r
	  fact2=.125*sqrt(5.*pi)
	  fact4=.0625*sqrt(pi)
          k=fact2*amp(2)/r2
          h2=4.*fact2*amp(3)*sfact/r2
          ko=fact4*amp(7)/r2
          h2o=6.*fact4*amp(8)*sfact/r2
	
c  calculate auxiliary functions h1, h2, g, k
          h1=0.
          h1o=0.
c  calculate the radial derivatives of G and K times r
          rgr=fact2*amp(4)/r2
          rkr=fact2*amp(5)/r2
          rgro=fact4*amp(9)/r2
          rkro=fact4*amp(10)/r2
       
c  calculate the functions k1 and k2
          k1=k+sfact*rgr-2.*sfact*h1/r
          k2=.5*(h2-k-rkr+k*smass/(r-2.*smass))/sfact
          k1o=ko+sfact*rgro-2.*sfact*h1o/r
          k2o=.5*(h2o-ko-rkro+ko*smass/(r-2.*smass))/sfact
       
c  calculate the Zerilli wave-function
          zpsi2(i)=(4.*r*sfact*sfact*k2+6.*r*k1)/(1.+1.5*smass/r)
          zpsi4(i)=(4.*r*sfact*sfact*k2o+20.*r*k1o)/(1.+smass/(3.*r))

c         if (r.lt.4.)then
c           rtort(i)=0.
c          else
c           rtort(i)=r+2.*smass*log(.5*r/smass-1.)
c         end if
830     continue

c Calculate the value of eta for given value of r in schwatrzchild
c coordinates. This is given by:
c
c      eta = ln(2rp/m)
c       rp = isotropic radial coordinate
c        r = schwatrzchild radial coordinate
c        r = ( 1 + m/(2rp))*rp
c
c Then find the index of that value of eta, i.e., eta/delx + 2
c need to add 2 since eta=0 is i=2
c
c Detectors sit a distances of 15M, 30M and 60M in Schwartzchild
c coordinates. Find the eta index for these distances.
c
        i15=nint( (log((15.-1.)+sqrt(15.*15.-2.))) / delx ) + 2
        i30=nint( (log((30.-1.)+sqrt(30.*30.-2.))) / delx ) + 2
        i60=nint( (log((60.-1.)+sqrt(60.*60.-2.))) / delx ) + 2

c
c One may find it useful to print out data at several shells.
c rtort is the so-called tortoise coordinate useful in black
c hole perturbation theory. The declaration for rtort at the top of
c this routine must first be uncommented. Also, the if block just
c preceeding the line labeled 830 must be uncommented.

c       write(34,*)time,rsch(i15),rtort(i15),
c    $	  	   zpsi2(i15),zpsi4(i15)
c       write(35,*)time,rsch(i30),rtort(i30),
c    $		   zpsi2(i30),zpsi4(i30)
c       write(36,*)time,rsch(i60),rtort(i60),
c    $		   zpsi2(i60),zpsi4(i60)

        if(idet.eq.15)i=i15
        if(idet.eq.30)i=i30
        if(idet.eq.60)i=i60

        write(90,101) time,zpsi2(i)/(sqrt(12.)*adm)
        write(91,101) time,zpsi4(i)/(sqrt(180.)*adm)

 101    format(1x,e14.7,1x,e14.7)
	if(iscrewup.eq.1)write(37,*)'r negative in zerilli'
       return
       end





