*
*  The following parameters define the coordinate grid and time spacing.
************************************************************************
*
*    xmax = value of etamax (the value of eta at the outermost zone)
*    imax = number of spatial steps between eta=0 and eta=xmax
*    delx = coordinate spacing in eta direction
*           (i=2       => eta=0)
*           (i=imax+2  => eta=xmax)
*
*    ymax = value of thetamax (the maximum angular value in the grid)
*    jmax = number of spatial steps between theta=0 and theta=ymax
*           (Note: Coordinate grid straddles theta=0)
*    dely = coordinate spacing in theta direction
*           (j=2       => theta=+dely/2)
*           (j=jmax+1  => theta = pi/2 + dely/2)
*
*    tmax = total time for run in units of ADM mass
*         = -1 if no evolution 
*         =  0 if one iteration
*    delt = time increment size
*
*    ich  = defines time slicing condition
*         = 2 : initially, at t=0, g00=tanh(.5*eta) (else g00=1)
*         = 3 : maximal slicing (evolution) with symmetric bc at throat
*         = 4 : use determinant of 3-metric to define g00
*
* 
*-----------------------------------------------------------------------
*
*               For Low resolution set imax=100
*               For Med resolution set imax=200
*
*-----------------------------------------------------------------------
*
*
*  The following parameters define the frequency and type of data dumps
************************************************************************
*
*    nprnt2d   = time increment for computing 2d data  (in iterations #)
*    prnt3d    = time increment for computing 3d data  (in units of ADM)
*    i3dprint0 = 1 if 3D data to be printed at t=0 (0 otherwise)
*
*    ihor      = 1 to find apparent horizon 
*              = 0 otherwise
*    ibrinvnp  = 1 to compute bel-robinson vectors, invariants 
*                  and Newman-Penrose scalars from "electric" and
*                  "magnetic" pieces of the Weyl tensor
*              = 0 otherwise
*    iyor      = 1 to compute Cotton-York tensor components 
*              = 0 otherwise
*    idet      = location (in M) of gravitational wave detector
*
************************************************************************
*
      real nprnt2d,prnt3d,rimax,tmax,pi,xmax,ymax,rjmax
      integer imax,jmax,idet,ihor,ibrinvnp,iyor,i3dprint0,ich

c  variable parameters
c-----------------------------------------------------------------------
c
c     parameter(imax = 100)
      parameter(imax = 200)

      parameter(tmax      = 10.0              )
      parameter(prnt3d    = 10.0              )
      parameter(idet      = 15                )
      parameter(ihor      = 1                 )
      parameter(ibrinvnp  = 1                 )
      parameter(iyor      = 0                 )
      parameter(i3dprint0 = 1                 )
      parameter(ich       = 3                 )
      parameter(pi        = 3.1415926535897932)
      parameter(xmax      = 6.0               )
      parameter(ymax      = pi/2.0            )

      parameter(jmax=imax*ymax/xmax+2)

      parameter(rimax     = imax              )
      parameter(rjmax     = jmax              )
      parameter(nprnt2d   = 10.0*rimax/100.   )

*
*  Define the initial conditions for the Brill wave
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*  x0    = initial position of Brill wave
*  s0    = initial width of Brill wave
*  bwamp = initial amplitude of Brill wave
*  pow   = exponent of sin in q function
*  i3bh  = 1 to simulate 3 black holes (does not work) 
*        = 0 otherwise
*
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* 
      real x0,s0,bwamp
      integer pow,i3bh

      parameter(x0=1.0, s0=1.0, bwamp=1.0,  pow=2, i3bh=0)
c     parameter(x0=2.5, s0=.5,  bwamp=10.0, pow=2, i3bh=1)

 
c
c  fixed parameters
c-----------------------------------------------------------------------
c  (should NOT have to adjust)
c
      real delx,dely,dxdx,dydy,dx2dy2,dxdy,delt,hfpi,deta,dtheta
      real xi,nprnt3d
      integer im,jm,idim,jdim,nthtamax,netamax
      integer nem1,nep1,ntm1,ntp1,ntp2,ntm2,nep2,ngridpts

      parameter(    im = imax+3,       jm = jmax+1         )
      parameter(  idim = imax+1,     jdim = jmax-1         )
      parameter(  delx = xmax/rimax, dely = ymax/(rjmax-1.))
      parameter(  dxdx = delx**2,    dydy = dely**2        ) 
      parameter(dx2dy2 = dxdx*dydy,  dxdy = delx*dely      )
      parameter(  delt = delx                              )

      parameter(    hfpi = pi/2.,         deta = delx      )
      parameter(  dtheta = dely,            xi = delx/dely )
      parameter(nthtamax = jmax,       netamax = imax+2    )
      parameter(    nem1 = netamax-1,     nep1 = netamax+1 )
      parameter(    ntm1 = nthtamax-1,    ntp1 = nthtamax+1)
      parameter(    ntp2 = nthtamax+2,    ntm2 = nthtamax-2)
      parameter(    nep2 = netamax+2                       )
      parameter(ngridpts = netamax*nthtamax                )
      parameter( nprnt3d = prnt3d/delt                     )

