program main

C     Numerical solution of the 3-D wave equation
C     using stereographic coordinate patches for the 
C     calculation of the angular momentum operator. 
C_______________________________________________________________________
C
C     Begin: user - configurable parameters
C
C     Grid size

      integer          nd,     nx
      parameter       (nd = 8, nx = 32)

C     Safety margin factor, du = margin * dx * dd ** 2
C     The Courant condition might be violated if margin .gt. 1

      double precision margin
      parameter       (margin = 1.)

C     Starting and final evolution time

      double precision begintime,      endtime
      parameter       (begintime = 0., endtime = 1.)

C     Number of dumps to make - a value .lt. 0 indicates no output,
C     not even from the "norm" routine. Useful only to perform timing 
C     comparisons.

      integer          dumps
      parameter       (dumps = 10)

C     Particular analytic solution to check against. The solution
C     is a pure Y_{ll,mm} spherical harmonic.

      integer          ll,      mm
      parameter       (ll = 2, mm = 2)

C     End: user - configurable parameters
C_______________________________________________________________________

C     Time, grid spacings, iteration counters

      double precision time
      double precision dd, dx, du
      integer          i, it, nt, scoop

C     Grid arrays, field values, mask array

      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      integer          mask(-(nd+1):(nd+1),-(nd+1):(nd+1))

C     Create grid functions, calculate timestep

      call setup(mask, x, q, p, nd, nx, dd, dx, du, ll, mm, 
     &           margin, time, begintime, endtime, nt, dumps, scoop)

C     Load initial data on the (new) field array
C     Then interpolate and copy to (old) array

      call initial(gn, x, q, p, nd, nx, time, ll, mm)
      do i = 0, nx
        call interpolate(gn, q, p, mask, nd, nx, dd, i)
      end do
      call copy(go, gn, nd, nx)

C     Do some output of the initial data if requested

      if (dumps .gt. 0) then
        call monitor(gn, x, q, p, nd, nx, it, ll, mm, time)
      end if

C     Now do the numerical calculation

      do it = 1, nt

C       Advance time

        time = time + du

C       Do the radial march outwards
 
        i = 1
        call gstart(go, gn, x, q, p, nd, nx, dd, dx, du)
        call interpolate(gn, q, p, mask, nd, nx, dd, i)

        do i = 2, nx - 1
          call gmiddle(go, gn, x, q, p, nd, nx, dd, dx, du, i) 
          call interpolate(gn, q, p, mask, nd, nx, dd, i)
        end do

        i = nx
        call gend(go, gn, x, q, p, nd, nx, dd, dx, du)
        call interpolate(gn, q, p, mask, nd, nx, dd, i)
  
