subroutine charac ( data, m, ifail )

c*********************************************************************72
c
c  this routine advances one grid step the solution of the system
c
c    a1 * u(x) + a2 * u(y) + a3 * v(x) + a4 * v(y) = h1
c    b1 * u(x) + b2 * u(y) + b3 * v(x) + b4 * v(y) = h2
c
c  where u(x) means partial derivative of u with respect to x, and
c  so on, and a1 = a1(x,y,u,v), and so on.
c
c  the initial data is given in the matrix "data", each column of
c  four elements containing a value of x, y, u, v.
c
c  m is the number of data points on the initial curve.
c
      implicit none

      integer m

      real d0(4)
      real d1(4)
      real d2(4)
      real d3(4)
      real d4(4)
      real data(4,m)
      real dx
      real dx1
      real dx2
      real dy
      integer fail
      integer flag
      integer ifail
      integer i
      integer j
      real s1(4)
      real s2(4)
      real s3(4)
      real t1(4)
      real t2(4)
      real t3(4)
      real t4(4)
      real t5(4)
      real t6(4)
      real t7(4)
      real v1(8)
      real v2(8)
      real v3(8)
      real v4(8)
      real v5(8)
      real v6(8)
      real v7(8)
      real v8(8)
      real v9(8)
      real v10(8)

      common / chfail / fail, flag

      fail = 0
      m = m - 4
      if ( m .le. 0 ) return

      do j = 1, 4
        d1(j) = data(j,1)
        d2(j) = data(j,2)
        d3(j) = data(j,3)
        d4(j) = data(j,4)
      end do

      call chvar ( d1, v2 )
      if ( ifail .ne. 0 ) go to 250
      call chvar ( d2, v3 )
      if ( ifail .ne. 0 ) go to 250
      call chvar ( d3, v4 )
      if ( ifail .ne. 0 ) go to 250
      call chvar ( d4, v5 )
      if ( ifail .ne. 0 ) go to 250

      flag = 0
      dy = d2(2) - d1(2)
      dx = d2(1) - d1(1)
      if ( dy ) 148, 149, 150
148   dy = -dy
      dx = -dx
      go to 150
149   dy = 1.0e+00
      dx = 1.0e+30
150   dx1 = dy / v2(1)
      dx2 = dy / v2(2)
      if ( dx1 .lt. dx .and. dx2 .ge. dx ) go to 170
      if ( dx2 .lt. dx .and. dx1 .ge. dx ) go to 175
      if ( dx2 .ge. dx1 ) go to 175

