c fishpk41 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     this program illustrates the use of the subroutine pois3d to solve
c     the equation
c
c     (d/dx)(du/dx) + (d/dy)(du/dy) + (1+z)**2*(d/dz)(du/dz)
c
c     - 2*(1+z)*(du/dz) = 2*sin(x)*sin(y)*(1+z)**4                   (1)
c
c     on the parallelepiped -pi .lt. x .lt. pi, -pi .lt. y .lt. pi,
c     0 .lt. z .lt. 1 with boundary conditions
c
c     u periodic in x
c
c     u periodic in y
c
c     (du/dz)(x,y,0) = 4*sin(x)*sin(y)    -pi .lt. x,y .lt. pi       (2)
c
c     u(x,y,1) = 16*sin(x)*sin(y)         -pi .lt. x,y .lt. pi       (3)
c
c     using a finite difference grid with deltax (= dx) = 2*pi/30 ,
c     deltay (= dy) = 2*pi/30, and deltaz (= dz) = 1/10.
c         to set up the finite difference equations we define the grid
c     points
c
c         x(i) = -pi + (i-1)*dx         i=1,2,...,31
c
c         y(j) = -pi + (j-1)*dy         j=1,2,...,31
c
c         z(k) = (k-1)*dz               k=1,2,...,11
c
c     and let v(i,j,k) be an approximation to u(x(i),y(j),z(k)).
c     numbering the grid points in this fashion gives the set of
c     unknowns as v(i,j,k) for i=1,2,...,30, j=1,2,...,30, k=1,2,...,10.
c     hence, in the program l=30, m = 30, and n = 10.  at the interior
c     grid point (x(i),y(j),z(k)), we replace all derivatives in
c     equation (1) by second order central finite differences and
c     collect coefficients of v(i,j,k) to get the finite difference
c     equation
c
c        (v(i-1,j,k) - 2v(i,j,k) + v(i+1,j,k))/dx**2
c
c      + (v(i,j-1,k) - 2v(i,j,k) + v(i,j+1,k))/dy**2
c
c      + a(k)v(i,j,k-1) + b(k)v(i,j,k) + c(k)v(i,j,k+1) = f(i,j,k)   (4)
c
c     where for k=2,3,...,9
c
c     a(k) = (1+z(k))**2/dz**2 + (1+z(k))/dz
c
c     b(k) = -2(1+z(k))**2/dz**2
c
c     c(k) = (1+z(k))**2/dz**2 - (1+z(k))/dz
c
c     f(i,j,k) = 2sin(x(i))*sin(y(j))*(1+z(k))**4  for i,j=1,2,...,30.
c
c         to obtain equations for k=1, we replace the derivative in
c     equation (2) by a second order central finite difference approx-
c     imation, use this equation to eliminate the virtual unknown
c     v(i,j,0) in equation (4) and arrive at the equation
c
c              (v(i-1,j,1) -2v(i,j,1) + v(i+1,j,1))/dx**2
c
c            + (v(i,j-1,1) -2v(i,j,1) + v(i,j+1,1))/dy**2
c
c            + b(1)v(i,j,1) + c(1)v(i,j,2) = f(i,j,1)
c
c     where
c            b(1) = -c(1) = -2(1+z(1))**2/dz**2 = -2/dz**2
c
c            f(i,j,1) = (10 + 8/dz)sin(x(i))*sin(y(j))
c
c            for i,j=1,2,...,30.  for completeness we set a(1) = 0.
c
c         to obtain equations for k=10, we incorporate equation (3) into
c     equation (4) by setting
c
c            v(i,j,11) = u(x(i),y(j),1) = 16sin(x(i))*sin(y(j))
c
c     and arrive at the equation
c
c              (v(i-1,j,10) - 2v(i,j,10) + v(i+1,j,10))/dx**2
c
c            + (v(i,j-1,10) - 2v(i,j,10) + v(i,j+1,10))/dy**2
c
c            + a(10)v(i,j,9) + b(10)v(i,j,10) = f(i,j,10)
c
c     where
c
c            a(10) = (1+z(10))**2/dz**2 + (1+z(10))/dz
c
c            b(10) = -2(1+z(10))**2/dz**2
c
c            f(i,j,10) = 2sin(x(i))*sin(y(j))*((1+z(10))**4
c                        -8*((1+z(10))**2/dz**2 - (1+z(10))/dz))
c
c                        for i,j=1,2,...,30.
c
c     for completeness, we set c(10) = 0.  hence, in the program,
c     nperod = 1.
c         the periodicity conditions on u give the conditions
c
c             v(0,j,k) = v(30,j,k) and v(31,j,k) = v(1,j,k)
c                        for j=1,2,...,30 and k=1,2,...,10,
c             and
c             v(i,0,k) = v(i,30,k) and v(i,31,k) = v(i,1,k)
c                        for i=1,2,...,30 and k=1,2,...,10.
c
c     hence, in the program lperod = mperod = 0.
c
      dimension       f(32,33,10),a(10)      ,b(10)      ,c(10)      ,
     1                w(350)     ,x(30)      ,y(30)      ,z(10)
