c fishpk39 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     program to illustrate the use of subroutine blktri to
c     solve the equation
c
c     .5/s*(d/ds)(.5/s*du/ds)+.5/t*(d/dt)(.5/t*du/dt)
c                                                          (1)
c                   = 15/4*s*t*(s**4+t**4)
c
c     on the rectangle 0 .lt. s .lt. 1 and 0 .lt. t .lt. 1
c     with the boundary conditions
c
c     u(0,t) = 0
c                            0 .le. t .le. 1
c     u(1,t) = t**5
c
c     and
c
c     u(s,0) = 0
c                            0 .le. s .le. 1
c     u(s,1) = s**5
c
c     the exact solution of this problem is u(s,t) = (s*t)**5
c
c     define the integers m = 50 and n = 63. then define the
c     grid increments deltas = 1/(m+1) and deltat = 1/(n+1).
c
c     the grid is then given by s(i) = i*deltas for i = 1,...,m
c     and t(j) = j*deltat for j = 1,...,n.
c
c     the approximate solution is given as the solution to
c     the following finite difference approximation of equation (1).
c
c     .5/(s(i)*deltas)*((u(i+1,j)-u(i,j))/(2*s(i+.5)*deltas)
c                     -(u(i,j)-u(i-1,j))/(2*s(i-.5)*deltas))
c     +.5/(t(i)*deltat)*((u(i,j+1)-u(i,j))/(2*t(i+.5)*deltat) (2)
c                     -(u(i,j)-u(i,j-1))/(2*t(i-.5)*deltat))
c               = 15/4*s(i)*t(j)*(s(i)**4+t(j)**4)
c
c             where s(i+.5) = .5*(s(i+1)+s(i))
c                   s(i-.5) = .5*(s(i)+s(i-1))
c                   t(i+.5) = .5*(t(i+1)+t(i))
c                   t(i-.5) = .5*(t(i)+t(i-1))
c
c     the approach is to write equation (2) in the form
c
c     am(i)*u(i-1,j)+bm(i,j)*u(i,j)+cm(i)*u(i+1,j)
c       +an(j)*u(i,j-1)+bn(j)*u(i,j)+cn(j)*u(i,j+1)      (3)
c           = y(i,j)
c
c     and then call subroutine blktri to determine u(i,j)
c
c
c
      dimension       y(75,105)  ,am(75)     ,bm(75)     ,cm(75)     ,
     1                an(105)    ,bn(105)    ,cn(105)    ,w(823)     ,
     2                s(75)      ,t(105)
c
      iflg = 0
      np = 1
      n = 63
      mp = 1
      m = 50
      idimy = 75
c
c     generate and store grid points for the purpose of computing the
c     coefficients and the array y.
c
      deltas = 1./float(m+1)
      do 101 i=1,m
         s(i) = float(i)*deltas
  101 continue
      deltat = 1./float(n+1)
      do 102 j=1,n
         t(j) = float(j)*deltat
  102 continue
c
c     compute the coefficients am,bm,cm corresponding to the s direction
c
      hds = deltas/2.
      tds = deltas+deltas
      do 103 i=1,m
         temp1 = 1./(s(i)*tds)
         temp2 = 1./((s(i)-hds)*tds)
         temp3 = 1./((s(i)+hds)*tds)
         am(i) = temp1*temp2
         cm(i) = temp1*temp3
         bm(i) = -(am(i)+cm(i))
  103 continue
c
c     compute the coefficients an,bn,cn corresponding to the t direction
c
      hdt = deltat/2.
      tdt = deltat+deltat
      do 104 j=1,n
         temp1 = 1./(t(j)*tdt)
         temp2 = 1./((t(j)-hdt)*tdt)
         temp3 = 1./((t(j)+hdt)*tdt)
         an(j) = temp1*temp2
         cn(j) = temp1*temp3
         bn(j) = -(an(j)+cn(j))
  104 continue
c
c     compute right side of equation
c
      do 106 j=1,n
         do 105 i=1,m
            y(i,j) = 3.75*s(i)*t(j)*(s(i)**4+t(j)**4)
  105    continue
  106 continue
c
c     the nonzero boundary conditions enter the linear system via
c     the right side y(i,j). if the equations (3) given above
c     are evaluated at i=m and j=1,...,n then the term cm(m)*u(m+1,j)
c     is known from the boundary condition to be cm(m)*t(j)**5.
c     therefore this term can be included in the right side y(m,j).
c     the same analysis applies at j=n and i=1,..,m. note that the
c     corner at j=n,i=m includes contributions from both boundaries.
c
      do 107 j=1,n
         y(m,j) = y(m,j)-cm(m)*t(j)**5
  107 continue
      do 108 i=1,m
         y(i,n) = y(i,n)-cn(n)*s(i)**5
  108 continue
c
c     determine the approximate solution u(i,j)
c
  109 call blktri (iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w)
      iflg = iflg+1
      if (iflg-1) 109,109,110
c
c     compute discretization error
c
  110 err = 0.
      do 112 j=1,n
         do 111 i=1,m
            z = abs(y(i,j)-(s(i)*t(j))**5)
            if (z .gt. err) err = z
  111    continue
  112 continue
      iw = int(w(1))
      print 1001 , ierror,err,iw
      stop
c
 1001 format (1h1,20x,25hsubroutine blktri example///
     1        10x,46hthe output from the ncar control data 7600 was//
     2        32x,10hierror = 0/
     3        18x,34hdiscretization error = 1.64478e-05/
     4        12x,32hrequired length of w array = 823//
     5        10x,32hthe output from your computer is//
     6        32x,8hierror =,i2/18x,22hdiscretization error =,e12.5/
     7        12x,28hrequired length of w array =,i4)
c
      end
c