170   flag = 1
175   continue

      call chstep ( d1, v2, d3, v4, t1 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t1, v6 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( d2, v3, d4, v5, t6 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t6, v10 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( d3, v4, d4, v5, t2 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t2, v7 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( d2, v3, d3, v4, t4 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t4, v9 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( t4, v9, t2, v7, t3 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t3, v8 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( d1, v2, d2, v3, t5 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t5, v1 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( t5, v1, t4, v9, t4 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t4, v1 )
      if ( fail .ne. 0 ) go to 250

      call chstep ( t4, v1, t3, v8, t4 )
      if ( fail .ne. 0 ) go to 250
      call chvar ( t4, v9 )
      if ( fail .ne. 0 ) go to 250

      do i = 1, m

        do j = 1, 8
          v1(j) = v2(j)
          v2(j) = v3(j)
          v3(j) = v4(j)
          v4(j) = v5(j)
        end do

        do j = 1, 4
          d0(j) = d1(j)
          d1(j) = d2(j)
          d2(j) = d3(j)
          d3(j) = d4(j)
          d4(j) = data(j,i+4)
        end do

        call chvar ( d4, v5 )
        if ( fail .ne. 0 ) go to 250

        call chstep ( d0, v1, d4, v5, s1 )
        if ( fail .ne. 0 ) go to 250
        call chstep ( d2, v3, d4, v5, t5 )
        if ( fail .ne. 0 ) go to 250

        call chvar ( t5, v1 )
        if ( fail .ne. 0 ) go to 250

        call chstep ( t1, v6, t5, v1, s2 )
        if ( fail .ne. 0 ) go to 250
        call chstep ( d3, v4, d4, v5, t1 )
        if ( fail .ne. 0 ) go to 250

        call chvar ( t1, v6 )
        if ( fail .ne. 0 ) go to 250

        call chstep ( t2, v7, t1, v6, t7 )
        if ( fail .ne. 0 ) go to 250

        do j = 1, 4
          t2(j)    = t1(j)
          v7(j)    = v6(j)
          v7(j+4)  = v6(j+4)
          t1(j)    = t6(j)
          v6(j)    = v10(j)
          v6(j+4)  = v10(j+4)
          t6(j)    = t5(j)
          v10(j)   = v1(j)
          v10(j+4) = v1(j+4)
        end do

        call chvar ( t7, v1 )
        if ( fail .ne. 0 ) go to 250

        call chstep ( t3, v8, t7, v1, t5 )
        if ( fail .ne. 0 ) go to 250

        do j = 1, 4
          t3(j) = t7(j)
          v8(j) = v1(j)
          v8(j+4) = v1(j+4)
        end do

        call chvar ( t5, v1 )
        if ( fail .ne. 0 ) go to 250

        call chstep ( t4, v9, t5, v1, s3 )
        if ( fail .ne. 0 ) go to 250

        do j = 1, 4
          t4(j) = t5(j)
          v9(j) = v1(j)
          v9(j+4) = v1(j+4)
        end do
c
c  extrapolate the three successive approximations.
c
        do j = 1, 4
          data(j,i) = 0.3333333333e+00 
     &      * ( 8.0e+00 * s3(j) - 6.0e+00 * s2(j) + s1(j) )
        end do

      end do

      ifail = 0
      return
c
c  error exit.
c
250   ifail = fail
      return
      end
      subroutine chvar ( xyuv, var )

c*********************************************************************72
c
c  computes the values sigma1, sigma2, r1, r2, s1, s2, t1, t2,
c  stored in the list var, from the coefficient functions and
c  the values x, y, u, v in the list xyuv.
c
      implicit none

      real a
      real b
      real c
      real d
      integer i
      integer fail
      integer flag
      real t(10)
      real var(8)
      real xyuv(4)

      common / chfail / fail, flag

      call chcoef ( t, xyuv )
c
c  coefficients of system are stored in the list t.
c
      a = t(1) * t(8) - t(3) * t(6)
      b = 0.5e+00 
     &  * ( t(1) * t(9) - t(4) * t(6) - t(3) * t(7) + t(2) * t(8) )
      c = t(2) * t(9) - t(4) * t(7)

      if ( a .ne. 0.0e+00 ) go to 150
      if ( b .eq. 0.0e+00 ) go to 500
      var(1) = 1.0e+15
      var(2) = 0.5e+00 * c / b
      if ( b .gt. 0.0e+00 ) go to 400
      var(1) = var(2)
      var(2) = 1.0e+15
      go to 400

150   d = b * b - a * c
      if ( d .le. 0.0e+00 ) go to 500
      d = sqrt ( d )
      if ( b .lt. 0.0e+00 ) go to 300
      var(1) = ( b + d ) / a
      var(2) = c / ( a * var(1) )
      go to 400

300   var(2) = ( b - d ) / a
      var(1) = c / ( a * var(2) )

400   do i = 1, 2
        t(4) = t(1) * var(i) - t(2)
        t(9) = t(6) * var(i) - t(7)
        var(i+2) = t(1) * t(9) - t(6) * t(4)
        var(i+4) = t(3) * t(9) - t(8) * t(4)
        var(i+6) = t(5) * t(9) - t(10) * t(4)
      end do

      return
c
c  error exit.
c
500   fail = 1
      return
      end
      subroutine chstep ( d1, v1, d2, v2, d3 )

c*********************************************************************72
c
c  this routine computes the values of x3, y3, u3, v3 at the point
c  determined by the intersection of the characteristics through (x1,y1)
c  and (x2,y2).
c
      implicit none

      real aa
      real bb
      real d1(4)
      real d2(4)
      real d3(4)
      real d31
      real d32
      real dem1
      real dem2
      real eps
      integer fail
      integer flag
      real r1
      real r2
      real rho
      real s1
      real s2
      real sig
      real t1
      real t2
      real ta
      real tb
      real tc
      real td
      real v1(8)
      real v2(8)

      common / chfail / fail, flag

      eps = 1.0e-05
      if ( flag .ne. 0 ) go to 150
      sig = v1(2)
      r1 = v1(4)
      s1 = v1(6)
      t1 = v1(8)
      go to 180

150   sig = v1(1)
      r1 = v1(3)
      s1 = v1(5)
      t1 = v1(7)

180   continue

      if ( flag .ne. 0 ) go to 250

      rho = v2(1)
      r2 = v2(3)
      s2 = v2(5)
      t2 = v2(7)
      go to 280

250   rho = v2(2)
      r2 = v2(4)
      s2 = v2(6)
      t2 = v2(8)

280   continue
c
c  compute x3, y3, u3, v3.
c
      dem1 = sig - rho
      if ( abs ( dem1 ) .lt. eps * abs ( sig ) ) go to 900
      aa = d1(2) - sig * d1(1)
      bb = d2(2) - rho * d2(1)
      d32 = ( sig * bb - rho * aa ) / dem1
      d31 = ( bb - aa ) / dem1
      ta = t1 * ( d31 - d1(1) ) + r1 * d1(3) + s1 * d1(4)
      tb = t2 * ( d31 - d2(1) ) + r2 * d2(3) + s2 * d2(4)
      tc = r1 * s2
      td = r2 * s1
      dem2 = tc - td
      tc = amax1 ( abs ( tc ), abs ( td ) )
      if ( abs ( dem2 ) .le. eps * tc ) go to 900
      d3(3) = ( ta * s2 - tb * s1 ) / dem2
      d3(4) = ( tb * r1 - ta * r2 ) / dem2
      d3(2) = d32
      d3(1) = d31
      return

900   fail = 2
      return
      end