c
c     from the dimension statement we get that ldimf = 32, mdimf = 33,
c     and note that w has been dimensioned according to its description.
c
      ldimf = 32
      mdimf = 33
      pi = pimach(dum)
      lperod = 0
      l = 30
      dx = 2.*pi/float(l)
      c1 = 1./dx**2
      mperod = 0
      m = 30
      dy = 2.*pi/float(m)
      c2 = 1./dy**2
      nperod = 1
      n = 10
      dz = 1./float(n)
      dzsq = 1./dz**2
c
c     generate grid points for later use.
c
      do 101 i=1,l
         x(i) = -pi+float(i-1)*dx
  101 continue
      do 102 j=1,m
         y(j) = -pi+float(j-1)*dy
  102 continue
c
c     generate coefficients
c
      a(1) = 0.
      b(1) = -2.*dzsq
      c(1) = -b(1)
      z(1) = 0.
      do 103 k=2,n
         z(k) = float(k-1)*dz
         t = 1.+z(k)
         a(k) = t**2*dzsq+t/dz
         b(k) = -2.*t**2*dzsq
         c(k) = t**2*dzsq-t/dz
  103 continue
c
c     generate right side of equation
c
      do 106 i=1,l
         do 105 j=1,m
            do 104 k=2,n
               f(i,j,k) = 2.*sin(x(i))*sin(y(j))*(1.+z(k))**4
  104       continue
  105    continue
  106 continue
      do 108 i=1,l
         do 107 j=1,l
            f(i,j,1) = (10.+8./dz)*sin(x(i))*sin(y(j))
            f(i,j,n) = f(i,j,n)-c(n)*16.*sin(x(i))*sin(y(j))
  107    continue
  108 continue
      c(n) = 0.
c
c     call pois3d to solve equations.
c
      call pois3d (lperod,l,c1,mperod,m,c2,nperod,n,a,b,c,ldimf,mdimf,
     1             f,ierror,w)
c
c     compute discretization error.  the exact solution is
c
c              u(x,y,z) = sin(x)*sin(y)*(1+z)**4
c
      err = 0.
      do 111 i=1,l
         do 110 j=1,m
            do 109 k=1,n
               t = abs(f(i,j,k)-sin(x(i))*sin(y(j))*(1.+z(k))**4)
               if (t .gt. err) err = t
  109       continue
  110    continue
  111 continue
      print 1001 , ierror,err
      stop
c
 1001 format (1h1,20x,25hsubroutine pois3d example///
     1        10x,46hthe output from the ncar control data 7600 was//
     2        32x,10hierror = 0/
     3        18x,34hdiscretization error = 2.93277e-02//
     4        10x,32hthe output from your computer is//
     5        32x,8hierror =,i2/18x,22hdiscretization error =,e12.5)
c
      end
c     compute discretization error.  the exact solution is