c fishpk37 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     an example showing the use of sepeli to solve the separable
c     elliptic partial differential equation . . .
c     (x+1)**2*uxx+2*(x+1)*ux+exp(y)*uyy-(x+y)*u = g(x,y) on
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 1118 words.
c     this was determined by a previous call to sepeli  and print
c     out of w(1).
c
      dimension       usol(33,33),grhs(33,33),bda(33)    ,bdb(33)    ,
     1                w(1118)
c
c     declare coefficient subroutines external
c
      external        cofx       ,cofy
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 cofx (x,af,bf,cf)
         do 101 j=1,ny
            y = c+float(j-1)*dly
            call cofy (y,df,ef,ff)
c
c     set right hand side
c
            grhs(i,j) = af*uxxe(x,y)+bf*uxe(x,y)+cf*ue(x,y)+
     1                  df*uyye(x,y)+ef*uye(x,y)+ff*ue(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) = 1118.
c
c     set work space length in first word
c     set initial call parameter to zero
c
      intl = 0
c
c     obtain second order approximation
c
      iorder = 2
      call sepeli (intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,
     1             nbdcnd,dum,dum,dum,dum,cofx,cofy,grhs,usol,idmn,w,
     2             pertrb,ierror)
      err = 0.0
      do 105 i=1,nx
         x = a+float(i-1)*dlx
         do 104 j=1,ny
            y = c+float(j-1)*dly
            err = amax1(err,abs((usol(i,j)-ue(x,y))/ue(x,y)))
  104    continue
  105 continue
      err2 = err
c
c     obtain fourth order approximation
c
      iorder = 4
c
c     non-initial call
c
      intl = 1
      call sepeli (intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,
     1             nbdcnd,dum,dum,dum,dum,cofx,cofy,grhs,usol,idmn,w,
     2             pertrb,ierror)
c
c     compute discretization error
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))
      print 1001 , ierror,err2,err4,iw
c
c
 1001 format (1h1,20x,25hsubroutine sepeli example///
     1        20x,46hthe output from the ncar control data 7600 was//
     2        20x,10hierror = 0/
     3        20x,47hsecond order discretization error = 9.78910e-05/
     4        20x,47hfourth order discretization error = 1.47351e-06/
     5        20x,33hrequired length of w array = 1118//
     6        20x,32hthe output from your computer is//
     7        20x,8hierror =i2/
     8        20x,36hsecond order discretization error = e12.5/
     9        20x,36hfourth order discretization error = e12.5/
     +        20x,29hrequired length of w array = i4)
c
      end