C       Shift the field values and do some output

        call norm(gn, x, q, p, nd, nx, dd, dx, time, dumps, ll, mm)
        call copy(go, gn, nd, nx)
        if (dumps .gt. 0) then
          if ((mod(it,scoop) .eq. 0) .or. (it .eq. nt)) then
            call monitor(gn, x, q, p, nd, nx, it, ll, mm, time)
          end if
        end if

      end do 

      end
      subroutine setup(mask, x, q, p, nd, nx, dd, dx, du, ll, mm, 
     &                 margin, time, begintime, endtime, nt, 
     &                 dumps, scoop)

      integer mask(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer nd, nx
      double precision dd, dx, du
      integer ll, mm
      double precision margin, time, begintime, endtime
      integer nt, dumps, scoop

      integer i

C     Fill-in grid coordinate functions
C     ND is the real grid dimension, the grid is enlarged by one,
C     so that it holds the interpolated values as well.

      dd = 1. / dble(nd)
      dx = 1. / dble(nx)

      do i = 0, nx 
         x(i)  = i * dx
      end do 

      do i = -(nd+1), (nd+1)
         q(i)  = i * dd
         p(i)  = i * dd
      end do

C     Figure out the number of iterations (nt)
C     Calculate the timestep (du), and the sampling density (scoop).

      time = begintime
      du = margin * dx * dd ** 2

      nt = int((endtime - du) / du) + 1
      if (nt .lt. dumps) dumps = nt
      if (dumps .gt. 0) then
        scoop = int(nt / dumps)
      else
        scoop = 1
      end if

      if (dumps .gt. 0) then
        open (11,file='dat/run.log' ,status='new')
        write (11,*) 'Run Parameters'
        write (11,*) '-------------------------------------'
        write (11,*) 'Grid size            : ', nx, ' x ', nd, ' x ', nd
        write (11,*) 'Initial time         : ', time
        write (11,*) 'Final time           : ', endtime
        write (11,*) 'Number of iterations : ', nt
        write (11,*) 'Number of data dumps : ', dumps
        write (11,*) ' '
        call flush(11) 
      end if

C     Set mask array for interpolation between stereographic patches

      call setmask(mask, q, p, nd)

      return
      end
      subroutine initial(gn, x, q, p, nd, nx, time, ll, mm)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      double precision time
      integer          nd, nx, ll, mm

      integer          i, k, l
      double precision u, fg

C     Fill in initial data for all x's, both patches, all angles

      u = time

      do i = 0, nx 
        do l = -nd, nd
          do k = -nd, nd
            gn(k,l,i,0) = fg(ll,mm,u,q(k),p(l),x(i),0)
            gn(k,l,i,1) = fg(ll,mm,u,q(k),p(l),x(i),1)
         end do
       end do
      end do

      return
      end

      double precision function fg(ll, mm, u, q, p, x, D)
      double precision u, q, p, x, y
      double precision lambda
      double precision cos1, ylm, cos2, phi
      integer ll,mm, D, l, m

C     Choose overall factor

      lambda = 1.
      l = ll
      m = mm

C     Calculate the cosine for this point and this patch

      y = cos1(q,p,D)
      phi = acos(cos2(q,p,D))
      if(abs(cos2(q,p,D)).gt.1) then
      print*, q, p, cos2(q,p,D)
      end if

      fg = lambda * 2**(l+1) *  x ** (l+1) * ylm(l, m, y, phi)
     &   / ((u + 1.) ** (l + 1) * (u * (1 - x) + x + 1.) ** (l+1))

      return 
      end 

      double precision function legendre(l,y)

C     Legendre Polynomials

      integer l, lc
      double precision y
      double precision pnext, tmp0, tmp1

C     l = 0 and 1 polynomials 

      tmp0 = 1.
      tmp1 = y

      if(l.eq.0) then
         legendre = tmp0
      else if (l.eq.1) then 
         legendre = tmp1
      else 

C     Higher order recursive formula

         do lc = 2, l
            pnext = (y * (2*lc-1) * tmp1 - (lc-1) * tmp0) / lc
            tmp0 = tmp1
            tmp1 = pnext
         end do

	 legendre = pnext

      endif
      return
      end

      double precision function ylm(l, m, x, phi)

C     Spherical Harmonics
C     Re(Ylm) (The real part)

      integer l, m
      double precision x, phi
      double precision pi, norm
      double precision plgndr, factrl

      pi = 3.14159265358979323846

      norm = sqrt((2. * l + 1.) / (4. * pi) 
     &       * factrl(l - m) / factrl(l + m))

      ylm = norm * plgndr(l, m, x) * cos(m * phi)

      return
      end

      double precision function plgndr(l,m,x)

C     Associated Legendre Polynomials

      integer l,m
      double precision x
      integer i,ll
      double precision fact,pll,pmm,pmmp1,somx2

      if(m.lt.0.or.m.gt.l.or.abs(x).gt.1.) then
         stop 'bad arguments to plgndr'
      end if

      pmm=1.
      if(m.gt.0) then
	somx2=sqrt((1.-x)*(1.+x))
	fact=1.
	do i=1,m
		pmm=-pmm*fact*somx2
		fact=fact+2.
	enddo
      endif

      if(l.eq.m) then
         plgndr=pmm
      else
         pmmp1=x*(2*m+1)*pmm
         if(l.eq.m+1) then
            plgndr=pmmp1
	 else
            do ll=m+2,l
		pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
		pmm=pmmp1
		pmmp1=pll
            enddo
	plgndr=pll
	endif
      endif
      return
      end

      double precision function factrl(n)

C     Factorial

      integer n
      double precision gammln
      integer j,ntop
      double precision a(33)
      save ntop,a
      data ntop,a(1)/0,1./

      if (n.lt.0) then
          stop 'negative factorial'
      else if (n.le.ntop) then
          factrl=a(n+1)
      else if (n.le.32) then
          do j=ntop+1,n
            a(j+1)=j*a(j)
          end do
          ntop=n
          factrl=a(n+1)
      else
      factrl = dexp(gammln(dble(n)+1.0D0))
      endif
      return
      end

      double precision function gammln(xx)
      double precision xx
      double precision half,one,fpf
      parameter (half=0.5d0, one=1.0d0, fpf=5.5d0)
      integer j
      double precision ser,stp,tmp,x,cof(6)
      save cof,stp
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     & -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/

      x = xx - one
      tmp = x + fpf
      tmp = (x+half) * log(tmp) - tmp
      ser = one

      do j=1,6
         x=x+one
         ser=ser+cof(j)/x
      end do
      gammln = tmp+log(stp*ser)

      return
      end

      double precision function p1(u,q,p,x,D)
      double precision u, q, p, x
      double precision lambda
      double precision cos2, sin1
      integer D
      
      lambda = 1.

      p1 = lambda * 4 * sin1(q,p,D) * cos2(q,p,D) * x**2 /
     &    ( (-x - u + u * x - 1)**2 * (u + 1)**2)


      return 
      end 

      double precision function sin1(q,p,D)
      double precision q, p  
      integer D
      sin1 = 2. * sqrt(q*q + p*p) / (1. + q*q + p*p) 
      if (sin1.gt.1) print*, 'Error in sin1'
      end

      double precision function sin2(q,p,D)
      double precision q, p  
      integer D
      if(q.eq.0.and.p.eq.0) then
        sin2 = 0.
      else
          if(D.eq.0) then 
               sin2 =   p / sqrt(q*q + p*p)
          else if(D.eq.1) then
               sin2 = - p / sqrt(q*q + p*p)
          end if 
      end if
      if (sin2.gt.1) print*, 'Error in sin2'
      end

      double precision function cos1(q,p,D)
      double precision q, p  
      integer D
      if(D.eq.0) then 
         cos1 = (q*q + p*p -1.) / (1. + q*q + p*p)
      else if(D.eq.1) then
         cos1 = (1. - q*q - p*p) / (1. + q*q + p*p) 
      end if
      if (cos1.gt.1) print*, 'Error in cos1'
      end

      double precision function cos2(q,p,D)
      double precision q, p  
      integer D
      if(q.eq.0.and.p.eq.0) then
        cos2 = 0.
      else if(p.eq.0.and.q.gt.0) then
        cos2 = 1.
      else if(p.eq.0.and.q.lt.0) then
        cos2 = - 1.
      else 
        cos2 = q / sqrt(q*q + p*p)
      end if
      if (cos2.gt.1) cos2 = 1.
      end
      subroutine copy(go, gn, nd, nx)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      integer          nd, nx

      integer          i, k, l

C     Copy the "new" array into the "old" one. We could use a
C     "canned" routine from a library to do this chore.

      do i = 0, nx
        do l = -(nd+1), (nd+1)
          do k = -(nd+1), (nd+1)
             go(k,l,i,0) = gn(k,l,i,0)
             go(k,l,i,1) = gn(k,l,i,1)
          end do
        end do
      end do

      return
      end
      subroutine gstart(go, gn, x, q, p, nd, nx, dd, dx, du)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, du

      integer          D, i, k, l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs
      double precision rp, rq, rr, rs, rc
      double precision gp, gq, gr, gs
      double precision l2, l2im1, l2i, factor, area
      double precision gmean

C     Calculate the characteristic deviation 

      i = 1
      dold = 0.5  * du * (1. - x(i-1)) ** 2 
      dnew = 0.25 * du * (1. - x(i))   ** 2 
      d1 = dold / dx 
      d2 = dnew / dx 

      xp = x(i-1) 
      xq = x(i)   - dnew
      xr = x(i-1) + dold
      xs = x(i)   + dnew

      rp = xp / (1. - xp)
      rq = xq / (1. - xq)
      rr = xr / (1. - xr)
      rs = xs / (1. - xs)

      rc = 0.5 * (rp + rs)

C     area = \int_\Sigma du dr

      area = 0.5 * du * (rq - rp + rs - rr)
 
C     Calculate gn for both patches independently

      do D = 0, 1
        gmean = 0.
        do l = -nd, nd
          do k = -nd, nd

            factor = - 0.25 * (1 + q(k) ** 2 + p(l) ** 2) ** 2
            l2im1 = 0.
            l2i = factor * (1. / dd ** 2) * 
     &            (go(k+1,l,i,D)
     &           + go(k-1,l,i,D) + go(k,l+1,i,D)
     &           + go(k,l-1,i,D) - 4. * go(k,l,i,D)) 
            l2 = 0.5 * (l2im1 + l2i)

C           Interpolate the g = r * f values

            gp = 0.

            gr = (1. - d1) * go(k,l,i-1,D) + d1 * go(k,l,i,D)

            gs = (1. - d2) * go(k,l,i,D)   + d2 * go(k,l,i+1,D)

            gq =  gp + gs - gr - 0.5 * area * l2 / rc ** 2

C           Extrapolate to obtain the new gamma grid value

            gn(k,l,i,D) = (gq - gn(k,l,i-1,D) * d2) / (1. - d2)
            gmean = gmean + gn(k,l,i,D)
          end do
        end do

        gmean = gmean / (2.* nd + 1.) ** 2

C       Average out the angular dependence

        do l = -nd, nd
          do k = -nd, nd
            gn(k,l,i,D) = gmean
          end do
        end do
      end do

      return
      end 
      subroutine gmiddle(go, gn, x, q, p, nd, nx, dd, dx, du, i)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, du
      integer          i

      integer          D, k, l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs, xc
      double precision rp, rq, rr, rs
      double precision gp, gq, gr, gs
      double precision l2, l2im1, l2i, factor, integral

C     Calculate the characteristic deviation 

        dold = 0.25 * du * (1. - x(i-1)) ** 2
        dnew = 0.25 * du * (1. - x(i))   ** 2 

        d1 = dold / dx 
        d2 = dnew / dx 

        xp = x(i-1) - dold
        xq = x(i)   - dnew
        xr = x(i-1) + dold
        xs = x(i)   + dnew

        rp = xp / (1. - xp)
        rr = xr / (1. - xr)
        rq = xq / (1. - xq)
        rs = xs / (1. - xs)
      
        xc = 0.5 * (xp + xs)

C       integral = \int_\Sigma du dr 1/r^2

        integral = 2. * log((rq * rr) / (rp * rs))
 
C       Calculate gn for both patches independently

        do D = 0, 1
          do l = -nd, nd
            do k = -nd, nd

C             Calculate the Laplacian
       
              factor = - 0.25 * (1 + q(k) ** 2 + p(l) ** 2) ** 2

              l2im1 = factor * (1. / dd ** 2) * 
     &                (gn(k+1,l,i-1,D) + gn(k-1,l,i-1,D)
     &               + gn(k,l+1,i-1,D) + gn(k,l-1,i-1,D)
     &               - 4. * gn(k,l,i-1,D)) 

              l2i = factor * (1. / dd ** 2) * 
     &              (go(k+1,l,i,D) + go(k-1,l,i,D)
     &             + go(k,l+1,i,D) + go(k,l-1,i,D)
     &             - 4. * go(k,l,i,D)) 

              l2 = 0.5 * (l2im1 / x(i-1) ** 2 +  l2i / x(i) ** 2)

C             Interpolate the g = r * f values

              gp = (1. - d1) * gn(k,l,i-1,D) + d1 * gn(k,l,i-2,D)

              gr = (1. - d1) * go(k,l,i-1,D) + d1 * go(k,l,i,D)

              gs = (1. - d2) * go(k,l,i,D)   + d2 * go(k,l,i+1,D)

C             Calculate the gq box value

              gq = gp + gs - gr - 0.5 * integral * l2 * xc ** 2

C             Extrapolate to obtain the new gamma grid value

              gn(k,l,i,D) = (gq - gn(k,l,i-1,D) * d2) / (1. - d2)

            end do
          end do
        end do

      return
      end 
      subroutine gend(go, gn, x, q, p, nd, nx, dd, dx, du)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, du

      integer          D, i, k, l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs
      double precision rp,  rr
      double precision gp, gq, gr, gs
      double precision l2, l2im1, l2i, factor, integral

C     Calculate the characteristic deviation 

      i = nx
      dold = 0.25 * du * (1. - x(i-1)) ** 2
      dnew = 0.
      d1 = dold / dx 
      d2 = dnew / dx 

      xp = x(i-1) - dold
      xq = 1.
      xr = x(i-1) + dold
      xs = 1.

      rp = xp / (1. - xp)
      rr = xr / (1. - xr)

C     integral = \int_\Sigma du dr 1/r^2

      integral = 2. * log(rr / rp)
 
C     Calculate gn for both patches independently

      do D = 0, 1
        do l = -nd, nd
          do k = -nd, nd

C           Calculate the Laplacian
       
            factor = - 0.25 * (1 + q(k) ** 2 + p(l) ** 2) ** 2

            l2im1 = factor * (1. / dd ** 2)
     &            * (gn(k+1,l,i-1,D) + gn(k-1,l,i-1,D) 
     &             + gn(k,l+1,i-1,D) + gn(k,l-1,i-1,D) 
     &             - 4. * gn(k,l,i-1,D)) 

            l2i = factor * (1. / dd ** 2)  
     &          * (go(k+1,l,i,D) + go(k-1,l,i,D) 
     &           + go(k,l+1,i,D) + go(k,l-1,i,D)
     &           - 4. * go(k,l,i,D)) 

            l2 = 0.5 * (l2im1 + l2i)

C           Interpolate the g = r * f values

            gp = (1. - d1) * gn(k,l,i-1,D) + d1 * gn(k,l,i-2,D)

            gr = (1. - d1) * go(k,l,i-1,D) + d1 * go(k,l,i,D)

            gs = go(k,l,i,D)

C           Calculate the gp box value

            gq = gp + gs - gr - 0.5 * integral * l2

C           Extrapolate to obtain the new gamma grid value

            gn(k,l,i,D) = (gq - gn(k,l,i-1,D) * d2) / (1. - d2)

          end do
        end do
      end do

      return
      end 
      subroutine monitor(gn, x, q, p, nd, nx, it, ll, mm, time)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx, it, ll, mm
      double precision time

      integer          i, k, l
      double precision fg

      character*19 gefile, gnfile
      character*12 ctime

      open(60, file = 'dat/ge.gnu', status = 'unknown') 
      open(61, file = 'dat/gn.gnu', status = 'unknown') 

      write (ctime,'(D12.6E2)') time
      gefile = 'dat/ge_' // ctime
      gnfile = 'dat/gn_' // ctime
      open (30, file = gefile, status = 'unknown') 
      open (31, file = gnfile, status = 'unknown') 

      i = nx
      l = 0

      do i = 0, nx 
        do k = -nd, nd
          write (30,10) fg(ll,mm,time,q(k),p(l),x(i),0) - gn(k,l,i,0)
          write (31,10) gn(k,l,i,0) 
        end do
        write(30,20)
        write(31,20)
      end do

      close (30)
      close (31)

      write (60, '(a)') 'splot "' // gefile(5:19) // '"'
      write (61, '(a)') 'splot "' // gnfile(5:19) // '"'
      call flush (60)
      call flush (61)

   10 format(3E23.12E3)
   20 format('')

      return
      end
      subroutine interpolate(gn, q, p, mask, nd, nx, dd, i)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          mask(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          nd, nx, i
      double precision dd

      double precision z
      parameter       (z = 1.)
      integer          D, U, k, l
      integer          jg, ig, ig2, jg2
      integer          k1, k2, l1, l2
      double precision qf, pf, cp, cq, cp1, cp2, cq1, cq2
      double precision q1, q2, p1, p2, v1, v2, q1p, q2p, p1p, p2p
      double precision r1, r2, r3, r4, rx, qd1, qd2, pd1, pd2
      double precision q_s2n, p_s2n, cubic

C     Fourth order interpolation at the patch boundary
C     Upper is the grid we are looking at (rectangular)
C     Lower is the other patch, with circular arc coordinate lines

C     Do both grids 

      do D = 0,1

        if(D.eq.0) U = 1
        if(D.eq.1) U = 0

        do k= -(nd+1), (nd+1)
          do l= -(nd+1), (nd+1)

C           Look for the boundary

            if(mask(k,l).ne.0) then

C             Find the coordinate of the fictitious point in the lower 
C             grid.

              qf = q_s2n(q(k),p(l))
              pf = p_s2n(q(k),p(l))
 
C             Find the guard point index (see notes for formulae)
C             Also find the outer neighboor index.

              if(mask(k,l).eq.11.or.mask(k,l).eq.21) then
                ig = int(qf/dd)
                jg = int(pf/dd) + int((sign(z,pf) - 1)/2)
              else if(mask(k,l).eq.12.or.mask(k,l).eq.22) then
                ig = int(qf/dd) + int((sign(z,qf) - 1)/2)
                jg = int(pf/dd)
              else if(mask(k,l).eq.3) then
                ig = int(qf/dd) 
                jg = int(pf/dd)
              end if
              ig2 = ig - int(sign(z,q(ig)))
              jg2 = jg - int(sign(z,p(jg)))

C             Find the upper grid points used in the cubic interpolation

              k1 = k - 1 * int(sign(z,q(k)))
              k2 = k - 2 * int(sign(z,q(k)))
              l1 = l - 1 * int(sign(z,p(l)))
              l2 = l - 2 * int(sign(z,p(l)))

C             Give symbolic names to the constant coordinate lines
C             involved

              cp   = p(l)
              cq   = q(k)
              cq1 = q(ig)
              cq2 = q(ig2)
              cp1 = p(jg)
              cp2 = p(jg2)

C             Calculate the intersection points depending on which 
C             combination of crossing lines we use

              if(mask(k,l).eq.11) then 
                q1  = (1. + sqrt(1. - 4. * cp**2 * cq1**2))/(2.*cq1)
                p1p = -2.*cp*cq1**2/(1.+ sqrt(1. - 4. * cp**2 *cq1**2))
                q2  = (1. + sqrt(1. - 4. * cp**2 * cq2**2))/(2.*cq2)
                p2p = -2.*cp*cq2**2/(1.+ sqrt(1. - 4. * cp**2 *cq2**2))
              else if(mask(k,l).eq.21) then
                p1  = sign(z,p(l)) * sqrt(cq/cq1 - cq**2)
                p1p = - (cq1/cq) * p1
                p2  = sign(z,p(l)) * sqrt(cq/cq2 - cq**2)
                p2p = - (cq2/cq) * p2
              else if(mask(k,l).eq.12) then
                q1  = sign(z,q(k)) * sqrt(-cp/cp1 - cp**2)
                q1p = - (cp1/cp) * q1
                q2  = sign(z,q(k)) * sqrt(-cp/cp2 - cp**2)
                q2p = - (cp2/cp) * q2
              else if(mask(k,l).eq.22) then
                p1  = (-1. - sqrt(1. - 4. * cq**2 * cp1**2))/(2.*cp1)
                q1p = 2.*cq*cp1**2/(1. + sqrt(1. - 4. * cq**2 *cp1**2))
                p2  = (-1. - sqrt(1. - 4. * cq**2 * cp2**2))/(2.*cp2)
                q2p = 2.*cq*cp2**2/(1. + sqrt(1. - 4. * cq**2 *cp2**2))
              end if 

C             Do the lower interpolation depending on which direction 
C             is used

              if(mask(k,l).eq.11.or.mask(k,l).eq.21) then 

                v1 = cubic(p1p, p(jg-1), p(jg), p(jg+1), p(jg+2)
     &                    ,gn(ig,jg-1,i,U), gn(ig,jg,i,U)
     &                    ,gn(ig,jg+1,i,U), gn(ig,jg+2,i,U))

                v2 = cubic(p2p, p(jg-1), p(jg), p(jg+1), p(jg+2)
     &                    ,gn(ig2,jg-1,i,U), gn(ig2,jg,i,U)
     &                    ,gn(ig2,jg+1,i,U), gn(ig2,jg+2,i,U))

              else if(mask(k,l).eq.12.or.mask(k,l).eq.22) then

                v1 = cubic(q1p, q(ig-1), q(ig), q(ig+1), q(ig+2)
     &                    ,gn(ig-1,jg,i,U), gn(ig,jg,i,U)
     &                    ,gn(ig+1,jg,i,U), gn(ig+2,jg,i,U))

                v2 = cubic(q2p, q(ig-1), q(ig), q(ig+1), q(ig+2)
     &                    ,gn(ig-1,jg2,i,U), gn(ig,jg2,i,U)
     &                    ,gn(ig+1,jg2,i,U), gn(ig+2,jg2,i,U))

              else if(mask(k,l).eq.3) then

                v1 = gn(ig,jg,i,U)
                v2 = gn(ig2,jg2,i,U)

              end if

C             Do the upper interpolation depending on which direction 
C             is used

              if(mask(k,l).eq.11.or.mask(k,l).eq.12) then

                gn(k,l,i,D) = cubic(q(k), q(k2), q(k1), q1, q2
     &                             ,gn(k2,l,i,D), gn(k1,l,i,D), v1, v2)
             
              else if(mask(k,l).eq.21.or.mask(k,l).eq.22) then

                gn(k,l,i,D) = cubic(p(l), p(l2), p(l1), p1, p2
     &                   ,gn(k,l2,i,D), gn(k,l1,i,D), v1, v2)

              else if(mask(k,l).eq.3) then

                qd1 = q_s2n(q(ig),p(jg))
                pd1 = p_s2n(q(ig),p(jg))
                qd2 = q_s2n(q(ig2),p(jg2))
                pd2 = p_s2n(q(ig2),p(jg2))
                    
                r1 = sqrt(q(k1)**2+p(l1)**2)
                r2 = sqrt(q(k2)**2+p(l2)**2)
                rx = sqrt(q(k)**2+p(l)**2)
                r3 = sqrt(qd1**2+pd1**2)
                r4 = sqrt(qd2**2+pd2**2)

                gn(k,l,i,D) = cubic(rx,r1,r2,r3,r4
     &                             ,gn(k1,l1,i,D),gn(k2,l2,i,D),v1,v2)

              end if
            end if
          end do
        end do
      end do

      return
      end

      subroutine setmask(mask, q, p, nd)
      integer mask(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer nd

      integer i, j

C     Set mask = 0 for interior points
C     Set mask = 11 for x edge and interpolation along constant x'
C     Set mask = 12 for x edge and interpolation along constant y'
C     Set mask = 21 for y edge and interpolation along constant x'
C     Set mask = 22 for y edge and interpolation along constant y'
C     Set mask = 3 for vertices(?)

      do i = -(nd+1), (nd+1)
        do j = -(nd+1), (nd+1)
          if(abs(j).eq.(nd+1).or.abs(i).eq.(nd+1)) then
            if(abs(j)*abs(i).eq.((nd+1)*(nd+1))) then 
              mask(i,j) = 3
            else
              if(abs(i).eq.(nd+1).and.abs(j).le.(nd+1)/2)  then 
                mask(i,j) = 11 
              else if(abs(i).eq.(nd+1).and.abs(j).gt.(nd+1)/2) then
                mask(i,j) = 12
              end if
              if(abs(j).eq.(nd+1).and.abs(i).le.(nd+1)/2) then
                mask(i,j) = 22
              else if (abs(j).eq.(nd+1).and.abs(i).gt.(nd+1)/2) then
                mask(i,j) = 21
              end if 
            end if
          else
            mask(i,j) = 0
          end if
        end do
      end do

      return
      end

      double precision function q_s2n(qs,ps)
      double precision qs, ps
      q_s2n = qs / (qs**2 + ps**2)
      end 

      double precision function p_s2n(qs,ps)
      double precision qs, ps
      p_s2n = - ps / (qs**2 + ps**2)
      end 

      double precision function cubic(yp,yp1,yp2,yp3,yp4,f1,f2,f3,f4)
      double precision yp,yp1,yp2,yp3,yp4,f1,f2,f3,f4
      double precision c1, c2, c3, c4

      c1 = (yp  - yp2) * (yp  - yp3) * (yp  - yp4) / (
     &     (yp1 - yp2) * (yp1 - yp3) * (yp1 - yp4) )


      c2 = (yp  - yp1) * (yp  - yp3) * (yp  - yp4) / (
     &     (yp2 - yp1) * (yp2 - yp3) * (yp2 - yp4) )
 

      c3 = (yp  - yp1) * (yp  - yp2) * (yp  - yp4) / (
     &     (yp3 - yp1) * (yp3 - yp2) * (yp3 - yp4) )

      c4 = (yp  - yp1) * (yp  - yp2) * (yp  - yp3) / (
     &     (yp4 - yp1) * (yp4 - yp2) * (yp4 - yp3) )

      cubic = c1 * f1 + c2 * f2 + c3 * f3 + c4 * f4

      end
      subroutine norm(gn, x, q, p, nd, nx, dd, dx, time, dumps, ll, mm)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer          nd, nx
      double precision dd, dx, time
      integer          dumps, ll, mm

      integer          i, k, l, D
      double precision fg ,ex, nu, diff, l1, l2, linf
      integer          mk, ml, mi

C     Calculate the global l1, l2 and linfinity norms.

      l1 = 0.
      l2 = 0.
      linf = 0.

      mk = 0
      ml = 0
      mi = 0

      D = 0
      do i = 0, nx
        do l = -nd, nd
          do k = -nd, nd
           ex = fg(ll,mm,time,q(k),p(l),x(i),D)
           nu = gn(k,l,i,D) 
           diff = abs(ex - nu)
           l1 = l1 + diff
           l2 = l2 + diff * diff
           if (diff.ge.linf) then
               mk = k
               mi = i
               ml = l
               linf = diff
           end if
         end do 
        end do 
      end do

      l2 = sqrt(l2) / sqrt(dble(nx * nd * nd))
      l1 = l1 / (dble(nx * nd * nd))

      if (dumps .gt. 0) then
        open(14, file = 'dat/norml1.dat', status = 'unknown')
        open(15, file = 'dat/norml2.dat', status = 'unknown')
        open(16, file = 'dat/normli.dat', status = 'unknown')
        write (14,10) time, l1
        write (15,10) time, l2
        write (16,20) time, linf, mk, ml, mi
        call flush (14)
        call flush (15)
        call flush (16)
      end if

   10 format(2E23.12E3)
   20 format(2E23.12E3,3I5)

      return
      end