c fishpk40 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     program to illustrate the use of subroutine poistg to
c     solve the equation
c
c     (1/cos(x))(d/dx)(cos(x)(du/dx)) + (d/dy)(du/dy) =
c
c           2*y**2*(6-y**2)*sin(x)
c
c     on the rectangle -pi/2 .lt. x .lt. pi/2 and
c     0 .lt. y .lt. 1 with the boundary conditions
c
c     (du/dx) (-pi/2,y) = (du/dx)(pi/2,y) = 0 , 0 .le. y .le. 1  (2)
c
c     u(x,0) = 0                                           (3)
c                                 -pi/2 .le. x .le. pi/2
c     (du/dy)(x,1) = 4sin(x)                               (4)
c
c     using finite differences on a staggered grid with
c     deltax (= dx) = pi/40 and deltay (= dy) = 1/20 .
c        to set up the finite difference equations we define
c     the grid points
c
c     x(i) = -pi/2 + (i-0.5)dx            i=1,2,...,40
c
c     y(j) = (j-o.5)dy                    j=1,2,...,20
c
c     and let v(i,j) be an approximation to u(x(i),y(j)).
c     numbering the grid points in this fashion gives the set
c     of unknowns as v(i,j) for i=1,2,...,40 and j=1,2,...,20.
c     hence, in the program m = 40 and n = 20.  at the interior
c     grid point (x(i),y(j)), we replace all derivatives in
c     equation (1) by second order central finite differences,
c     multiply by dy**2, and collect coefficients of v(i,j) to
c     get the finite difference equation
c
c     a(i)v(i-1,j) + b(i)v(i,j) + c(i)v(i+1,j)
c
c     + v(i,j-1) - 2v(i,j) + v(i,j+1) = f(i,j)            (5)
c
c     where s = (dy/dx)**2, and for i=2,3,...,39
c
c     a(i) = s*cos(x(i)-dx/2)
c
c     b(i) = -s*(cos(x(i)-dx/2)+cos(x(i)+dx/2))
c
c     c(i) = s*cos(x(i)+dx/2)
c
c     f(i,j) = 2dy**2*y(j)**2*(6-y(j)**2)*sin(x(i)) , j=1,2,...,19.
c
c        to obtain equations for i = 1, we replace equation (2)
c     by the second order approximation
c
c     (v(1,j)-v(0,j))/dx = 0
c
c     and use this equation to eliminate v(0,j) in equation (5)
c     to arrive at the equation
c
c     b(1)v(1,j) + c(1)v(2,j) + v(1,j-1) - 2v(1,j) + v(1,j+1)
c
c                       = f(1,j)
c
c     where
c
c     b(1) = -s*(cos(x(1)-dx/2)+cos(x(1)+dx/2))
c
c     c(1) = -b(1)
c
c     for completeness, we set a(1) = 0.
c        to obtain equations for i = 40, we replace the derivative
c     in equation (2) at x=pi/2 in a similar fashion, use this
c     equation to eliminate the virtual unknown v(41,j) in equation
c     (5) and arrive at the equation
c
c     a(40)v(39,j) + b(40)v(40,j)
c
c     + v(40,j-1) - 2v(40,j) + v(40,j+1) = f(40,j)
c
c     where
c
c     a(40) = -b(40) = -s*(cos(x(40)-dx/2)+cos(x(40)+dx/2))
c
c     for completeness, we set c(40) = 0.  hence, in the
c     program mperod = 1.
c        for j = 1, we replace equation (3) by the second order
c     approximation
c
c                (v(i,0) + v(i,1))/2 = 0
c
c     to arrive at the condition
c
c                v(i,0) = -v(i,1) .
c
c     for j = 20, we replace equation (4) by the second order
c     approximation
c
c                (v(i,21) - v(i,20))/dy = 4*sin(x)
c
c     and combine this equation with equation (5) to arrive at
c     the equation
c
c     a(i)v(i-1,20) + b(i)v(i,20) + c(i)v(i+1,20)
c
c     + v(i,19) - 2v(i,20) + v(i,21) = f(i,20)
c
c     where
c
c     v(i,21) = v(i,20)  and
c
c     f(i,20) = 2*dy**2*y(j)**2*(6-y(j)**2)*sin(x(i)) - 4*dy*sin(x(i))
c
c     hence, in the program nperod = 2 .
c        the exact solution to this problem is
c
c        u(x,y) = y**4*cos(x) .
c
      dimension       f(42,20)   ,a(40)      ,b(40)      ,c(40)      ,
     1                w(600)     ,x(40)      ,y(20)
c
c     from dimension statement we get value of idimf = 42.  also
c     note that w has been dimensioned
c     9m + 4n + m(int(log2(n))) = 360 + 80 + 160 = 600 .
c
      idimf = 42
      mperod = 1
      m = 40
      pi = pimach(dum)
      dx = pi/float(m)
      nperod = 2
      n = 20
      dy = 1./float(n)
c
c     generate and store grid points for computation.
c
      do 101 i=1,m
         x(i) = -pi/2.+(float(i)-0.5)*dx
  101 continue
      do 102 j=1,n
         y(j) = (float(j)-0.5)*dy
  102 continue
c
c     generate coefficients .
c
      s = (dy/dx)**2
      a(1) = 0.
      b(1) = -s*cos(-pi/2.+dx)/cos(x(1))
      c(1) = -b(1)
      do 103 i=2,m
         a(i) = s*cos(x(i)-dx/2.)/cos(x(i))
         c(i) = s*cos(x(i)+dx/2.)/cos(x(i))
         b(i) = -(a(i)+c(i))
  103 continue
      a(40) = -b(40)
      c(40) = 0.
c
c     generate right side of equation.
c
      do 105 i=1,m
         do 104 j=1,n
            f(i,j) = 2.*dy**2*y(j)**2*(6.-y(j)**2)*sin(x(i))
  104    continue
  105 continue
      do 106 i=1,m
         f(i,n) = f(i,n)-4.*dy*sin(x(i))
  106 continue
      call poistg (nperod,n,mperod,m,a,b,c,idimf,f,ierror,w)
c
c     compute discretization error.  the exact solution is
c
c     u(x,y) = y**4*sin(x)
c
      err = 0.
      do 108 i=1,m
         do 107 j=1,n
            t = abs(f(i,j)-y(j)**4*sin(x(i)))
            if (t .gt. err) err = t
  107    continue
  108 continue
      print 1001 , ierror,err,w(1)
      stop
c
 1001 format (1h1,20x,25hsubroutine poistg example///
     1        10x,46hthe output from the ncar control data 7600 was//
     2        32x,10hierror = 0/
     3        18x,34hdiscretization error = 5.64171e-04/
     4        12x,32hrequired length of w array = 560//
     5        10x,32hthe output from your computer is//
     6        32x,8hierror =,i2/18x,22hdiscretization error =,e12.5/
     7        12x,28hrequired length of w array =,f4.0)
c
      end
c     compute discretization error.  the exact solution is