      
      subroutine radhor(nt,g11,g22,g33,h11,h22,h33,h12,fsi,phi,g00,
     $                  g01,g02,ric11,ric22,ric12,ric33,i3dprint)
c*********************************************************************
c
c  main subroutine to direct calculations of
c     Riemann invariants
c     Bel-Robinson vectors 
c     Newman-Penrose scalars 
c     Cotton-York tensor components 
c     apparent horizon (location, mass, embedding and curvature)
c
c
c---------------------------------------------------------------------
c  definitions of key variables:
c
c     h      = coordinate defining location of trapped surface
c     adm    = ADM mass
c     ahm    = apparent horizon mass
c     aharea = apparent horizon area
c     bre    = Bel-Ronbinson vector component in eta direction
c     brt    = Bel-Robinson vector component in theta direction
c     brn    = magnitude of Bel-Robinson vector
c     crvi   = Riemann invariant I
c     crvj   = Riemann invariant J
c     y13    = Cotton-York tensor component
c     y23    = Cotton-York tensor component
c     psi0   = Newman-Penrose scalar from "E" and "B" components
c     psi1   = Newman-Penrose scalar from "E" and "B" components
c     psi2   = Newman-Penrose scalar from "E" and "B" components
c     psi3   = Newman-Penrose scalar from "E" and "B" components
c     psi4   = Newman-Penrose scalar from "E" and "B" components
c     psi0d  = Newman-Penrose scalar from direct construction
c     psi1d  = Newman-Penrose scalar from direct construction
c     psi2d  = Newman-Penrose scalar from direct construction
c     psi3d  = Newman-Penrose scalar from direct construction
c     psi4d  = Newman-Penrose scalar from direct construction
c     psi4i  = integrated flux quantity to compute losses from psi4 
c     brei   = integrated flux quantity to compute losses from bre
c     zpsi2  = zerilli function (l=2)
c     zpsi4  = zerilli function (l=4)
c     zpsit  = zerilli function (sum total of l=2 and l=4)
c     igdet  = integer locating detector along the eta grid 
c
c
c----------------------------------------------------------------------
c  description of output files:
c
c     $$$hdf      = hdf files of variable $$$ on the 2D grid
c     $$$.xdat    = radial profile of variable $$$ along the equatorial
c                   axis at different times
c     $$$.tdat    = value of variable $$$ at specific radial distance
c                   (defined by idet) along the equator at diferent 
c                   times.
c     dedt$$$.dat = time vs. (dE/dt) for variable $$$
c     te$$$.dat   = time vs. total radiated energy (E) for variable $$$
c     zpsi2.dat   = time vs. l=2 zerilli function at distance idet
c     zpsi4.dat   = time vs. l=4 zerilli function at distance idet
c     embd.dat    = embedding of apparent horizon
c                   (rho, z, gaussian curvature)
c     cpoe.dat    = time vs. ratio of polar to equatorial 
c                   circumference of apparent horizon
c     hor.dat     = theta vs. eta location of apparent horizon
c     ahm.dat     = time vs. apparent horizon mass
c     info.dat    = grid and physical parameters of run
c
c*******************************************************************
c
      include 'param.h'
      common /metrcblk/ a(im,jm),b(im,jm),d(im,jm)
      common /excrvblk/ ha(im,jm),hb(im,jm),hc(im,jm),hd(im,jm)
      common /gdervblk/ a1n(im,jm),b1n(im,jm),d1t(im,jm),
     &                  d1n(im,jm),a1t(im,jm),b1t(im,jm), 
     &                  b2n(im,jm),d2n(im,jm),a2n(im,jm),
     &                  a2t(im,jm),b2t(im,jm),d2m(im,jm),
     &                  d2t(im,jm),a2m(im,jm),b2m(im,jm)
      common /hdervblk/ ha1n(im,jm),ha1t(im,jm),hb1n(im,jm),
     &                  hb1t(im,jm),hc1n(im,jm),hc1t(im,jm),
     &                  hd1n(im,jm),hd1t(im,jm)
      common /cnfmlblk/ psim4(im,jm),psi(im,jm),phi1n(im,jm),
     &                  phi2n(im,jm),phi1t(im,jm),
     &                  phi2t(im,jm),phi2m(im,jm)
      common /gaugeblk/ alpha(im,jm),beta1u(im,jm),
     &                  beta2u(im,jm),alpha1t(im,jm),
     &                  beta1u1t(im,jm),beta2u1n(im,jm),
     &                  beta1u2t(im,jm),beta2u2m(im,jm)
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /gridblck/ eta(im),theta(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /radiablk/ bre(im,jm),brt(im,jm),brn(im,jm),ssj(im,jm),
     &                  crvi(im,jm),crvj(im,jm),ssi(im,jm),
     &                  y13(im,jm),y23(im,jm),psi0(im,jm),psi1(im,jm),
     &                  psi2(im,jm),psi3(im,jm),psi4(im,jm)
      common /weylblck/ psi0d(im,jm),psi1d(im,jm),psi2d(im,jm),
     &                  psi3d(im,jm),psi4d(im,jm)
      common /energblk/ psi4i(im,jm),brei(im,jm)
      common /timeblck/ scaleprm,jtime,times,timesold
      common /admmablk/ adm,igdet
      common/embdblk/ z0,zz(jm),rho(jm),gauss(jm)

      dimension g11(im,jm),g22(im,jm),g33(im,jm),fsi(im,jm)
      dimension h11(im,jm),h22(im,jm),h33(im,jm),h12(im,jm)
      dimension ric11(im,jm),ric22(im,jm),ric12(im,jm),ric33(im,jm)
      dimension tic11(im,jm),tic22(im,jm),tic12(im,jm),tic33(im,jm)
      dimension g00(im,jm),g01(im,jm),g02(im,jm),phi(im,jm)
      dimension hadm(2*jmax-1),p4i(2*jmax-1)
      dimension bri(2*jmax-1)

c  setup data and parameters at the initial time step
c------------------------------------------------------------------
      jtime=nt
      if(nt.eq.0)then
        times=0.
        scaleprm=2.0
        sqrt2=sqrt(2.)
        sqrthfm=sqrt(.5*scaleprm)

c  setup guess for locating the apparent horizon on initial slice
        h0=0.01
        do 33 j=1,nthtamax
          h(j)=h0
  33    continue
      
        do 73 j=1,nthtamax
        do 74 i=1,netamax-1
          psi4i(i,j)=0.
          brei(i,j)=0.
 74     continue
 73     continue
  
        do 10 i=1,netamax
          t=float(i-1)*deta
          eta(i)=t
   10   continue

        do 20 j=1,nthtamax
          t=(float(j)-1.5)*dtheta
          theta(j)=t
          sn(j)=sin(t)
          cs(j)=cos(t)
          ct(j)=cs(j)/sn(j)
   20   continue

        if(ihor.eq.1)then
          write(70,*) int(tmax/prnt3d +0.5)+1
          write(70,*) nthtamax
          write(72,*) int(tmax/prnt3d +0.5)+1
          write(72,*) nthtamax
          do 43 j=1,nthtamax
            write(70,*) theta(j)
            write(72,*) theta(j)
 43       continue
        endif

c  place "detector" on grid (at r = 15, 30, 60 ADM)
c--------------------------------------------------------------------
c Calculate the value of eta for given value of r in Schwarzschild
c coordinates. This is given by:
c
c      eta = ln(2rp/m)
c       rp = isotropic radial coordinate
c        r = Schwarzschild 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 at distances of 15M, 30M and 60M in Schwarzschild
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

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

        write(6,*) 'detector at ',idet,' M, ',' i= ',igdet
      endif


c define metric, extrinsic curvature, lapse and shift variables
c and their derivatives
c----------------------------------------------------------------------
      do 31 j=1,nthtamax+1
      do 30 i=2,netamax+1
        im1=i-1
        jm1=j-1
        ip1=i+1
        jp1=j+1
        if(jm1.eq.0)jm1=3
        if(jp1.eq.nthtamax+2)jp1=nthtamax-1
        if(ip1.eq.netamax+2)ip1=netamax+1

        a(im1,j)=g11(i,j)
        b(im1,j)=g22(i,j)
        d(im1,j)=g33(i,j)
        psi(im1,j)=fsi(i,j)
        psim4(im1,j)=fsi(i,j)**(-4.)
        alpha(im1,j)=g00(i,j)
        beta1u(im1,j)=g01(i,j)
        beta2u(im1,j)=g02(i,j)
        tic11(im1,j)=ric11(i,j)
        tic12(im1,j)=ric12(i,j)
        tic22(im1,j)=ric22(i,j)
        tic33(im1,j)=ric33(i,j)

        ha(im1,j)=h11(i,j)
        hb(im1,j)=h22(i,j)
        hd(im1,j)=h33(i,j)
        hc(im1,j)=h12(i,j)

        a1n(im1,j)=(g11(ip1,j)-g11(im1,j))/(2.*deta)
        b1n(im1,j)=(g22(ip1,j)-g22(im1,j))/(2.*deta)
        d1n(im1,j)=(g33(ip1,j)-g33(im1,j))/(2.*deta)
        phi1n(im1,j)=(phi(ip1,j)-phi(im1,j))/(2.*deta)
        beta2u1n(im1,j)=(g02(ip1,j)-g02(im1,j))/(2.*deta)

        a1t(im1,j)=(g11(i,jp1)-g11(i,jm1))/(2.*dtheta)
        b1t(im1,j)=(g22(i,jp1)-g22(i,jm1))/(2.*dtheta)
        d1t(im1,j)=(g33(i,jp1)-g33(i,jm1))/(2.*dtheta)
        phi1t(im1,j)=(phi(i,jp1)-phi(i,jm1))/(2.*dtheta)
        alpha1t(im1,j)=(g00(i,jp1)-g00(i,jm1))/(2.*dtheta)
        beta1u1t(im1,j)=(g01(i,jp1)-g01(i,jm1))/(2.*dtheta)

        a2n(im1,j)=(g11(ip1,j)+g11(im1,j)-2.*g11(i,j))/deta**2
        b2n(im1,j)=(g22(ip1,j)+g22(im1,j)-2.*g22(i,j))/deta**2
        d2n(im1,j)=(g33(ip1,j)+g33(im1,j)-2.*g33(i,j))/deta**2
        phi2n(im1,j)=(phi(ip1,j)+phi(im1,j)-2.*phi(i,j))/deta**2

        a2t(im1,j)=(g11(i,jp1)+g11(i,jm1)-2.*g11(i,j))/dtheta**2
        b2t(im1,j)=(g22(i,jp1)+g22(i,jm1)-2.*g22(i,j))/dtheta**2
        d2t(im1,j)=(g33(i,jp1)+g33(i,jm1)-2.*g33(i,j))/dtheta**2
        phi2t(im1,j)=(phi(i,jp1)+phi(i,jm1)-2.*phi(i,j))/dtheta**2
        beta1u2t(im1,j)=(g01(i,jp1)+g01(i,jm1)-2.*g01(i,j))/dtheta**2

        a2m(im1,j)=(g11(ip1,jp1)-g11(ip1,jm1)-g11(im1,jp1)
     $             +g11(im1,jm1))/(4.*deta*dtheta)
        b2m(im1,j)=(g22(ip1,jp1)-g22(ip1,jm1)-g22(im1,jp1)
     $             +g22(im1,jm1))/(4.*deta*dtheta)
        d2m(im1,j)=(g33(ip1,jp1)-g33(ip1,jm1)-g33(im1,jp1)
     $             +g33(im1,jm1))/(4.*deta*dtheta)
        phi2m(im1,j)=(phi(ip1,jp1)-phi(ip1,jm1)-phi(im1,jp1)
     $               +phi(im1,jm1))/(4.*deta*dtheta)
        beta2u2m(im1,j)=(g02(ip1,jp1)-g02(ip1,jm1)-g02(im1,jp1)
     $                  +g02(im1,jm1))/(4.*deta*dtheta)

        ha1n(im1,j)=(h11(ip1,j)-h11(im1,j))/(2.*deta)
        hb1n(im1,j)=(h22(ip1,j)-h22(im1,j))/(2.*deta)
        hd1n(im1,j)=(h33(ip1,j)-h33(im1,j))/(2.*deta)
        hc1n(im1,j)=(h12(ip1,j)-h12(im1,j))/(2.*deta)

        ha1t(im1,j)=(h11(i,jp1)-h11(i,jm1))/(2.*dtheta)
        hb1t(im1,j)=(h22(i,jp1)-h22(i,jm1))/(2.*dtheta)
        hd1t(im1,j)=(h33(i,jp1)-h33(i,jm1))/(2.*dtheta)
        hc1t(im1,j)=(h12(i,jp1)-h12(i,jm1))/(2.*dtheta)
   30 continue
   31 continue

c  fix boundary values for derivatives of antisymmetric variables
      do 32 i=1,netamax
        hc1n(i,1)=-hc1n(i,2)
        hc1t(i,1)=hc1t(i,2)
        beta2u1n(i,1)=-beta2u1n(i,2)
        beta2u2m(i,1)=-beta2u2m(i,2)

        hc1n(i,ntp1)=-hc1n(i,nthtamax)
        hc1t(i,ntp1)=hc1t(i,nthtamax)
        beta2u1n(i,ntp1)=-beta2u1n(i,nthtamax)
        beta2u2m(i,ntp1)=-beta2u2m(i,nthtamax)
 32   continue
        

c  compute ADM mass and apparent horizon mass on initial time slice
c------------------------------------------------------------------
      if(nt.eq.0)then
        i=netamax-1
        do 11 j=1,jmax-1
          hadm(j)=.5*psi(i,j)*(.5*phi1n(i,j)-1.)*sn(j)
          hadm(2*jmax-j)=.5*psi(i,j+1)*(.5*phi1n(i,j+1)-1.)*sn(j+1)
   11   continue
        hadm(jmax)=.5*psi(i,jmax)*(.5*phi1n(i,jmax)-1.)*sn(jmax)

        call simpson(hadm,dtheta,2*jmax-1,ans)
c       ans=ans-2.*(.5*hadm(1)*dtheta/2.)
        adm=-sqrthfm*exp(0.5*eta(i))*ans
        write(6,*) 'adm mass = ',adm

        call aphor
        call ahmass
        write(6,*) 'ahm, ahm/adm, aharea = ',ahm*adm,ahm,aharea
        write(6,*) '1-ahm/adm = ',1.-ahm
        write(6,*)
      endif


c  define variables needed to compute radiative losses
      i=igdet
      rr=(sqrthfm)**2*exp(eta(i))/adm
      timesold=times
      times=float(nt)*delt/adm
      delts=times-timesold


c  compute Bel-Robinson vector, Newman-Penrose scalars and Riemann
c  invariants from "E" and "B" construction of Weyl tensor
c----------------------------------------------------------------------
      call brinvnp(tic11,tic22,tic12,tic33)

c  compute time integrated (i suffix) energy flux integral from psi4
c  and bel-robinson velocity at a fixed radial distance defined by idet
c----------------------------------------------------------------------
      do 54 j=1,nthtamax
        dum=sqrt(psi(i,j)**4*a(i,j))/adm
        psi4i(i,j)=psi4(i,j)*delts+psi4i(i,j)
        brei(i,j)=sign(1.,psi4(i,j))*sqrt(dum*abs(bre(i,j))/2.)*delts
     $           +brei(i,j)
 54   continue
     
      do 12 j=1,jmax-1
        p4i(j)=psi4i(i,j)*psi4i(i,j)*sn(j)
        p4i(2*jmax-j)=psi4i(i,j+1)*psi4i(i,j+1)*sn(j+1)

        bri(j)=brei(i,j)*brei(i,j)*sn(j)
        bri(2*jmax-j)=brei(i,j+1)*brei(i,j+1)*sn(j+1)
 12   continue
      p4i(jmax)=psi4i(i,jmax)*psi4i(i,jmax)*sn(jmax)
      bri(jmax)=brei(i,jmax)*brei(i,jmax)*sn(jmax)

c  perform spatial angular integrations using simpson integration
      call simpson(p4i,dtheta,2*jmax-1,ans)
c     ans=ans-2.*(.5*p4i(1)*dtheta/2.)
      dedtp4i=.5*rr**2*ans

      call simpson(bri,dtheta,2*jmax-1,ans)
c     ans=ans-2.*(.5*bri(1)*dtheta/2.)
      dedtbri=.5*rr**2*ans


c  generate output data files
c--------------------------------------------------------------------
c
c  compute energy loss (dE/dt) from bri and psi4
        write(96,101) times,dedtbri
        write(98,101) times,dedtp4i

c  compute horizon, mass and polar to equatorial circumference
        if(ihor.eq.1)then
          call aphor
          call ahmass
          call circumf(cpoe)
          write(64,101) times,ahm
          write(71,101) times,cpoe
        endif

c  write out radiation indicators at radius of detector along 
c  equatorial axis at different times
        if(ibrinvnp.eq.1)then
          write(50,101) times,bre(i,nthtamax)
          write(51,101) times,brt(i,nthtamax)
          write(52,101) times,brn(i,nthtamax)
          write(54,101) times,crvi(i,nthtamax)
          write(55,101) times,crvj(i,nthtamax)
          write(59,101) times,psi0(i,nthtamax)
          write(61,101) times,psi2(i,nthtamax)
          write(63,101) times,psi4(i,nthtamax)
        endif

c  compute Cotton-York tensor components
        if(iyor.eq.1)then
          call york
          write(57,101) times,y13(i,nthtamax)
          write(58,101) times,y23(i,nthtamax)
        endif


        if(i3dprint.eq.1)then

c  output HDF files
          call outrad

c  output radial spatial profiles along equatorial axis at fixed times
c  psi4 and psi0 are scaled by powers of the radius
          write(46,*)
          write(47,*)
          write(48,*)
          write(49,*)
          do 39 i=1,netamax-1
            rr=sqrthfm**2*exp(eta(i))/adm
            write(46,*) eta(i),psi4(i,nthtamax)*rr
            write(47,*) eta(i),psi0(i,nthtamax)*rr**3
            write(48,*) eta(i),psi4d(i,nthtamax)*rr
            write(49,*) eta(i),psi0d(i,nthtamax)*rr**3
 39       continue

c  compute apparent horizon embedding
          if(ihor.eq.1)then
            call embd
            write(70,*) times
            write(72,*) times
            do 44 j=1,nthtamax
              write(70,101) rho(j),zz(j),gauss(j)
              write(72,101) theta(j),h(j)
 44         continue
          endif
        endif

 101  format(5(e13.6,1x))
      return
      end




      subroutine outrad
c**********************************************************************
c
c  Output HDF files
c
c**********************************************************************
c
      include 'param.h'
      common /radiablk/ bre(im,jm),brt(im,jm),brn(im,jm),ssj(im,jm),
     &                  crvi(im,jm),crvj(im,jm),ssi(im,jm),
     &                  y13(im,jm),y23(im,jm),psi0(im,jm),psi1(im,jm),
     &                  psi2(im,jm),psi3(im,jm),psi4(im,jm)
      common /weylblck/ psi0d(im,jm),psi1d(im,jm),psi2d(im,jm),
     &                  psi3d(im,jm),psi4d(im,jm)
      common /gridblck/ eta(im),theta(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /timeblck/ scaleprm,jtime,times,timesold
      common /admmablk/ adm,igdet
c     dimension brer3(im,jm)
      character*10 var,file

c  scale by the spherical radius to some power
c     do 10 j=1,jm
c     do 11 i=1,im
c       rr=(sqrthfm)**2*exp(eta(i))/adm
c       brer3(i,j)=bre(i,j)*rr**3
c11   continue
c10   continue

c  output HDF files
      if(ibrinvnp.eq.1)then
        var=' bre   '
        file= 'brehdf'
        call writehdf(var,file,bre,times,im,jm)
        var=' brt '
        file= 'brthdf'
        call writehdf(var,file,brt,times,im,jm)
        var=' brn '
        file= 'brnhdf'
        call writehdf(var,file,brn,times,im,jm)
        var=' crvi '
        file= 'crvihdf'
        call writehdf(var,file,crvi,times,im,jm)
        var=' crvj '
        file= 'crvjhdf'
        call writehdf(var,file,crvj,times,im,jm)
        var=' psi0 '
        file= 'psi0hdf'
        call writehdf(var,file,psi0,times,im,jm)
        var=' psi2 '
        file= 'psi2hdf'
        call writehdf(var,file,psi2,times,im,jm)
        var=' psi4 '
        file= 'psi4hdf'
        call writehdf(var,file,psi4,times,im,jm)
        var=' psi0d'
        file= 'psi0dhdf'
        call writehdf(var,file,psi0d,times,im,jm)
        var=' psi2d'
        file= 'psi2dhdf'
        call writehdf(var,file,psi2d,times,im,jm)
        var=' psi4d'
        file= 'psi4dhdf'
        call writehdf(var,file,psi4d,times,im,jm)
      endif

      if(iyor.eq.1)then
        var=' y13 '
        file= 'y13hdf'
        call writehdf(var,file,y13,times,im,jm)
        var=' y23 '
        file= 'y23hdf'
        call writehdf(var,file,y23,times,im,jm)
      endif
      return
      end




      subroutine aphor
c**********************************************************************
c
c  main subroutine driver to locate trapped surfaces
c
c**********************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /diagblck/ dd(jm),du(jm),dl(jm)
      common /residblk/ res(jm)
      real rhs(jm)

      call diffh
      do 90 inr=1,50
        call lincoeff
        call residual
 
        do 10 j=2,nthtamax
          rhs(j)=-res(j)
  10    continue
        call tridiag(dd,du,dl,rhs,dh)
 
        call newh
        call check(iflag)
        if (iflag.eq.1) then
          go to 99
        endif
   90 continue
 
      print*,'apparent horizon finder did not converge'
   99 continue
      return
      end
 
 


      subroutine newh
c*******************************************************************
c
c  adjust h(j): the trapped surface coordinate
c
c*******************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      isum=0
      do 10 j=2,nthtamax
        h(j)=h(j)+dh(j)
        i=int(sign(1.,h(j)))
        isum=i+isum
   10 continue
      if (isum.eq.-(nthtamax-1)) then
        do 20 j=2,nthtamax
          h(j)=-h(j)
   20   continue
      endif
      h(1)=h(2)
      h(nthtamax+1)=h(nthtamax)
      call diffh
      return
      end
 
 
 

      subroutine lincoeff
c***********************************************************************
c
c compute linearized coefficients needed for Newton-Raphson iteration
c on the equations for trapped surfaces
c
c***********************************************************************
c
      include 'param.h'
      common /diagblck/ dd(jm),du(jm),dl(jm)
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
 
      call interpl3
 
c main diagonal
      do 10 j=2,nthtamax
        h2=h1t(j)*h1t(j)
        h3=h2*h1t(j)
        trk=hai(j)/ai(j)+hbi(j)/bi(j)+hdi(j)/di(j)
        trk1n=ha1ni(j)/ai(j)-hai(j)*a1ni(j)/(ai(j)*ai(j))
     &       +hb1ni(j)/bi(j)-hbi(j)*b1ni(j)/(bi(j)*bi(j))
     &       +hd1ni(j)/di(j)-hdi(j)*d1ni(j)/(di(j)*di(j))
        f=ai(j)*ai(j)*h2/bi(j)+ai(j)
        fh=sqrt(f)
        f1n=2.*ai(j)*a1ni(j)*h2/bi(j)-b1ni(j)*h2*(ai(j)/bi(j))**2.
     &     +a1ni(j)
 
        s1 =-2./(dtheta*dtheta)
        if(j.eq.2)s1=-1./(dtheta*dtheta)
        if(j.eq.nthtamax)s1=-1./(dtheta*dtheta)
 
        r1 =a1ni(j)*ct(j)/bi(j)
        r2 =-ai(j)*b1ni(j)*ct(j)/(bi(j)*bi(j))
        r3 =0.5*a1ni(j)*d1ti(j)/(bi(j)*di(j))
        r4 =0.5*ai(j)*d2mi(j)/(bi(j)*di(j))
        r5 =-0.5*ai(j)*d1ti(j)*b1ni(j)/(bi(j)*bi(j)*di(j))
        r6 =-0.5*ai(j)*d1ti(j)*d1ni(j)/(bi(j)*di(j)*di(j))
        r7 =0.5*a2mi(j)/bi(j)
        r8 =-0.5*a1ti(j)*b1ni(j)/(bi(j)*bi(j))
        r9 =phi2mi(j)*ai(j)/bi(j)
        r10=phi1ti(j)*a1ni(j)/bi(j)
        r11=-phi1ti(j)*ai(j)*b1ni(j)/(bi(j)*bi(j))
        s2=h3*(r1+r2+r3+r4+r5+r6+r7+r8+r9+r10+r11)
 
        r1=-0.5*d2ni(j)/di(j)
        r2=0.5*d1ni(j)*d1ni(j)/(di(j)*di(j))
        r3=-b2ni(j)/bi(j)
        r4=b1ni(j)*b1ni(j)/(bi(j)*bi(j))
        r5=0.5*a2ni(j)/ai(j)
        r6=-0.5*a1ni(j)*a1ni(j)/(ai(j)*ai(j))
        r7=-phi2ni(j)
        s3=h2*(r1+r2+r3+r4+r5+r6+r7)
 
        r1=0.5*d2mi(j)/di(j)
        r2=-0.5*d1ti(j)*d1ni(j)/(di(j)*di(j))
        r3=-0.5*b2mi(j)/bi(j)
        r4=0.5*b1ti(j)*b1ni(j)/(bi(j)*bi(j))
        r5=a2mi(j)/ai(j)
        r6=-a1ti(j)*a1ni(j)/(ai(j)*ai(j))
        r7=phi2mi(j)
        s4=h1t(j)*(r1+r2+r3+r4+r5+r6+r7)
 
        r1=-0.5*b1ni(j)*d1ni(j)/(ai(j)*di(j))
        r2=-0.5*bi(j)*d2ni(j)/(ai(j)*di(j))
        r3=0.5*bi(j)*d1ni(j)*a1ni(j)/(ai(j)*ai(j)*di(j))
        r4=0.5*bi(j)*d1ni(j)*d1ni(j)/(ai(j)*di(j)*di(j))
        r5=-0.5*b2ni(j)/ai(j)
        r6=0.5*b1ni(j)*a1ni(j)/(ai(j)*ai(j))
        r7=-phi1ni(j)*b1ni(j)/ai(j)
        r8=-phi2ni(j)*bi(j)/ai(j)
        r9=phi1ni(j)*bi(j)*a1ni(j)/(ai(j)*ai(j))
        s5=r1+r2+r3+r4+r5+r6+r7+r8+r9
 
        r1=-bi(j)*hai(j)/(ai(j)*ai(j))
        r2=-h2*hbi(j)/bi(j)
        r3=2.*hci(j)*h1t(j)/ai(j)
        s6=0.5*psii(j)*psii(j)*phi1ni(j)*fh*(r1+r2+r3)
 
        s7=0.5*psii(j)*psii(j)*(f1n/fh)*(r1+r2+r3)
 
        r1=-b1ni(j)*hai(j)/(ai(j)*ai(j))
        r2=-bi(j)*ha1ni(j)/(ai(j)*ai(j))
        r3=2.*bi(j)*hai(j)*a1ni(j)/(ai(j)*ai(j)*ai(j))
        r4=-h2*hb1ni(j)/bi(j)
        r5=h2*hbi(j)*b1ni(j)/(bi(j)*bi(j))
        r6=2.*hc1ni(j)*h1t(j)/ai(j)
        r7=-2.*hci(j)*a1ni(j)*h1t(j)/(ai(j)*ai(j))
        s8=psii(j)*psii(j)*fh*(r1+r2+r3+r4+r5+r6+r7)
 
        r1=b1ni(j)*psii(j)*psii(j)/(ai(j)*ai(j))
        r2=0.5*psii(j)*psii(j)*phi1ni(j)*bi(j)/(ai(j)*ai(j))
        r3=-2.*psii(j)*psii(j)*bi(j)*a1ni(j)/(ai(j)**3.)
        s9=(fh**3.)*trk*(r1+r2+r3)
 
        s10=1.5*bi(j)*psii(j)*psii(j)*f1n*fh*trk/(ai(j)*ai(j))
 
        s11=bi(j)*psii(j)*psii(j)*(fh**3.)*trk1n/(ai(j)*ai(j))
 
        dd(j)=s1+s2+s3+s4+s5+s6+s7+s8+s9+s10+s11

c  boundary elements
        alp=ai(j)*ct(j)/bi(j)+0.5*ai(j)*d1ti(j)/(bi(j)*di(j))+
     $      0.5*a1ti(j)/bi(j)+phi1ti(j)*ai(j)/bi(j)
        bet=-0.5*d1ni(j)/di(j)-b1ni(j)/bi(j)
     $      +0.5*a1ni(j)/ai(j)-phi1ni(j)
        gam=ct(j)+0.5*d1ti(j)/di(j)-0.5*b1ti(j)/bi(j)+
     $      a1ti(j)/ai(j)+phi1ti(j)
        rmu=bi(j)*psii(j)**2*trk/(ai(j)**2)
 
        s12=-3.*alp*h2/(2.*dtheta)
        s13=-bet*h1t(j)/dtheta
        s14=-gam/(2.*dtheta)
        s15=-1.5*rmu*fh*ai(j)**2*h1t(j)/(bi(j)*dtheta)
        s16=-0.5*psii(j)**2*ai(j)**2*h1t(j)/(fh*bi(j)*dtheta)
     $      *(-bi(j)*hai(j)/ai(j)**2-hbi(j)*h2/bi(j)
     $      +2.*hci(j)*h1t(j)/ai(j))
        s17=-psii(j)**2*fh*(-hbi(j)*h1t(j)/(bi(j)*dtheta)
     $      +hci(j)/(ai(j)*dtheta))
     
        if(j.eq.2)dd(j)=dd(j)+s12+s13+s14+s15+s16+s17
        if(j.eq.nthtamax)dd(j)=dd(j)-s12-s13-s14-s15-s16-s17
   10 continue
 
c upper and lower diagonals
c
      dtm1=1./dtheta
      do 20 j=2,nthtamax
        h2=h1t(j)*h1t(j)
        h3=h2*h1t(j)
        trk=hai(j)/ai(j)+hbi(j)/bi(j)+hdi(j)/di(j)
        f=ai(j)*ai(j)*h2/bi(j)+ai(j)
        fh=sqrt(f)
 
        s1=1./(dtheta*dtheta)
 
        r1=ai(j)*ct(j)/bi(j)
        r2=0.5*ai(j)*d1ti(j)/(bi(j)*di(j))
        r3=0.5*a1ti(j)/bi(j)
        r4=phi1ti(j)*ai(j)/bi(j)
        s2=1.5*h2*dtm1*(r1+r2+r3+r4)
 
        r1=-0.5*d1ni(j)/di(j)
        r2=-b1ni(j)/bi(j)
        r3=0.5*a1ni(j)/ai(j)
        r4=-phi1ni(j)
        s3=h1t(j)*dtm1*(r1+r2+r3+r4)
 
        r1=ct(j)
        r2=0.5*d1ti(j)/di(j)
        r3=-0.5*b1ti(j)/bi(j)
        r4=a1ti(j)/ai(j)
        r5=phi1ti(j)
        s4=0.5*dtm1*(r1+r2+r3+r4+r5)
 
        r1=-bi(j)*hai(j)/(ai(j)*ai(j))
        r2=-h2*hbi(j)/bi(j)
        r3=2.*hci(j)*h1t(j)/ai(j)
        s5=(0.5*ai(j)*ai(j)*psii(j)*psii(j)*h1t(j)*dtm1/(bi(j)*fh))
     &    *(r1+r2+r3)
 
        r1=-hbi(j)*h1t(j)/bi(j)
        r2=hci(j)/ai(j)
        s6=psii(j)*psii(j)*fh*dtm1*(r1+r2)
 
        s7=1.5*psii(j)*psii(j)*fh*h1t(j)*dtm1*trk
 
        du(j)=s1+s2+s3+s4+s5+s6+s7
        dl(j)=s1-s2-s3-s4-s5-s6-s7
   20 continue
      du(nthtamax)=0.
      dl(2)=0.
      return
      end
 
 
 
      subroutine tridiag(dg,dgu,dgl,rhs,ans)
c**********************************************************************
c
c  this subroutine inverts a tridiagonal metric with
c  diagonal dg, upper diagonal dgu, lower diagonal dgl,
c  and right hand side vector rhs.  the result is returned
c  in array ans.
c
c**********************************************************************
 
      include 'param.h'
      dimension dg(jm),dgu(jm),dgl(jm),rhs(jm),ans(jm)
      real gam(jm)
 
      bet=dg(2)
      ans(2)=rhs(2)/bet
      do 10 j=3,nthtamax
        gam(j)=dgu(j-1)/bet
        bet=dg(j)-dgl(j)*gam(j)
        ans(j)=(rhs(j)-dgl(j)*ans(j-1))/bet
   10 continue
      do 20 j=nthtamax-1,2,-1
        ans(j)=ans(j)-gam(j+1)*ans(j+1)
   20 continue
      return
      end
 
 
 
      subroutine diffh
c**********************************************************************
c
c  differentiate h(j): the trapped surface coordinate
c
c**********************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      real hx(0:nthtamax+2)
      do 5 j=2,nthtamax
        hx(j)=h(j)
    5 continue
      hx(1)=h(2)
      hx(0)=h(3)
      hx(ntp1)=h(nthtamax)
      hx(ntp2)=h(nthtamax-1)
      do 10 j=1,nthtamax
        h1t(j)=(hx(j+1)-hx(j-1))/(2.*dtheta)
        h2t(j)=(hx(j+1)-2.*hx(j)+hx(j-1))/(dtheta*dtheta)
   10 continue
      return
      end
 
 
 
      subroutine residual
c**********************************************************************
c
c  compute residual in Newton-Raphson method of solving for
c  trapped surfaces
c
c**********************************************************************
c
      include 'param.h'
      common /residblk/ res(jm)
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
 
      call interpl3
 
      do 10 j=2,nthtamax
        h2=h1t(j)*h1t(j)
        h3=h1t(j)*h2
        f=ai(j)*ai(j)*h2/bi(j)+ai(j)
        fh=sqrt(f)
        fh3=fh*fh*fh
        trk=hai(j)/ai(j)+hbi(j)/bi(j)+hdi(j)/di(j)
 
        ss0=h2t(j)
 
        r1=ai(j)*ct(j)/bi(j)
        r2=0.5*ai(j)*d1ti(j)/(bi(j)*di(j))
        r3=0.5*a1ti(j)/bi(j)
        r4=phi1ti(j)*ai(j)/bi(j)
        s1=h3*(r1+r2+r3+r4)
 
        r1=-0.5*d1ni(j)/di(j)
        r2=-b1ni(j)/bi(j)
        r3=0.5*a1ni(j)/ai(j)
        r4=-phi1ni(j)
        s2=h2*(r1+r2+r3+r4)
 
        r1=ct(j)
        r2=0.5*d1ti(j)/di(j)
        r3=-0.5*b1ti(j)/bi(j)
        r4=a1ti(j)/ai(j)
        r5=phi1ti(j)
        s3=h1t(j)*(r1+r2+r3+r4+r5)
 
        r1=-0.5*bi(j)*d1ni(j)/(ai(j)*di(j))
        r2=-0.5*b1ni(j)/ai(j)
        r3=-bi(j)*phi1ni(j)/ai(j)
        s4=r1+r2+r3
 
        r1=-bi(j)*hai(j)/(ai(j)*ai(j))
        r2=-h2*hbi(j)/bi(j)
        r3=2.*hci(j)*h1t(j)/ai(j)
        s5=psii(j)*psii(j)*fh*(r1+r2+r3)
 
        s6=bi(j)*psii(j)*psii(j)*fh3*trk/(ai(j)*ai(j))
 
        res(j)=ss0+s1+s2+s3+s4+s5+s6
   10 continue
      return
      end
 
 
 
      subroutine check(iflag)
c*********************************************************************
c  
c  check tolerance for convergence of Newton-Raphson
c
c*********************************************************************
c
      include 'param.h'
      common /residblk/ res(jm)
      call residual
      rmax=-1.
      do 20 j=2,nthtamax
        r=abs(res(j))
        rmax=max(r,rmax)
   20 continue
      epsilon=1.0e-8
      if (rmax.lt.epsilon) then
      iflag=1.
      else
      iflag=0.
      endif
      return
      end
 
 
 
      subroutine interpl3
c***********************************************************************
c
c  routine to set up for cubic spline interpolation
c
c***********************************************************************
c
      include 'param.h'
      common /metrcblk/ a(im,jm),b(im,jm),d(im,jm)
      common /excrvblk/ ha(im,jm),hb(im,jm),hc(im,jm),hd(im,jm)
      common /cnfmlblk/ psim4(im,jm),psi(im,jm),phi1n(im,jm),
     &                  phi2n(im,jm),phi1t(im,jm),
     &                  phi2t(im,jm),phi2m(im,jm)
      common /gdervblk/ a1n(im,jm),b1n(im,jm),d1t(im,jm),
     &                  d1n(im,jm),a1t(im,jm),b1t(im,jm), 
     &                  b2n(im,jm),d2n(im,jm),a2n(im,jm),
     &                  a2t(im,jm),b2t(im,jm),d2m(im,jm),
     &                  d2t(im,jm),a2m(im,jm),b2m(im,jm)
      common /hdervblk/ ha1n(im,jm),ha1t(im,jm),hb1n(im,jm),
     &                  hb1t(im,jm),hc1n(im,jm),hc1t(im,jm),
     &                  hd1n(im,jm),hd1t(im,jm)
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
 
      call spline3( a, ai,1., 1., 1., 1.)
      call spline3( b, bi,1., 1., 1., 1.)
      call spline3( d, di,1., 1., 1., 1.)
      call spline3(ha,hai,0., 1., 1., 1.)
      call spline3(hb,hbi,0., 1., 1., 1.)
      call spline3(hd,hdi,0., 1., 1., 1.)
      call spline3(hc,hci,0.,-1.,-1.,-1.)
 
      call spline3( a1n, a1ni,0., 1., 1.,-1.)
      call spline3( b1n, b1ni,0., 1., 1.,-1.)
      call spline3( d1n, d1ni,0., 1., 1.,-1.)
      call spline3(ha1n,ha1ni,0., 1., 1.,-1.)
      call spline3(hb1n,hb1ni,0., 1., 1.,-1.)
      call spline3(hd1n,hd1ni,0., 1., 1.,-1.)
      call spline3(hc1n,hc1ni,0.,-1.,-1., 1.)
 
      call spline3( a1t, a1ti,0.,-1.,-1., 1.)
      call spline3( b1t, b1ti,0.,-1.,-1., 1.)
      call spline3( d1t, d1ti,0.,-1.,-1., 1.)
      call spline3(ha1t,ha1ti,0.,-1.,-1., 1.)
      call spline3(hb1t,hb1ti,0.,-1.,-1., 1.)
      call spline3(hd1t,hd1ti,0.,-1.,-1., 1.)
      call spline3(hc1t,hc1ti,0., 1., 1.,-1.)
 
      call spline3(a2n,a2ni,0., 1., 1., 1.)
      call spline3(b2n,b2ni,0., 1., 1., 1.)
      call spline3(d2n,d2ni,0., 1., 1., 1.)
 
      call spline3(a2t,a2ti,0., 1., 1., 1.)
      call spline3(b2t,b2ti,0., 1., 1., 1.)
      call spline3(d2t,d2ti,0., 1., 1., 1.)
 
      call spline3(a2m,a2mi,0.,-1.,-1.,-1.)
      call spline3(b2m,b2mi,0.,-1.,-1.,-1.)
      call spline3(d2m,d2mi,0.,-1.,-1.,-1.)
 
      p1=psi(netamax,1)
      p2=phi1n(netamax,1)
      p3=phi1t(netamax,1)
      p4=phi2n(netamax,1)
      p5=phi2t(netamax,1)
      p6=phi2m(netamax,1)
      call spline3(psi  ,psii  ,p1, 1., 1., 1.)
      call spline3(phi1n,phi1ni,p2, 1., 1.,-1.)
      call spline3(phi1t,phi1ti,p3,-1.,-1., 1.)
      call spline3(phi2n,phi2ni,p4, 1., 1., 1.)
      call spline3(phi2t,phi2ti,p5, 1., 1., 1.)
      call spline3(phi2m,phi2mi,p6,-1.,-1.,-1.)
      return
      end
 
 
 
      subroutine spline3(u,ui,bnd,sa,se,st)
c**********************************************************************
c
c  a cubic spline algorithm for interpolations
c
c**********************************************************************
c
      include 'param.h'
      dimension u(im,jm),ui(jm)
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /dervablk/ ub(-1:im+2,-1:jm+2)
      common /gridblck/ eta(im),theta(jm)
      real sbnd(jm)
      real sa,se,st
 
      do 1 j=1,nthtamax
        sbnd(j)=bnd
    1 continue
 
      call boundry2(u,sbnd,sa,se,st)
 
      do 10 j=1,nthtamax
        i=int(h(j)/deta)+1
        ip1=i+1
        im1=i-1
        ip2=i+2
 
c  form radial derivatives at i and i+1
        ub1ni=(ub(ip1,j)-ub(im1,j))/(2.*deta)
        ub1nip1=(ub(ip2,j)-ub(i,j))/(2.*deta)
 
c  form cubic spline coefficients at i 
        sp=abs((ub(ip1,j)-ub(i,j))/deta)
        sm=abs((ub(i,j)-ub(im1,j))/deta)
        rmin=3.*min(sp,sm)
        fmin=min(0.,ub1ni)
        fmax=max(0.,ub1ni)
        filtmin=min(fmax,rmin)
        filtmax=max(fmin,-rmin)
        sigma=sign(1.,ub1ni)
        c2i=filtmin
        if (sigma.lt.0.) c2i=filtmax
 
        sp=abs((ub(ip2,j)-ub(ip1,j))/deta)
        sm=abs((ub(ip1,j)-ub(i,j))/deta)
        rmin=3.*min(sp,sm)
        fmin=min(0.,ub1nip1)
        fmax=max(0.,ub1nip1)
        filtmin=min(fmax,rmin)
        filtmax=max(fmin,-rmin)
        sigma=sign(1.,ub1nip1)
        c2ip1=filtmin
        if (sigma.lt.0.) c2ip1=filtmax
 
        sp=(ub(ip1,j)-ub(i,j))/deta
        c3i=(3.*sp-c2ip1-2.*c2i)/deta
        c4i=(-2.*sp+c2ip1+c2i)/(deta*deta)
 
        d=h(j)-eta(i)
        d2=d*d
        d3=d2*d
 
        ui(j)=ub(i,j)+d*c2i+d2*c3i+d3*c4i
   10 continue
      return
      end
 
 
 
      subroutine ahmass
c**********************************************************************
c
c  compute mass of apparent horizon
c
c**********************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /admmablk/ adm,igdet
      real arg(2*jmax-1)

      call diffh
      call interpl3

      do 10 j=1,jmax-1
        r1=(psii(j)**4.)*sqrt(bi(j)*di(j))/adm**2
        r2=sqrt(1.+(ai(j)/bi(j))*h1t(j)*h1t(j))
        arg(j)=r1*r2*sn(j)

        r3=(psii(j+1)**4.)*sqrt(bi(j+1)*di(j+1))/adm**2
        r4=sqrt(1.+(ai(j+1)/bi(j+1))*h1t(j+1)*h1t(j+1))
        arg(2*jmax-j)=r3*r4*sn(j+1)
   10 continue
      r1=(psii(jmax)**4.)*sqrt(bi(jmax)*di(jmax))/adm**2
      r2=sqrt(1.+(ai(jmax)/bi(jmax))*h1t(jmax)*h1t(jmax))
      arg(jmax)=r1*r2*sn(jmax)

      call simpson(arg,dtheta,2*jmax-1,ans)
c     ans=ans-2.*(.5*arg(1)*dtheta/2.)
      aharea=2.*pi*ans
      ahm=0.25*sqrt(aharea/pi)
      return
      end
 
 
 
      subroutine circumf(ans)
c***********************************************************************
c
c  compute polar and equatorial circumference of trapped surface
c
c***********************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
      common /admmablk/ adm,igdet
      real arg(2*jmax-1)
 
      call diffh
      call interpl3
 
c  equatorial circumference of trapped surface
      circe=2.*pi*psii(nthtamax)*psii(nthtamax)*sqrt(di(nthtamax))
      circe=circe/adm
 
c  polar circumference of trapped surface
      do 10 j=1,jmax-1
        arg(j)=psii(j)*psii(j)*sqrt(ai(j)*h1t(j)*h1t(j)+bi(j))/adm
        arg(2*jmax-j)=psii(j+1)*psii(j+1)*sqrt(ai(j+1)*h1t(j+1)
     $               *h1t(j+1)+bi(j+1))/adm
   10 continue
      arg(jmax)=psii(jmax)*psii(jmax)*sqrt(ai(jmax)*h1t(jmax)
     $         *h1t(jmax)+bi(jmax))/adm
      call simpson(arg,dtheta,2*jmax-1,ans)
c     ans=ans-2.*(.5*arg(1)*dtheta/2.)
      circp=2.*ans
 
      ans=circp/circe
      return
      end
 
 
 

      subroutine boundry2(ux,bnd,sa,se,st)
c***********************************************************************
c  
c  define boundary values outside computational domain, used for
c  interpolating data when locating trapped surfaces
c
c***********************************************************************
c
      include 'param.h'
      dimension ux(im,jm),bnd(jm)
      common /dervablk/ ub(-1:im+2,-1:jm+2)
 
      do 5 j=1,nthtamax
      do 5 i=1,netamax
    5 ub(i,j)=ux(i,j)
 
      do 10 i=1,netamax
        ub(i,0)=sa*ux(i,3)
        ub(i,-1)=sa*ux(i,4)
        ub(i,ntp1)=se*ux(i,nthtamax)
        ub(i,ntp2)=se*ux(i,ntm1)
   10 continue
 
      do 20 j=1,nthtamax
        ub(0,j)=st*ux(2,j)
        ub(-1,j)=st*ux(3,j)
   20 continue
 
      ub(0,0)=ux(2,3)
      ub(0,-1)=ux(2,4)
      ub(-1,0)=ux(3,3)
      ub(-1,-1)=ux(3,4)
      ub(0,ntp1)=ux(2,nthtamax)
      ub(0,ntp2)=ux(2,ntm1)
      ub(-1,ntp1)=ux(3,nthtamax)
      ub(-1,ntp2)=ux(3,ntm1)

      do 30 j=1,nthtamax
        ub(nep1,j)=bnd(j)
        ub(nep2,j)=bnd(j)
   30 continue
 
      ub(nep1, 0)=sa*bnd(3)
      ub(nep1,-1)=sa*bnd(4)
      ub(nep2, 0)=sa*bnd(3)
      ub(nep2,-1)=sa*bnd(4)
 
      ub(nep1,ntp1)=se*bnd(nthtamax)
      ub(nep1,ntp2)=se*bnd(ntm1)
      ub(nep2,ntp1)=se*bnd(nthtamax)
      ub(nep2,ntp2)=se*bnd(ntm1)
      return
      end
 
 
 
      subroutine brinvnp(ric11,ric22,ric12,ric33)
c**********************************************************************
c
c  compute Bel-Robinson vector, Riemann invariants and Newman-Penrose
c  scalars from "E" and "B" construction of Weyl tensor
c
c**********************************************************************
c
      include 'param.h'
      common /metrcblk/ a(im,jm),b(im,jm),d(im,jm)
      common /excrvblk/ ha(im,jm),hb(im,jm),hc(im,jm),hd(im,jm)
      common /gdervblk/ a1n(im,jm),b1n(im,jm),d1t(im,jm),
     &                  d1n(im,jm),a1t(im,jm),b1t(im,jm), 
     &                  b2n(im,jm),d2n(im,jm),a2n(im,jm),
     &                  a2t(im,jm),b2t(im,jm),d2m(im,jm),
     &                  d2t(im,jm),a2m(im,jm),b2m(im,jm)
      common /hdervblk/ ha1n(im,jm),ha1t(im,jm),hb1n(im,jm),
     &                  hb1t(im,jm),hc1n(im,jm),hc1t(im,jm),
     &                  hd1n(im,jm),hd1t(im,jm)
      common /cnfmlblk/ psim4(im,jm),psi(im,jm),phi1n(im,jm),
     &                  phi2n(im,jm),phi1t(im,jm),
     &                  phi2t(im,jm),phi2m(im,jm)
      common /gridblck/ eta(im),theta(jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /radiablk/ bre(im,jm),brt(im,jm),brn(im,jm),ssj(im,jm),
     &                  crvi(im,jm),crvj(im,jm),ssi(im,jm),
     &                  y13(im,jm),y23(im,jm),psi0(im,jm),psi1(im,jm),
     &                  psi2(im,jm),psi3(im,jm),psi4(im,jm)
      common /admmablk/ adm,igdet
      dimension ric11(im,jm),ric22(im,jm),ric12(im,jm),ric33(im,jm)

      do 31 j=1,nthtamax
      do 11 i=1,netamax-1
        fac=1./sqrt(a(i,j)*b(i,j)*d(i,j))

c  construct Christoffel symbols 
c  (g133 and g233 are divided by sin(theta)**2, like ric33)
c--------------------------------------------------------------------
        g111=0.5*(a1n(i,j)/a(i,j)+phi1n(i,j))
        g112=0.5*(a1t(i,j)/a(i,j)+phi1t(i,j))
        g122=-0.5*(b1n(i,j)+b(i,j)*phi1n(i,j))/a(i,j)
        g211=-0.5*(a1t(i,j)+a(i,j)*phi1t(i,j))/b(i,j)
        g212=0.5*(b1n(i,j)/b(i,j)+phi1n(i,j))
        g222=0.5*(b1t(i,j)/b(i,j)+phi1t(i,j))
        g313=0.5*(d1n(i,j)/d(i,j)+phi1n(i,j))
        g323=0.5*(d1t(i,j)/d(i,j)+phi1t(i,j))+ct(j)
        g133=-0.5*(d1n(i,j)+d(i,j)*phi1n(i,j))/a(i,j)
        g233=-0.5*(d1t(i,j)+d(i,j)*phi1t(i,j))/b(i,j)
     &            -d(i,j)*ct(j)/b(i,j)

c  construct electric and magnetic components of the Riemann tensor
c  (a factor of psi^4 is taken out of all the electric components
c   and an additional sin(theta)^2 is taken out of e33. 
c   a factor of sin(theta)*psi^2 is taken out of all the 
c   magnetic components)
c---------------------------------------------------------------------
        e11=-ric11(i,j)*psim4(i,j)
     &           -ha(i,j)*(hb(i,j)/b(i,j)+hd(i,j)/d(i,j))
     &           +hc(i,j)*hc(i,j)/b(i,j)
        e12=-ric12(i,j)*psim4(i,j)-hc(i,j)*hd(i,j)/d(i,j)
        e22=-ric22(i,j)*psim4(i,j)
     &           -hb(i,j)*(ha(i,j)/a(i,j)+hd(i,j)/d(i,j))
     &           +hc(i,j)*hc(i,j)/a(i,j)
        e33=-ric33(i,j)*psim4(i,j)
     &           -hd(i,j)*(ha(i,j)/a(i,j)+hb(i,j)/b(i,j))

        r1=d(i,j)*(hc1n(i,j)+hc(i,j)*phi1n(i,j)
     &            -ha1t(i,j)-ha(i,j)*phi1t(i,j))
        r2=a(i,j)*(hd1t(i,j)+hd(i,j)*phi1t(i,j)+2.*hd(i,j)*ct(j))
        r3=d(i,j)*(-ha(i,j)*g112+hb(i,j)*g211+hc(i,j)*(g111-g212))
        r4=a(i,j)*(-hc(i,j)*g133-hb(i,j)*g233+hd(i,j)*g323)
        b13=-0.5*fac*(r1+r2-r3-r4)

        r1=d(i,j)*(-hc1t(i,j)-hc(i,j)*phi1t(i,j)
     &             +hb1n(i,j)+hb(i,j)*phi1n(i,j))
        r2=-b(i,j)*(hd1n(i,j)+hd(i,j)*phi1n(i,j))
        r3=d(i,j)*(-ha(i,j)*g122+hb(i,j)*g212+hc(i,j)*(g112-g222))
        r4=b(i,j)*( ha(i,j)*g133+hc(i,j)*g233-hd(i,j)*g313)
        b23=-0.5*fac*(r1+r2-r3-r4)
 
c  construct radiation quantities from E & B fields
c----------------------------------------------------------------------
        r1=-b13*e12/a(i,j)
        r2=-b23*e22/b(i,j)
        r3= b23*e33/d(i,j)
        bre(i,j)=fac*(r1+r2+r3)*psim4(i,j)

        r1= b13*e11/a(i,j)
        r2= b23*e12/b(i,j)
        r3=-b13*e33/d(i,j)
        brt(i,j)=fac*(r1+r2+r3)*psim4(i,j)

        r1=sqrt(bre(i,j)*bre(i,j)*a(i,j)+brt(i,j)*brt(i,j)*b(i,j))
        brn(i,j)=psi(i,j)*psi(i,j)*r1

        bre(i,j)=bre(i,j)*adm**5
        brt(i,j)=brt(i,j)*adm**5
        brn(i,j)=brn(i,j)*adm**4

        psi4(i,j)=.25*(e33/d(i,j)-e22/b(i,j))
     $           +.5*b23/(psi(i,j)**2*sqrt(b(i,j)*d(i,j)))
        psi0(i,j)=.25*(e33/d(i,j)-e22/b(i,j))
     $           -.5*b23/(psi(i,j)**2*sqrt(b(i,j)*d(i,j)))
        psi2(i,j)=-.5*e11/a(i,j)
        psi4(i,j)=psi4(i,j)*adm**2
        psi0(i,j)=psi0(i,j)*adm**2
        psi2(i,j)=psi2(i,j)*adm**2

        e2=(e11/a(i,j))**2.+2.*e12*e12/(a(i,j)*b(i,j))
     &    +(e22/b(i,j))**2.+(e33/d(i,j))**2.
        b2=2.*psim4(i,j)*(b13*b13/a(i,j)+b23*b23/b(i,j))/d(i,j)
        crvi(i,j)=4.*(e2-b2)
        crvi(i,j)=crvi(i,j)*adm**4
 
        ss1=(e11/a(i,j))**3.
        ss2=(e22/b(i,j))**3.
        ss3=(e33/d(i,j))**3.
        ss4=3.*e11*e12*e12/(a(i,j)*a(i,j)*b(i,j))
        ss5=3.*e12*e12*e22/(a(i,j)*b(i,j)*b(i,j))
        r1=ss1+ss2+ss3+ss4+ss5
 
        ss6 =-3.*b13*b13*e11/(a(i,j)*a(i,j))
        ss7 =-6.*b13*b23*e12/(a(i,j)*b(i,j))
        ss8 =-3.*b23*b23*e22/(b(i,j)*b(i,j))
        ss9 =-3.*b13*b13*e33/(a(i,j)*d(i,j))
        ss10=-3.*b23*b23*e33/(b(i,j)*d(i,j))
        r2=psim4(i,j)*(ss6+ss7+ss8+ss9+ss10)/d(i,j)

        crvj(i,j)=8.*(r1+r2)
        crvj(i,j)=crvj(i,j)*adm**6
 11   continue
 31   continue

c  apply boundary condition at equator
      j=nthtamax+1
      do 51 i=1,netamax-1
        bre(i,j)=bre(i,j-1)
        brt(i,j)=brt(i,j-1)
        brn(i,j)=brn(i,j-1)
        crvi(i,j)=crvi(i,j-1)
        crvj(i,j)=crvj(i,j-1)
        psi4(i,j)=psi4(i,j-1)
        psi0(i,j)=psi0(i,j-1)
        psi2(i,j)=psi2(i,j-1)
 51   continue

c  spherically symmetric invariants
      do 37 j=1,nthtamax+1
      do 17 i=1,netamax-1
        e=eta(i)
        r=2*adm*(cosh(0.5*e)**2.)
        ssi(i,j)=(24.*adm*adm)/(r**(6.0))
        ssj(i,j)=(48.*adm*adm*adm)/(r**(9.0))

        ssi(i,j)=ssi(i,j)*adm**4
        ssj(i,j)=ssj(i,j)*adm**6
 17   continue
 37   continue

c compute Newman-Penrose scalars from direct construction
c (see David Bernstein's thesis reference [33])
c----------------------------------------------------------------------
      call weyl
      return
      end
 




      subroutine york
c**********************************************************************
c
c  compute nonvanishing components of Cotton-York tensor
c  (see David Bernstein's thesis reference [33])
c
c**********************************************************************
c
      include 'param.h'
      common /metrcblk/ a(im,jm),b(im,jm),d(im,jm)
      common /gdervblk/ a1n(im,jm),b1n(im,jm),d1t(im,jm),
     &                  d1n(im,jm),a1t(im,jm),b1t(im,jm), 
     &                  b2n(im,jm),d2n(im,jm),a2n(im,jm),
     &                  a2t(im,jm),b2t(im,jm),d2m(im,jm),
     &                  d2t(im,jm),a2m(im,jm),b2m(im,jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /radiablk/ bre(im,jm),brt(im,jm),brn(im,jm),ssj(im,jm),
     &                  crvi(im,jm),crvj(im,jm),ssi(im,jm),
     &                  y13(im,jm),y23(im,jm),psi0(im,jm),psi1(im,jm),
     &                  psi2(im,jm),psi3(im,jm),psi4(im,jm)
      real a2t1n(im,jm),a3t(im,jm),b2n1t(im,jm),b3n(im,jm),
     &     d2t1n(im,jm),d2n1t(im,jm),d3t(im,jm),d3n(im,jm),
     &     b2t1n(im,jm),b3t(im,jm)
      real xbnd(jm)

c  construct yorkcrv

      do 22 j=1,nthtamax
 22   xbnd(j)=0.

      call frstdrvs(a2t,a2t1n,a3t,1.,xbnd,netamax,nthtamax)
      call frstdrvs(b2t,b2t1n,b3t,1.,xbnd,netamax,nthtamax)
      call frstdrvs(d2t,d2t1n,d3t,1.,xbnd,netamax,nthtamax)
      call frstdrvs(b2n,b3n,b2n1t,1.,xbnd,netamax,nthtamax)
      call frstdrvs(d2n,d3n,d2n1t,1.,xbnd,netamax,nthtamax)

      do 34 j=1,nthtamax
      do 24 i=1,netamax-1
      fac=(a(i,j)*b(i,j)*d(i,j)*sn(j)**2.)**(1./3.)
 
      r1 =-0.5*d2t(i,j)*ct(j)/(b(i,j)*d(i,j))
      r2 =0.5*d1t(i,j)*d1t(i,j)*ct(j)/(b(i,j)*d(i,j)*d(i,j))
      r3 =0.5*b1t(i,j)*d1t(i,j)*ct(j)/(b(i,j)*b(i,j)*d(i,j))
      r4 =-0.5*a1t(i,j)*d1t(i,j)*ct(j)/(a(i,j)*b(i,j)*d(i,j))
      r5 =-0.5*d2n(i,j)*ct(j)/(a(i,j)*d(i,j))
      r6 =0.5*d1n(i,j)*d1n(i,j)*ct(j)/(a(i,j)*d(i,j)*d(i,j))
      r7 =-0.25*b1n(i,j)*d1n(i,j)*ct(j)/(a(i,j)*b(i,j)*d(i,j))
      r8 =0.25*a1n(i,j)*d1n(i,j)*ct(j)/(a(i,j)*a(i,j)*d(i,j))
      r9 =0.25*b2t(i,j)*ct(j)/(b(i,j)*b(i,j))
      r10=-0.5*b1t(i,j)*b1t(i,j)*ct(j)/(b(i,j)*b(i,j)*b(i,j))
      r11=0.5*b2n(i,j)*ct(j)/(a(i,j)*b(i,j))
      r12=-0.25*b1n(i,j)*b1n(i,j)*ct(j)/(a(i,j)*b(i,j)*b(i,j))
      r13=-0.25*a1n(i,j)*b1n(i,j)*ct(j)/(a(i,j)*a(i,j)*b(i,j))
      r14=0.25*a2t(i,j)*ct(j)/(a(i,j)*b(i,j))
      r15=0.5*d1t(i,j)*ct(j)*ct(j)/(b(i,j)*d(i,j))
      r16=-0.25*b1t(i,j)*ct(j)*ct(j)/(b(i,j)*b(i,j))
      r17=-0.25*a1t(i,j)*ct(j)*ct(j)/(a(i,j)*b(i,j))
      r18=-0.25*d3t(i,j)/(b(i,j)*d(i,j))
      r19=0.5*d1t(i,j)*d2t(i,j)/(b(i,j)*d(i,j)*d(i,j))
      r20=0.375*b1t(i,j)*d2t(i,j)/(b(i,j)*b(i,j)*d(i,j))
      r21=-0.125*a1t(i,j)*d2t(i,j)/(a(i,j)*b(i,j)*d(i,j))
      r22=-0.25*d1t(i,j)*d1t(i,j)*d1t(i,j)/(b(i,j)*d(i,j)**3)
      r23=-0.25*b1t(i,j)*d1t(i,j)*d1t(i,j)/(b(i,j)**2.*d(i,j)**2)
      r24=-0.25*d1n(i,j)*d1n(i,j)*d1t(i,j)/(a(i,j)*d(i,j)**3)
      r25=0.125*b2t(i,j)*d1t(i,j)/(b(i,j)*b(i,j)*d(i,j))
      r26=-0.25*b1t(i,j)*b1t(i,j)*d1t(i,j)/(d(i,j)*b(i,j)**3)
      r27=0.25*b2n(i,j)*d1t(i,j)/(a(i,j)*b(i,j)*d(i,j))
      r28=-0.125*b1n(i,j)*b1n(i,j)*d1t(i,j)/(a(i,j)*d(i,j)*b(i,j)**2)
      r29=-0.125*a1n(i,j)*b1n(i,j)*d1t(i,j)/(b(i,j)*d(i,j)*a(i,j)**2)
      r30=0.125*a2t(i,j)*d1t(i,j)/(a(i,j)*b(i,j)*d(i,j))
      r31=0.5*d1t(i,j)/(b(i,j)*d(i,j))
      r32=-0.25*d2n1t(i,j)/(a(i,j)*d(i,j))
      r33=0.25*a1t(i,j)*d2n(i,j)/(a(i,j)*a(i,j)*d(i,j))
      r34=0.5*d1n(i,j)*d2m(i,j)/(a(i,j)*d(i,j)*d(i,j))
      r35=-0.125*b1n(i,j)*d2m(i,j)/(a(i,j)*b(i,j)*d(i,j))
      r36=0.125*a1n(i,j)*d2m(i,j)/(a(i,j)*a(i,j)*d(i,j))
      r37=-0.25*a1t(i,j)*d1n(i,j)*d1n(i,j)/(a(i,j)**2*d(i,j)**2)
      r38=0.125*b1n(i,j)*b1t(i,j)*d1n(i,j)/(a(i,j)*d(i,j)*b(i,j)**2)
      r39=-0.125*b2m(i,j)*d1n(i,j)/(a(i,j)*d(i,j)*b(i,j))
      r40=0.125*a1t(i,j)*b1n(i,j)*d1n(i,j)/(b(i,j)*d(i,j)*a(i,j)**2)
      r41=-0.25*a1n(i,j)*a1t(i,j)*d1n(i,j)/(d(i,j)*a(i,j)**3)
      r42=0.125*a2m(i,j)*d1n(i,j)/(a(i,j)*a(i,j)*d(i,j))
      r43=-0.125*a1t(i,j)*b2t(i,j)/(a(i,j)*b(i,j)*b(i,j))
      r44=0.25*a1t(i,j)*b1t(i,j)*b1t(i,j)/(a(i,j)*b(i,j)**3)
      r45=-0.25*b2n(i,j)*b1t(i,j)/(a(i,j)*b(i,j)*b(i,j))
      r46=0.25*b1n(i,j)*b1n(i,j)*b1t(i,j)/(a(i,j)*b(i,j)**3)
      r47=0.125*a1n(i,j)*b1n(i,j)*b1t(i,j)/(a(i,j)**2.*b(i,j)**2)
      r48=-0.375*a2t(i,j)*b1t(i,j)/(a(i,j)*b(i,j)*b(i,j))
      r49=0.25*a1t(i,j)*a1t(i,j)*b1t(i,j)/(a(i,j)**2*b(i,j)**2)
      r50=-0.75*b1t(i,j)/(b(i,j)*b(i,j))
      r51=0.25*b2n1t(i,j)/(a(i,j)*b(i,j))
      r52=-0.25*a1t(i,j)*b2n(i,j)/(a(i,j)*a(i,j)*b(i,j))
      r53=-0.25*b1n(i,j)*b2m(i,j)/(a(i,j)*b(i,j)*b(i,j))
      r54=-0.125*a1n(i,j)*b2m(i,j)/(a(i,j)*a(i,j)*b(i,j))
      r55=0.125*a1t(i,j)*b1n(i,j)*b1n(i,j)/(a(i,j)**2*b(i,j)**2)
      r56=0.25*a1n(i,j)*a1t(i,j)*b1n(i,j)/(b(i,j)*a(i,j)**3)
      r57=-0.125*a2m(i,j)*b1n(i,j)/(a(i,j)*a(i,j)*b(i,j))
      r58=0.25*a3t(i,j)/(a(i,j)*b(i,j))
      r59=-0.5*a1t(i,j)*a2t(i,j)/(a(i,j)*a(i,j)*b(i,j))
      r60=0.25*((a1t(i,j)/a(i,j))**3)/b(i,j)
      r61=0.25*a1t(i,j)/(a(i,j)*b(i,j))
 
      y13(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &        +r10+r11+r12+r13+r14+r15+r16+r17+r18+r19
     &        +r20+r21+r22+r23+r24+r25+r26+r27+r28+r29
     &        +r30+r31+r32+r33+r34+r35+r36+r37+r38+r39
     &        +r40+r41+r42+r43+r44+r45+r46+r47+r48+r49
     &        +r50+r51+r52+r53+r54+r55+r56+r57+r58+r59
     &        +r60+r61
 
      y13(i,j)=y13(i,j)*fac
 
      t1 =-0.25*b1t(i,j)*d1n(i,j)*ct(j)/(b(i,j)*b(i,j)*d(i,j))
      t2 =0.25*a1t(i,j)*d1n(i,j)*ct(j)/(a(i,j)*b(i,j)*d(i,j))
      t3 =0.5*b1n(i,j)*b1t(i,j)*ct(j)/(b(i,j)**3)
      t4 =-0.25*b2m(i,j)*ct(j)/(b(i,j)*b(i,j))
      t5 =-0.25*a1t(i,j)*b1n(i,j)*ct(j)/(a(i,j)*b(i,j)*b(i,j))
      t6 =-0.25*a1n(i,j)*a1t(i,j)*ct(j)/(a(i,j)*a(i,j)*b(i,j))
      t7 =0.25*a2m(i,j)*ct(j)/(a(i,j)*b(i,j))
      t8 =-0.5*d1n(i,j)*ct(j)*ct(j)/(b(i,j)*d(i,j))
      t9 =0.5*b1n(i,j)*ct(j)*ct(j)/(b(i,j)*b(i,j))
      t10=-0.25*b1n(i,j)*d2t(i,j)/(b(i,j)*b(i,j)*d(i,j))
      t11=0.25*d1n(i,j)*d1t(i,j)*d1t(i,j)/(b(i,j)*d(i,j)**3)
      t12=0.25*b1n(i,j)*d1t(i,j)*d1t(I,j)/(b(i,j)**2*d(i,j)**2)
      t13=-0.5*d2m(i,j)*d1t(i,j)/(b(i,j)*d(i,j)*d(i,j))
      t14=0.25*b1n(i,j)*b1t(i,j)*d1t(i,j)/(d(i,j)*b(i,j)**3)
      t15=-0.125*b2m(i,j)*d1t(i,j)/(b(i,j)*b(i,j)*d(i,j))
      t16=-0.125*a1t(i,j)*b1n(i,j)*d1t(i,j)/(a(i,j)*d(i,j)*b(i,j)**2)
      t17=-0.125*a1n(i,j)*a1t(i,j)*d1t(i,j)/(b(i,j)*d(i,j)*a(i,j)**2)
      t18=0.125*a2m(i,j)*d1t(i,j)/(a(i,j)*b(i,j)*d(i,j))
      t19=0.25*d3n(i,j)/(a(i,j)*d(i,j))
      t20=-0.5*d1n(i,j)*d2n(i,j)/(a(i,j)*d(i,j)*d(i,j))
      t21=0.125*b1n(i,j)*d2n(i,j)/(a(i,j)*b(i,j)*d(i,j))
      t22=-0.375*a1n(i,j)*d2n(i,j)/(a(i,j)*a(i,j)*d(i,j))
      t23=0.25*d2t1n(i,j)/(b(i,j)*d(i,j))
      t24=-0.125*b1t(i,j)*d2m(i,j)/(b(i,j)*b(i,j)*d(i,j))
      t25=0.125*a1t(i,j)*d2m(i,j)/(a(i,j)*b(i,j)*d(i,j))
      t26=0.25*d1n(i,j)*d1n(i,j)*d1n(i,j)/(a(i,j)*d(i,j)**3)
      t27=0.25*a1n(i,j)*d1n(i,j)*d1n(i,j)/(a(i,j)**2.*d(i,j)**2)
      t28=0.125*a1t(i,j)*b1t(i,j)*d1n(i,j)/(a(i,j)*d(i,j)*b(i,j)**2)
      t29=-0.125*b2n(i,j)*d1n(i,j)/(a(i,j)*b(i,j)*d(i,j))
      t30=-0.25*a2t(i,j)*d1n(i,j)/(a(i,j)*b(i,j)*d(i,j))
      t31=0.125*a1t(i,j)*a1t(i,j)*d1n(i,j)/(b(i,j)*d(i,j)*a(i,j)**2)
      t32=-0.5*d1n(i,j)/(b(i,j)*d(i,j))
      t33=-0.125*a2n(i,j)*d1n(i,j)/(a(i,j)*a(i,j)*d(i,j))
      t34=0.25*a1n(i,j)*a1n(i,j)*d1n(i,j)/(d(i,j)*a(i,j)**3)
      t35=-0.25*a1t(i,j)*b1n(i,j)*d1t(i,j)/(a(i,j)*b(i,j)**3)
      t36=-0.125*a1n(i,j)*a1t(i,j)*b1t(i,j)/(a(i,j)**2*b(i,j)**2)
      t37=0.125*a2m(i,j)*b1t(i,j)/(a(i,j)*b(i,j)*b(i,j))
      t38=-0.25*b3n(i,j)/(a(i,j)*b(i,j))
      t39=0.5*b1n(i,j)*b2n(i,j)/(a(i,j)*b(i,j)*b(i,j))
      t40=0.375*a1n(i,j)*b2n(i,j)/(a(i,j)*a(i,j)*b(i,j))
      t41=0.125*a1t(i,j)*b2m(i,j)/(a(i,j)*b(i,j)*b(i,j))
      t42=-0.25*b1n(i,j)*b1n(i,j)*b1n(i,j)/(a(i,j)*b(i,j)**3)
      t43=-0.25*a1n(i,j)*b1n(i,j)*b1n(i,j)/(a(i,j)**2*b(i,j)**2)
      t44=0.125*a2n(i,j)*b1n(i,j)/(a(i,j)*a(i,j)*b(i,j))
      t45=-0.25*a1n(i,j)*a1n(i,j)*b1n(i,j)/(b(i,j)*a(i,j)**3)
      t46=0.25*a2t(i,j)*b1n(i,j)/(a(i,j)*b(i,j)*b(i,j))
      t47=-0.125*a1t(i,j)*a1t(i,j)*b1n(i,j)/(a(i,j)**2*b(i,j)**2)
      t48=0.5*b1n(i,j)/(b(i,j)*b(i,j))
      t49=0.25*a1n(i,j)*a2t(i,j)/(a(i,j)*a(i,j)*b(i,j))
      t50=-0.25*a1n(i,j)*a1t(i,j)*a1t(i,j)/(b(i,j)*a(i,j)**3)
      t51=0.25*a2m(i,j)*a1t(i,j)/(a(i,j)*a(i,j)*b(i,j))
      t52=-0.25*a2t1n(i,j)/(a(i,j)*b(i,j))
 
      y23(i,j)=    t1 +t2 +t3 +t4 +t5 +t6 +t7 +t8 +t9
     &        +t10+t11+t12+t13+t14+t15+t16+t17+t18+t19
     &        +t20+t21+t22+t23+t24+t25+t26+t27+t28+t29
     &        +t30+t31+t32+t33+t34+t35+t36+t37+t38+t39
     &        +t40+t41+t42+t43+t44+t45+t46+t47+t48+t49
     &        +t50+t51+t52
 
      y23(i,j)=y23(i,j)*fac
 24   continue
 34   continue

c  apply boundary condition at equator
      j=nthtamax+1
      do 51 i=1,netamax-1
        y13(i,j)=y13(i,j-1)
        y23(i,j)=y23(i,j-1)
 51   continue
      return
      end
 



      subroutine simpson(u,dx,nxmax,ans)
c***********************************************************************
c
c integration algorithm using simpson's formula
c
c***********************************************************************
c
      include 'param.h'
      dimension u(2*jmax-1)
      ans=0.
      do 10 k=5,nxmax-4
   10 ans=ans+u(k)
      ans=ans+(u(1)+u(nxmax  ))*(17./48.)
     &       +(u(2)+u(nxmax-1))*(59./48.)
     &       +(u(3)+u(nxmax-2))*(43./48.)
     &       +(u(4)+u(nxmax-3))*(49./48.)
      ans=ans*dx
      return
      end
 
 
 

      subroutine frstdrvs(ux,u1n,u1t,sym,bnd,ne,nt)
c**********************************************************************
c
c  routine to compute first derivatives in radial and angular 
c  directions
c
c**********************************************************************
c
      include 'param.h'
      dimension ux(im,jm),u1n(im,jm),u1t(im,jm),bnd(jm)
      common /dervablk/ ub(-1:im+2,-1:jm+2)
      real sym

      call boundary(ux,sym,bnd,netamax,nthtamax)
      do 10 j=1,nthtamax
      jp1=j+1
      jm1=j-1
      do 10 i=1,netamax
      u1n(i,j)=(ub(i+1,j)-ub(i-1,j))/(2.*deta)
      u1t(i,j)=(ub(i,jp1)-ub(i,jm1))/(2.*dtheta)
   10 continue
      return
      end



      subroutine boundary(ux,sym,bnd,ne,nt)
c**********************************************************************
c
c  define boundary values outside the computational domain
c
c**********************************************************************
c
      include 'param.h'
      common /dervablk/ ub(-1:im+2,-1:jm+2)
      dimension ux(im,jm),bnd(jm)
      do 5 j=1,nthtamax
      do 5 i=1,netamax
    5 ub(i,j)=ux(i,j)
      do 10 i=1,netamax
      ub(i,0)=sym*ux(i,2)
      ub(i,ntp1)=sym*ux(i,ntm1)
   10 continue
      do 20 j=1,nthtamax
      ub(0,j)=sym*ux(2,j)
      ub(nep1,j)=bnd(j)
   20 continue
      ub(0,0)=ux(2,2)
      ub(0,ntp1)=ux(2,ntm1)
      ub(nep1,0)=sym*bnd(2)
      ub(nep1,ntp1)=sym*bnd(ntm1)
      return
      end





      subroutine embd
c***********************************************************************
c
c  embed the 2D trapped surface into a flat 3D space and compute
c  its gaussian curvature
c
c***********************************************************************
c
      include 'param.h'
      common /horizblk/ h(jm),h1t(jm),h2t(jm),dh(jm),ahm,aharea
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /admmablk/ adm,igdet
      common /intplblk/ ai(jm),bi(jm),di(jm),hai(jm),hbi(jm),
     &                  hci(jm),hdi(jm),psii(jm),a1ni(jm),
     &                  b1ni(jm),d1ni(jm),a1ti(jm),b1ti(jm),
     &                  d1ti(jm),a2ni(jm),b2ni(jm),d2ni(jm),
     &                  a2ti(jm),b2ti(jm),d2ti(jm),a2mi(jm),
     &                  b2mi(jm),d2mi(jm),
     &                  phi1ni(jm),phi1ti(jm),phi2ni(jm),
     &                  phi2ti(jm),phi2mi(jm),ha1ni(jm),ha1ti(jm),
     &                  hb1ni(jm),hb1ti(jm),hc1ni(jm),hc1ti(jm),
     &                  hd1ni(jm),hd1ti(jm)
      common/embdblk/ z0,zz(jm),rho(jm),gauss(jm)

      call diffh
      call interpl3
      area=aharea*adm**2

c  compute gaussian curvature
      do 20 j=1,nthtamax
        dt=d1ti(j)+h1t(j)*d1ni(j)
        phit=phi1ti(j)+h1t(j)*phi1ni(j)
        at=a1ti(j)+h1t(j)*a1ni(j)
        bt=b1ti(j)+h1t(j)*b1ni(j)
        dtt=d2ti(j)+h2t(j)*d1ni(j)+2.*h1t(j)*d2mi(j)+h1t(j)**2*d2ni(j)
        phitt=phi2ti(j)+h2t(j)*phi1ni(j)+2.*h1t(j)*phi2mi(j)
     $       +h1t(j)**2*phi2ni(j)
 
        g1=2.*ai(j)*di(j)**2*h1t(j)**2+2.*bi(j)*di(j)**2
        g2=-ai(j)*di(j)**2*h1t(j)**2-bi(j)*di(j)**2
        g3=-2.*ai(j)*di(j)**2*h1t(j)*h2t(j)
        g4=ai(j)*di(j)*dt-at*di(j)**2
        g5=bi(j)*di(j)*dt-bt*di(j)**2
        g6=-2.*ai(j)*di(j)*dt*h1t(j)*h2t(j)
        g7=2.*ai(j)*di(j)*dtt-ai(j)*dt**2-at*di(j)*dt-4.*ai(j)*di(j)**2
        g8=2.*bi(j)*di(j)*dtt-bi(j)*dt**2-bt*di(j)*dt-4.*bi(j)*di(j)**2

        g9=2.*ai(j)*di(j)**2*h1t(j)**2+2.*bi(j)*di(j)**2
        g10=-4.*ai(j)*di(j)**2*h1t(j)*h2t(j)
        g11=4.*ai(j)*di(j)*dt-2.*at*di(j)**2
        g12=4.*bi(j)*di(j)*dt-2.*bt*di(j)**2
        g13=4.*di(j)**2*(ai(j)*h1t(j)**2+bi(j))**2

        r1=g1*(phitt+phit**2/4.)
     $    +g2*phit**2/2.
     $    +(g3+g4*h1t(j)**2+g5)*phit
     $    +(g6+g7*h1t(j)**2+g8)
        r2=g9*phit
     $    +g10+g11*h1t(j)**2+g12
        r3=g13*psii(j)**4
 
        gauss(j)=-(r1/r3)-ct(j)*(r2/r3)

c  rescale gaussian curvature by the area to be 4pi when spherical
        gauss(j)=gauss(j)*area
  20  continue


c  compute embedding
      sum=0.0
      do 10 j=1,nthtamax-1
        ja=j
        jb=j+1

        x1a=(psii(ja)**4)*(bi(ja)+ai(ja)*h1t(ja)**2)
        x2a=cs(ja)*sqrt(di(ja))*psii(ja)**2
        x3a=sn(ja)/(2.*sqrt(di(ja))*psii(ja)**2)
        psi1t=psii(ja)*phi1ti(ja)/4.
        psi1n=psii(ja)*phi1ni(ja)/4.
        x4a=4.*di(ja)*psii(ja)**3*(psi1t+h1t(ja)*psi1n)
     $     +psii(ja)**4*(d1ti(ja)+h1t(ja)*d1ni(ja))
        xa =x1a-(x2a+x3a*x4a)**2
        if(j.eq.1)xa=0.0

        x1b=(psii(jb)**4)*(bi(jb)+ai(jb)*h1t(jb)**2)
        x2b=cs(jb)*sqrt(di(jb))*psii(jb)**2
        x3b=sn(jb)/(2.*sqrt(di(jb))*psii(jb)**2)
        psi1t=psii(jb)*phi1ti(jb)/4.
        psi1n=psii(jb)*phi1ni(jb)/4.
        x4b=4.*di(jb)*psii(jb)**3*(psi1t+h1t(jb)*psi1n)
     $     +psii(jb)**4*(d1ti(jb)+h1t(jb)*d1ni(jb))
        xb =x1b-(x2b+x3b*x4b)**2

        if(xa.lt.-dtheta)go to 999
        if(xb.lt.-dtheta)go to 999
        if(xa.lt.0.)xa=0.
        if(xb.lt.0.)xb=0.

        funca=-sqrt(xa)
        funcb=-sqrt(xb)
        sum=sum+0.5*dtheta*(funca+funcb)
        if(j.eq.1)sum=0.25*dtheta*(funca+funcb)
        zz(jb)=sum
        rho(jb)=sn(jb)*sqrt(di(jb))*psii(jb)**2
  10  continue
      zz(1)=zz(2)
      rho(1)=sn(1)*sqrt(di(1))*psii(1)**2

c  Now rescale the figure by the sqrt(area). This will keep the
c  figure the same size, preventing any spurious growth due to
c  the longitudinal error in the evolution. It also prevents
c  growth due to the black hole gaining mass.

      z0=zz(nthtamax)
      do 11 j=1,nthtamax
        zz(j)=(zz(j)-z0)/sqrt(area)
	rho(j)=rho(j)/sqrt(area)
  11  continue
      go to 777

 999  continue

      do 888 j=1,nthtamax
        zz(j)=0.
        rho(j)=0.
 888  continue

 777  return
      end




      subroutine weyl
c***********************************************************************
c
c this subroutine constructs the Newman-Penrose scalars,
c psi0 through psi4, by direct construction 
c (see David Bernstein's thesis reference [33])
c
c***********************************************************************
c
      include 'param.h'
      common /metrcblk/ a(im,jm),b(im,jm),d(im,jm)
      common /excrvblk/ ha(im,jm),hb(im,jm),hc(im,jm),hd(im,jm)
      common /gdervblk/ a1n(im,jm),b1n(im,jm),d1t(im,jm),
     &                  d1n(im,jm),a1t(im,jm),b1t(im,jm), 
     &                  b2n(im,jm),d2n(im,jm),a2n(im,jm),
     &                  a2t(im,jm),b2t(im,jm),d2m(im,jm),
     &                  d2t(im,jm),a2m(im,jm),b2m(im,jm)
      common /hdervblk/ ha1n(im,jm),ha1t(im,jm),hb1n(im,jm),
     &                  hb1t(im,jm),hc1n(im,jm),hc1t(im,jm),
     &                  hd1n(im,jm),hd1t(im,jm)
      common /cnfmlblk/ psim4(im,jm),psi(im,jm),phi1n(im,jm),
     &                  phi2n(im,jm),phi1t(im,jm),
     &                  phi2t(im,jm),phi2m(im,jm)
      common /gaugeblk/ alpha(im,jm),beta1u(im,jm),
     &                  beta2u(im,jm),alpha1t(im,jm),
     &                  beta1u1t(im,jm),beta2u1n(im,jm),
     &                  beta1u2t(im,jm),beta2u2m(im,jm)
      common /constblk/ sn(jm),cs(jm),ct(jm),sqrthfm,sqrt2
      common /weylblck/ psi0d(im,jm),psi1d(im,jm),psi2d(im,jm),
     &                  psi3d(im,jm),psi4d(im,jm)
      common /admmablk/ adm,igdet

      do 11 j=1,nthtamax
      do 11 i=1,netamax-1
      f1=psim4(i,j)
      fh=1./(psi(i,j)*psi(i,j))
      sqrta=sqrt(a(i,j))
      sqrtb=sqrt(b(i,j))
c
c psi0
      r1 =0.25*f1*phi1t(i,j)*ct(j)/b(i,j)
      r2 =0.25*beta2u1n(i,j)*fh*ct(j)/(sqrta*alpha(i,j))
      r3 =0.25*sqrta*beta1u1t(i,j)*fh*ct(j)/(alpha(i,j)*b(i,j))
      r4 =0.25*a1t(i,j)*f1*ct(j)/(a(i,j)*b(i,j))
      r5 =-0.25*f1*phi2t(i,j)/b(i,j)
      r6 =0.125*f1*phi1t(i,j)*phi1t(i,j)/b(i,j)
      r7 =-0.125*beta2u1n(i,j)*fh*phi1t(i,j)/(sqrta*alpha(i,j))
      r8 =-0.125*sqrta*beta1u1t(i,j)*fh*phi1t(i,j)/(alpha(i,j)*b(i,j))
      r9 =0.125*d1t(i,j)*f1*phi1t(i,j)/(b(i,j)*d(i,j))
      r10=0.125*b1t(i,j)*f1*phi1t(i,j)/(b(i,j)*b(i,j))
      r11=-0.25*hd(i,j)*fh*phi1n(i,j)/(sqrta*d(i,j))
      r12=0.25*hb(i,j)*fh*phi1n(i,j)/(sqrta*b(i,j))
      r13=0.125*d1n(i,j)*f1*phi1n(i,j)/(a(i,j)*d(i,j))
      r14=-0.125*b1n(i,j)*f1*phi1n(i,j)/(a(i,j)*b(i,j))
      r15=-0.5*hd1n(i,j)*fh/(sqrta*d(i,j))
      r16=0.25*d1n(i,j)*hd(i,j)*fh/(sqrta*d(i,j)*d(i,j))
      r17=0.5*hb1n(i,j)*fh/(sqrta*b(i,j))
      r18=-0.25*b1n(i,j)*hb(i,j)*fh/(sqrta*b(i,j)*b(i,j))
      r19=0.25*d1n(i,j)*ha(i,j)*fh/(sqrta*a(i,j)*d(i,j))
      r20=-0.25*b1n(i,j)*ha(i,j)*fh/(sqrta*a(i,j)*b(i,j))
      r21=0.125*beta2u1n(i,j)*d1t(i,j)*fh/(sqrta*alpha(i,j)*d(i,j))
      r22=0.125*sqrta*beta1u1t(i,j)*d1t(i,j)*fh
     &   /(alpha(i,j)*b(i,j)*d(i,j))
      r23=-0.25*beta2u2m(i,j)*fh/(sqrta*alpha(i,j))
      r24=-0.125*b1t(i,j)*beta2u1n(i,j)*fh/(sqrta*alpha(i,j)*b(i,j))
      r25=0.25*alpha1t(i,j)*beta2u1n(i,j)*fh/(sqrta*alpha(i,j)**2.)
      r26=-0.125*a1t(i,j)*beta2u1n(i,j)*fh/(sqrta*a(i,j)*alpha(i,j))
      r27=-0.25*sqrta*beta1u2t(i,j)*fh/(alpha(i,j)*b(i,j))
      r28=0.125*sqrta*b1t(i,j)*fh*beta1u1t(i,j)
     &   /(alpha(i,j)*b(i,j)*b(i,j))
      r29=0.25*sqrta*alpha1t(i,j)*beta1u1t(i,j)*fh
     &   /(alpha(i,j)*alpha(i,j)*b(i,j))
      r30=-0.375*a1t(i,j)*beta1u1t(i,j)*fh/(sqrta*alpha(i,j)*b(i,j))
      r31=0.125*a1t(i,j)*d1t(i,j)*f1/(a(i,j)*b(i,j)*d(i,j))
      r32=0.25*d2n(i,j)*f1/(a(i,j)*d(i,j))
      r33=-0.125*d1n(i,j)*d1n(i,j)*f1/(a(i,j)*d(i,j)*d(i,j))
      r34=-0.125*a1n(i,j)*d1n(i,j)*f1/(a(i,j)*a(i,j)*d(i,j))
      r35=0.125*a1t(i,j)*b1t(i,j)*f1/(a(i,j)*b(i,j)*b(i,j))
      r36=-0.25*b2n(i,j)*f1/(a(i,j)*b(i,j))
      r37=0.125*b1n(i,j)*b1n(i,j)*f1/(a(i,j)*b(i,j)*b(i,j))
      r38=0.125*a1n(i,j)*b1n(i,j)*f1/(a(i,j)*a(i,j)*b(i,j))
      r39=-0.25*a2t(i,j)*f1/(a(i,j)*b(i,j))
      r40=0.125*a1t(i,j)*a1t(i,j)*f1/(a(i,j)*a(i,j)*b(i,j))
      r41=-0.5*ha(i,j)*hd(i,j)/(a(i,j)*d(i,j))
      r42=-0.5*hc(i,j)*hc(i,j)/(a(i,j)*b(i,j))
      r43=0.5*beta1u1t(i,j)*hc(i,j)/(alpha(i,j)*b(i,j))
      r44=0.5*ha(i,j)*hb(i,j)/(a(i,j)*b(i,j))
      r45=-0.25*beta1u1t(i,j)*beta2u1n(i,j)/(alpha(i,j)*alpha(i,j))
      r46=-0.25*a(i,j)*(beta1u1t(i,j)**2.)
     &   /(alpha(i,j)*alpha(i,j)*b(i,j))
      psi0d(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &          +r10+r11+r12+r13+r14+r15+r16+r17+r18+r19
     &          +r20+r21+r22+r23+r24+r25+r26+r27+r28+r29
     &          +r30+r31+r32+r33+r34+r35+r36+r37+r38+r39
     &          +r40+r41+r42+r43+r44+r45+r46
      psi0d(i,j)=psi0d(i,j)*adm**2
c
c psi1
      r1 =-0.5*fh*hd(i,j)*ct(j)/(sqrtb*d(i,j))
      r2 =0.5*fh*hb(i,j)*ct(j)/(b(i,j)*sqrtb)
      r3 =0.25*d1n(i,j)*f1*ct(j)/(sqrta*sqrtb*d(i,j))
      r4 =-0.25*b1n(i,j)*f1*ct(j)/(sqrta*sqrtb*b(i,j))
      r5 =-0.125*f1*phi1n(i,j)*phi1t(i,j)/(sqrta*sqrtb)
      r6 =-0.25*fh*hd(i,j)*phi1t(i,j)/(sqrtb*d(i,j))
      r7 =0.25*fh*hb(i,j)*phi1t(i,j)/(sqrtb*b(i,j))
      r8 =-0.125*f1*b1n(i,j)*phi1t(i,j)/(sqrta*sqrtb*b(i,j))
      r9 =0.25*f1*phi2m(i,j)/(sqrta*sqrtb)
      r10=0.125*sqrtb*beta2u1n(i,j)*fh*phi1n(i,j)/(a(i,j)*alpha(i,j))
      r11=0.125*beta1u1t(i,j)*fh*phi1n(i,j)/(alpha(i,j)*sqrtb)
      r12=-0.125*a1t(i,j)*f1*phi1n(i,j)/(sqrta*a(i,j)*sqrtb)
      r13=-0.5*hd1t(i,j)*fh/(sqrtb*d(i,j))
      r14=0.25*d1t(i,j)*hd(i,j)*fh/(sqrtb*d(i,j)*d(i,j))
      r15=0.25*d1t(i,j)*hb(i,j)*fh/(sqrtb*b(i,j)*d(i,j))
      r16=0.125*sqrtb*beta2u1n(i,j)*d1n(i,j)*fh
     &   /(a(i,j)*alpha(i,j)*d(i,j))
      r17=0.125*beta1u1t(i,j)*d1n(i,j)*fh/(alpha(i,j)*sqrtb*d(i,j))
      r18=-0.125*d1n(i,j)*d1t(i,j)*f1/(sqrta*sqrtb*d(i,j)*d(i,j))
      r19=-0.125*b1n(i,j)*d1t(i,j)*f1/(sqrta*sqrtb*b(i,j)*d(i,j))
      r20=0.25*d2m(i,j)*f1/(sqrta*sqrtb*d(i,j))
      r21=-0.125*a1t(i,j)*d1n(i,j)*f1/(sqrta*a(i,j)*sqrtb*d(i,j))
      r22=-0.25*sqrtb*beta2u1n(i,j)*hd(i,j)/(sqrta*alpha(i,j)*d(i,j))
      r23=-0.25*sqrta*beta1u1t(i,j)*hd(i,j)/(alpha(i,j)*sqrtb*d(i,j))
      psi1d(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &          +r10+r11+r12+r13+r14+r15+r16+r17+r18+r19
     &          +r20+r21+r22+r23
      psi1d(i,j)=psi1d(i,j)*adm**2
c
c psi2
      r1 =0.25*f1*phi1t(i,j)*ct(j)/(b(i,j))
      r2 =0.5*d1t(i,j)*f1*ct(j)/(b(i,j)*d(i,j))
      r3 =-0.25*b1t(i,j)*f1*ct(j)/(b(i,j)*b(i,j))
      r4 =0.25*f1*phi2t(i,j)/(b(i,j))
      r5 =0.125*d1t(i,j)*f1*phi1t(i,j)/(b(i,j)*d(i,j))
      r6 =-0.125*b1t(i,j)*f1*phi1t(i,j)/(b(i,j)*b(i,j))
      r7 =0.125*f1*phi1n(i,j)*phi1n(i,j)/a(i,j)
      r8 =0.125*d1n(i,j)*f1*phi1n(i,j)/(a(i,j)*d(i,j))
      r9 =0.125*b1n(i,j)*f1*phi1n(i,j)/(a(i,j)*b(i,j))
      r10=0.25*d2t(i,j)*f1/(b(i,j)*d(i,j))
      r11=-0.125*d1t(i,j)*d1t(i,j)*f1/(b(i,j)*d(i,j)*d(i,j))
      r12=-0.125*b1t(i,j)*d1t(i,j)*f1/(b(i,j)*b(i,j)*d(i,j))
      r13=0.125*b1n(i,j)*d1n(i,j)*f1/(a(i,j)*b(i,j)*d(i,j))
      r14=-0.5*f1/b(i,j)
      r15=-0.5*hb(i,j)*hd(i,j)/(b(i,j)*d(i,j))
      psi2d(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &          +r10+r11+r12+r13+r14+r15
      psi2d(i,j)=psi2d(i,j)*adm**2
c
c psi3
      r1 =-0.5*fh*hd(i,j)*ct(j)/(sqrtb*d(i,j))
      r2 =0.5*fh*hb(i,j)*ct(j)/(b(i,j)*sqrtb)
      r3 =-0.25*d1n(i,j)*f1*ct(j)/(sqrta*sqrtb*d(i,j))
      r4 =0.25*b1n(i,j)*f1*ct(j)/(sqrta*sqrtb*b(i,j))
      r5 =0.125*f1*phi1n(i,j)*phi1t(i,j)/(sqrta*sqrtb)
      r6 =-0.25*fh*hd(i,j)*phi1t(i,j)/(sqrtb*d(i,j))
      r7 =0.25*fh*hb(i,j)*phi1t(i,j)/(sqrtb*b(i,j))
      r8 =0.125*f1*b1n(i,j)*phi1t(i,j)/(sqrta*sqrtb*b(i,j))
      r9 =-0.25*f1*phi2m(i,j)/(sqrta*sqrtb)
      r10=0.125*sqrtb*beta2u1n(i,j)*fh*phi1n(i,j)/(a(i,j)*alpha(i,j))
      r11=0.125*beta1u1t(i,j)*fh*phi1n(i,j)/(alpha(i,j)*sqrtb)
      r12=0.125*a1t(i,j)*f1*phi1n(i,j)/(sqrta*a(i,j)*sqrtb)
      r13=-0.5*hd1t(i,j)*fh/(sqrtb*d(i,j))
      r14=0.25*d1t(i,j)*hd(i,j)*fh/(sqrtb*d(i,j)*d(i,j))
      r15=0.25*d1t(i,j)*hb(i,j)*fh/(sqrtb*b(i,j)*d(i,j))
      r16=0.125*sqrtb*beta2u1n(i,j)*d1n(i,j)*fh
     &   /(a(i,j)*alpha(i,j)*d(i,j))
      r17=0.125*beta1u1t(i,j)*d1n(i,j)*fh/(alpha(i,j)*sqrtb*d(i,j))
      r18=0.125*d1n(i,j)*d1t(i,j)*f1/(sqrta*sqrtb*d(i,j)*d(i,j))
      r19=0.125*b1n(i,j)*d1t(i,j)*f1/(sqrta*sqrtb*b(i,j)*d(i,j))
      r20=-0.25*d2m(i,j)*f1/(sqrta*sqrtb*d(i,j))
      r21=0.125*a1t(i,j)*d1n(i,j)*f1/(sqrta*a(i,j)*sqrtb*d(i,j))
      r22=0.25*sqrtb*beta2u1n(i,j)*hd(i,j)/(sqrta*alpha(i,j)*d(i,j))
      r23=0.25*sqrta*beta1u1t(i,j)*hd(i,j)/(alpha(i,j)*sqrtb*d(i,j))
      psi3d(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &          +r10+r11+r12+r13+r14+r15+r16+r17+r18+r19
     &          +r20+r21+r22+r23
      psi3d(i,j)=psi3d(i,j)*adm**2
c
c psi4
      r1 =0.25*f1*phi1t(i,j)*ct(j)/b(i,j)
      r2 =-0.25*beta2u1n(i,j)*fh*ct(j)/(sqrta*alpha(i,j))
      r3 =-0.25*sqrta*beta1u1t(i,j)*fh*ct(j)/(alpha(i,j)*b(i,j))
      r4 =0.25*a1t(i,j)*f1*ct(j)/(a(i,j)*b(i,j))
      r5 =-0.25*f1*phi2t(i,j)/b(i,j)
      r6 =0.125*f1*phi1t(i,j)*phi1t(i,j)/b(i,j)
      r7 =0.125*beta2u1n(i,j)*fh*phi1t(i,j)/(sqrta*alpha(i,j))
      r8 =0.125*sqrta*beta1u1t(i,j)*fh*phi1t(i,j)/(alpha(i,j)*b(i,j))
      r9 =0.125*d1t(i,j)*f1*phi1t(i,j)/(b(i,j)*d(i,j))
      r10=0.125*b1t(i,j)*f1*phi1t(i,j)/(b(i,j)*b(i,j))
      r11=0.25*hd(i,j)*fh*phi1n(i,j)/(sqrta*d(i,j))
      r12=-0.25*hb(i,j)*fh*phi1n(i,j)/(sqrta*b(i,j))
      r13=0.125*d1n(i,j)*f1*phi1n(i,j)/(a(i,j)*d(i,j))
      r14=-0.125*b1n(i,j)*f1*phi1n(i,j)/(a(i,j)*b(i,j))
      r15=0.5*hd1n(i,j)*fh/(sqrta*d(i,j))
      r16=-0.25*d1n(i,j)*hd(i,j)*fh/(sqrta*d(i,j)*d(i,j))
      r17=-0.5*hb1n(i,j)*fh/(sqrta*b(i,j))
      r18=0.25*b1n(i,j)*hb(i,j)*fh/(sqrta*b(i,j)*b(i,j))
      r19=-0.25*d1n(i,j)*ha(i,j)*fh/(sqrta*a(i,j)*d(i,j))
      r20=0.25*b1n(i,j)*ha(i,j)*fh/(sqrta*a(i,j)*b(i,j))
      r21=-0.125*beta2u1n(i,j)*d1t(i,j)*fh/(sqrta*alpha(i,j)*d(i,j))
      r22=-0.125*sqrta*beta1u1t(i,j)*d1t(i,j)*fh
     &   /(alpha(i,j)*b(i,j)*d(i,j))
      r23=0.25*beta2u2m(i,j)*fh/(sqrta*alpha(i,j))
      r24=0.125*b1t(i,j)*beta2u1n(i,j)*fh/(sqrta*alpha(i,j)*b(i,j))
      r25=-0.25*alpha1t(i,j)*beta2u1n(i,j)*fh/(sqrta*alpha(i,j)**2.)
      r26=0.125*a1t(i,j)*beta2u1n(i,j)*fh/(sqrta*a(i,j)*alpha(i,j))
      r27=0.25*sqrta*beta1u2t(i,j)*fh/(alpha(i,j)*b(i,j))
      r28=-0.125*sqrta*b1t(i,j)*fh*beta1u1t(i,j)
     &   /(alpha(i,j)*b(i,j)*b(i,j))
      r29=-0.25*sqrta*alpha1t(i,j)*beta1u1t(i,j)*fh
     &   /(alpha(i,j)*alpha(i,j)*b(i,j))
      r30=0.375*a1t(i,j)*beta1u1t(i,j)*fh/(sqrta*alpha(i,j)*b(i,j))
      r31=0.125*a1t(i,j)*d1t(i,j)*f1/(a(i,j)*b(i,j)*d(i,j))
      r32=0.25*d2n(i,j)*f1/(a(i,j)*d(i,j))
      r33=-0.125*d1n(i,j)*d1n(i,j)*f1/(a(i,j)*d(i,j)*d(i,j))
      r34=-0.125*a1n(i,j)*d1n(i,j)*f1/(a(i,j)*a(i,j)*d(i,j))
      r35=0.125*a1t(i,j)*b1t(i,j)*f1/(a(i,j)*b(i,j)*b(i,j))
      r36=-0.25*b2n(i,j)*f1/(a(i,j)*b(i,j))
      r37=0.125*b1n(i,j)*b1n(i,j)*f1/(a(i,j)*b(i,j)*b(i,j))
      r38=0.125*a1n(i,j)*b1n(i,j)*f1/(a(i,j)*a(i,j)*b(i,j))
      r39=-0.25*a2t(i,j)*f1/(a(i,j)*b(i,j))
      r40=0.125*a1t(i,j)*a1t(i,j)*f1/(a(i,j)*a(i,j)*b(i,j))
      r41=-0.5*ha(i,j)*hd(i,j)/(a(i,j)*d(i,j))
      r42=-0.5*hc(i,j)*hc(i,j)/(a(i,j)*b(i,j))
      r43=0.5*beta1u1t(i,j)*hc(i,j)/(alpha(i,j)*b(i,j))
      r44=0.5*ha(i,j)*hb(i,j)/(a(i,j)*b(i,j))
      r45=-0.25*beta1u1t(i,j)*beta2u1n(i,j)/(alpha(i,j)*alpha(i,j))
      r46=-0.25*a(i,j)*(beta1u1t(i,j)**2.)
     &   /(alpha(i,j)*alpha(i,j)*b(i,j))
      psi4d(i,j)=    r1 +r2 +r3 +r4 +r5 +r6 +r7 +r8 +r9
     &           +r10+r11+r12+r13+r14+r15+r16+r17+r18+r19
     &           +r20+r21+r22+r23+r24+r25+r26+r27+r28+r29
     &           +r30+r31+r32+r33+r34+r35+r36+r37+r38+r39
     &           +r40+r41+r42+r43+r44+r45+r46
      psi4d(i,j)=psi4d(i,j)*adm**2
   11 continue
      return
      end
 
 
 
 
 
