c fishpk36 from portlib                                  12/30/83
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                        f i s h p a k                          *
c     *                                                               *
c     *                                                               *
c     *     a package of fortran subprograms for the solution of      *
c     *                                                               *
c     *      separable elliptic partial differential equations        *
c     *                                                               *
c     *                  (version 3.1 , october 1980)                 *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *        john adams, paul swarztrauber and roland sweet         *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the national center for atmospheric research          *
c     *                                                               *
c     *                boulder, colorado  (80307)  u.s.a.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the national science foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c
c
c     example showing the use of sepx4 to solve the elliptic pde
c     (x+1)**2*uxx+2*(x+1)*ux+uyy-x*u=g(x,y) on the region
c     0.le.x.le.1, 0.le.y.le.1  with specified boundary conditions
c     at y=0,1 and mixed boundary conditions of the form
c     ux(0,y)+u(0,y), ux(1,y)+u(1,y) at x=0,1.
c     the approximation is generated on a uniform 33 by 33 grid.
c     the exact solution u(x,y)=(x*y)**3+1 is used to set the
c     right hand side, boundary conditions, and compute  second and
c     fourth order discretization error
c     the exact work space length required is 1005 words.  this
c     was determined by a previous call to sepx4 and print out of
c     w(1).
c
      dimension       usol(33,33),grhs(33,33),bda(33)    ,bdb(33)    ,
     1                w(1024)
      external cofx4
c
c     define arithmetic functions giving exact solution
c
      ue(s,t)=(s*t)**3+1.0
      uxe(s,t)=3.0*s**2*t**3
      uxxe(s,t)=6.0*s*t**3
      uye(s,t)=3.0*s**3*t**2
      uyye(s,t)=6.0*s**3*t
c
c     set limits on region
c
      a = 0.0
      b = 1.0
      c = 0.0
      d = 1.0
c
c     set grid size
c
      m = 32
      n = 32
      dlx = (b-a)/float(m)
      dly = (d-c)/float(n)
      nx = m+1
      ny = n+1
          do 102 i=1,nx
          x = a+float(i-1)*dlx
c
c     set specified boundary conditions at y=c,d
c
          usol(i,1) = ue(x,c)
          usol(i,ny) = ue(x,d)
          call cofx4 (x,af,bf,cf)
              do 101 j=1,ny
              y = c+float(j-1)*dly
c
c     set right hand side
c
              grhs(i,j) = af*uxxe(x,y)+bf*uxe(x,y)+cf*ue(x,y)+uyye(x,y)
  101         continue
  102     continue
c
c     set mixed boundary conditions at x=a,b
c
      alpha = 1.0
      beta = 1.0
          do 103 j=1,ny
          y = c+float(j-1)*dly
          bda(j) = uxe(a,y)+alpha*ue(a,y)
          bdb(j) = uxe(b,y)+beta*ue(b,y)
  103     continue
c
c     set boundary swithces
c
      mbdcnd = 3
      nbdcnd = 1
c
c     set first dimension of usol,grhs and work space length
c
      idmn = 33
      w(1) = 1024.
c
c     obtain second order approximation
c
      iorder = 2
      call sepx4 (iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,
     1            dum,dum,cofx4,grhs,usol,idmn,w,pertrb,ierror)
c
c     compute second order discretization error (relative)
c     also reset specified boundaries and right hand side.
c
      err = 0.0
          do 105 i=1,nx
          x = a+float(i-1)*dlx
          usol(i,1) = ue(x,c)
          usol(i,ny) = ue(x,d)
          call cofx4 (x,af,bf,cf)
              do 104 j=1,ny
              y = c+float(j-1)*dly
              err = amax1(err,abs((usol(i,j)-ue(x,y))/ue(x,y)))
              grhs(i,j) = af*uxxe(x,y)+bf*uxe(x,y)+cf*ue(x,y)+uyye(x,y)
  104         continue
  105     continue
      err2=err
c
c     obtain fourth order approximation
c
      iorder = 4
      call sepx4 (iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,
     1            dum,dum,cofx4,grhs,usol,idmn,w,pertrb,ierror)
c
c     compute fourth order discretization error (relative)
c
      err = 0.0
          do 107 j=1,ny
          y = c+float(j-1)*dly
              do 106 i=1,nx
              x = a+float(i-1)*dlx
              err = amax1(err,abs((usol(i,j)-ue(x,y))/ue(x,y)))
  106         continue
  107     continue
      err4=err
      iw = int(w(1)+0.5)
      print 1001,ierror,err2,err4,iw
 1001 format(1h1,20x,25hsubroutine sepx4  example ///
     120x,46hthe output from the ncar control data 7600 was //
     220x,10hierror = 0 /
     320x,48hsecond order discretization error =  1.5985e-04  /
     420x,48hfourth order discretization error =  1.85749e-06  /
     520x,33hrequired length of w array = 1024 //
     620x, 32hthe output from your computer is //
     720x, 8hierror = i2 /
     820x,36hsecond order discretization error = e12.5 /
     920x,36hfourth order discretization error = e12.5 /
     920x,29hrequired length of w array = i5)
c
      end