c--#include <fintrf.h>
#define MW_POINTER_TYPE integer*4
#define MW_INT_TYPE integer*4

#define DBLE_PRECISION 
c--#define F77_TO_F95
#define VERBOSE 1

#ifdef DBLE_PRECISION  
# ifdef F77_TO_F95
#  define INT_TYPE integer(kind=4)
#  define REAL_TYPE real(kind=8)
#  define COMPLEX_TYPE complex(kind=8)
# else
#  define INT_TYPE integer*4
#  define REAL_TYPE real*8
#  define COMPLEX_TYPE complex*16
# endif
# define CMPLX dcmplx
# define REAL_PART dble
# define IMAG_PART dimag
# define CABS cdabs
# define CSQRT cdsqrt
#else
# define REAL_PART real
# define IMAG_PART aimag
#endif

!
! purpose:  mexsepeli solves for the 2nd order finite difference
!           approximation to the separable elliptic equation arising
!           from conformal mapping and grid generation.
!
! The calling sequence from MATLAB should be
!
!   >> [x, y, perturb, ierror] = mexsepeli ( x,    y, ...
!                                            l2,   m2, ...
!                                            seta, sxi, ...
!                                            nwrk, nx2 );
!
! output:  x,y:  represent "x" and "y" grid coordinates
!      perturb:  optional (don't think this is needed for our 
!                                                 application)
!       ierror:  optional. An error flag that indicates invalid 
!                input parameters or failure to find a solution.
!                          = 0 no error
!                          = 4 if attempt to find a solution fails.
!                              (the linear system generated is not
!                              diagonally dominant.)
!
! input:   x,y:  represent "x" and "y" grid coordinates
!           l2:  x independent variable ranges from 0 to l2
!           m2:  y independent variable ranges from 0 to m2
!         seta:  digitized points on boundaries 1,3
!          sxi:  digitized points on boundaries 2,4
!
!    so   1)  nlhs = 2
!         2)  plhs(1) = x
!             plhs(2) = y
!         3)  nrhs = 6 
!         4)  prhs(1) = x 
!             prhs(2) = y 
!             prhs(3) = l2
!             prhs(4) = m2
!             prhs(5) = seta
!             prhs(6) = sxi
!
! 
! $Id: mexsepeli.f,v 1.2 1997/10/03 19:43:04 jevans Exp jevans $
! Currently locked by $Locker: jevans $ (not locked if blank)
!  (Time in GMT, EST=GMT-5:00
!     $Log: mexsepeli.f,v $
! Revision 1.2  1997/10/03  19:43:04  jevans  Had to change compile
! time options to include "-align dcommons" and "-r8".
!
! the "mexCopyReal8ToPtr" and related functions suffer from the
! fact that they assume all arrays to be a single column.  Since
! my x and y are dimensioned to be larger than their logical
! size (i.e. they are 800x800, but might logically be only 100x80
! or so), those arrays had to be columnified before passing them
! back to MATLAB, and "uncolumnified" upon passing them into
! MATLAB.
!
! Revision 1.1  1997/10/03  17:50:43  jevans
! Initial revision
!
      subroutine mexfunction ( nlhs, plhs, nrhs, prhs )
      implicit none
      MW_POINTER_TYPE plhs(*), prhs(*)
      MW_INT_TYPE nlhs, nrhs

#include "mexsepeli.h"
                             ! Buffer variables to pass into matlab
      MW_POINTER_TYPE p      ! routines. Their type should be exactly
      real*8 bffr(0:nx2*ny2) ! as specified by matlab. Note: these
      MW_INT_TYPE size       ! and only these variables are allowed
                             ! to enter in matlab calls.
      MW_INT_TYPE mxIsNumeric,             ! these are MatLab
     &            mxGetM,  mxGetN,         ! interface library 
     &            mxGetPr, mxCreateFull,   ! functions
     &            size_m,size_n

!-- internal -->

      REAL_TYPE A,B,C,D, pertrb
      INT_TYPE i,j, ierr
#ifdef VERBOSE
      write(*,*) 'Enter mexfunction'
#endif
!
! Before doing anything else, check for proper number of arguments.
! There must be at least argument on the left hand size and no more
! than 3.
!
      if (nlhs .lt. 1) then
         call mexErrMsgTxt ( 'Not enough input args for mexsepeli.' )
      elseif (nlhs .gt.3 ) then
         call mexErrMsgTxt ( 'Too many input args for mexsepeli.' )
      endif
      if (nrhs .ne. 6) then
         call mexErrMsgTxt('mexsepeli requires 6 input arguments')
      endif
!
! Check that all input arguments were numeric.
!
      if (mxIsNumeric(prhs(1)) .ne. 1) then
         call mexErrMsgTxt( 'x input array should be numeric' )
      endif
      if (mxIsNumeric(prhs(2)) .ne. 1) then
         call mexErrMsgTxt( 'y input array should be numeric' )
      endif
      if (mxIsNumeric(prhs(3)) .ne. 1) then
         call mexErrMsgTxt( 'l2 input should be numeric' )
      endif
      if ( mxIsNumeric(prhs(4)) .ne. 1 ) then
         call mexErrMsgTxt( 'm2 input should be numeric' )
      endif
      if ( mxIsNumeric(prhs(5)) .ne. 1 ) then
         call mexErrMsgTxt( 'seta input array should be numeric' )
      endif
      if ( mxIsNumeric(prhs(6)) .ne. 1 ) then
         call mexErrMsgTxt( 'sxi input array should be numeric' )
      endif
!
! Check to see that the l2 and m2 arguments were scalars
!
      size_m=mxGetM(prhs(3))
      size_n=mxGetN(prhs(3))
      if (size_n.ne.1 .or. size_m.ne.1) then
         call mexErrMsgTxt( 'l2 argument should be a scalar' )
      endif

      size_m=mxGetM(prhs(4))
      size_n=mxGetN(prhs(4))
      if (size_n.ne.1 .or. size_m.ne.1 ) then
         call mexErrMsgTxt(' m2 argument should be a scalar' )
      endif
#ifdef VERBOSE
      write(*,*) 'mexfunction: passed initial clecks'
#endif
!
! Load the data into Fortran matrices.
!
      p=mxGetPr(prhs(1))                     ! X -->
      size_m=mxGetM(prhs(1))
      size_n=mxGetN(prhs(1))
      size=size_m*size_n
      call mxCopyPtrToReal8 (p, bffr, size)
      do j=0,size_n-1
        do i=0,size_m-1
          x(i,j)=bffr(i+j*size_m)
        enddo
      enddo
#ifdef VERBOSE
      write(*,*) 'x: size_m,size_n =', size_m,size_n
#endif

      p=mxGetPr(prhs(2))                     ! Y --> 
      size_m=mxGetM(prhs(2))
      size_n=mxGetN(prhs(2))
      size=size_m*size_n
      call mxCopyPtrToReal8 (p, bffr, size)
      do j=0,size_n-1
        do i=0,size_m-1
          y(i,j)=bffr(i+j*size_m)
        enddo
      enddo
#ifdef VERBOSE
      write(*,*) 'y: size_m,size_n =', size_m,size_n
#endif
                                             ! Ranges for
      size=1                                 ! independent 
      p=mxGetPr(prhs(3))                     ! variables 
      call mxCopyPtrToReal8 (p, bffr, size)  ! 
      L2=bffr(0)                             !   L2,M2
      p=mxGetPr(prhs(4))                     !
      call mxCopyPtrToReal8 (p, bffr, size)  ! 0<x<L2, 0<y<Y2
      M2=bffr(0)
#ifdef VERBOSE
      write(*,*) 'L2,M2 =', L2,M2
#endif
      p=mxGetPr(prhs(5))                     ! SETA -->
      size_m=mxGetM(prhs(5))
      size_n=mxGetN(prhs(5))
      size=size_m*size_n
      call mxCopyPtrToReal8 (p, bffr, size)
      do i=0,size-1
        SETA(i)=bffr(i)
      enddo
      NETA=size
#ifdef VERBOSE
      write(*,*) 'NETA =', NETA 
#endif
      p=mxGetPr(prhs(6))                     ! SXI -->
      size_m=mxGetM(prhs(6))
      size_n=mxGetN(prhs(6))
      size=size_m*size_n
      call mxCopyPtrToReal8 (p, bffr, size)
      do i=0,size-1
        SXI(i)=bffr(i)
      enddo
      NXI=size
#ifdef VERBOSE
      write(*,*) 'NXI =', NXI 
#endif
      do j=0,NETA-1                          ! Initialize
        do i=0,NXI-1                         ! right hand side
          rhs(i,j)=0.
        enddo
      enddo
!
! Call the computational subroutine.
!
      A=0.
      B=float(L2)
      C=0.
      D=float(M2)
      ewrk(1) = nwrk
      call sepeli (0,2, A,B, L2,1, wrk,wrk,wrk,wrk,
     &                  C,D, M2,1, wrk,wrk,wrk,wrk,
     &             rhs, x, nx2+1, ewrk,pertrb,ierr)
      if ( ierr .eq. 1 ) then
         call mexErrMsgTxt
     &      ( 'range of independent variables is out of whack' )
      else if ( ierr .eq. 2 ) then
         call mexErrMsgTxt
     &         ( 'boundary condition mbdcnd wrongly specified' )
      else if ( ierr .eq. 3 ) then
         call mexErrMsgTxt
     &         ( 'boundary condition nbdcnd wrongly specified' )
      else if ( ierr .eq. 4 ) then
         call mexErrMsgTxt
     &  ( 'linear system generated is not diagonally dominant' )
      else if ( ierr .eq. 5 ) then
         call mexErrMsgTxt( 'idmn is too small' )
      else if ( ierr .eq. 6 ) then
         call mexErrMsgTxt( 'm is too small or too large' )
      else if ( ierr .eq. 7 ) then
         call mexErrMsgTxt( 'n is too small or too large' )
      else if ( ierr .eq. 8 ) then
         call mexErrMsgTxt( 'iorder is not 2 or 4' )
      else if ( ierr .eq. 9 ) then
         call mexErrMsgTxt( 'intl is not 0 or 1' )
      else if ( ierr .eq. 10 ) then
         call mexErrMsgTxt( 'afun*dfun less than or equal to 0' )
      else if ( ierr .eq. 11 ) then
         call mexErrMsgTxt
     &          ('work space length input in w(1) is not right' )
      end if

      A=0.
      B=float(L2)
      C=0.
      D=float(M2)
      ewrk(1) = nwrk
      call sepeli (0,2, A,B, L2,1, wrk,wrk,wrk,wrk,
     &                  C,D, M2,1, wrk,wrk,wrk,wrk,
     &             rhs, y, nx2+1, ewrk,pertrb,ierr)
      if ( ierr .eq. 1 ) then
         call mexErrMsgTxt
     &        ( 'range of independent variables is out of whack' )
      else if ( ierr .eq. 2 ) then
         call mexErrMsgTxt
     &           ( 'boundary condition mbdcnd wrongly specified' )
      else if ( ierr .eq. 3 ) then
         call mexErrMsgTxt
     &           ( 'boundary condition nbdcnd wrongly specified' )
      else if ( ierr .eq. 4 ) then
         call mexErrMsgTxt
     &    ( 'linear system generated is not diagonally dominant' )
      else if ( ierr .eq. 5 ) then
         call mexErrMsgTxt( 'idmn is too small' )
      else if ( ierr .eq. 6 ) then
         call mexErrMsgTxt( 'm is too small or too large' )
      else if ( ierr .eq. 7 ) then
         call mexErrMsgTxt( 'n is too small or too large' )
      else if ( ierr .eq. 8 ) then
         call mexErrMsgTxt( 'iorder is not 2 or 4' )
      else if ( ierr .eq. 9 ) then
         call mexErrMsgTxt( 'intl is not 0 or 1' )
      else if ( ierr .eq. 10 ) then
         call mexErrMsgTxt( 'afun*dfun less than or equal to 0' )
      else if ( ierr .eq. 11 ) then
         call mexErrMsgTxt
     &         ( 'work space length input in w(1) is not right' )
      end if
!
! Create matrices for output arguments Since the mxCopyReal8ToPtr
! routine just assumes one long column matrix, we have to "columnify"
! both x and y.
!
      size_m=mxGetM(prhs(1))             ! X -->
      size_n=mxGetN(prhs(1))
      size=size_m*size_n
      plhs(1)=mxCreateFull (size_m, size_n, 0)
      p=mxGetPr(plhs(1))
      do j=0,size_n-1
        do i=0,size_m-1
          bffr(i+j*size_m)=x(i,j)
        enddo
      enddo
      call mxCopyReal8ToPtr (bffr, p, size)

      size_m=mxGetM(prhs(2))             ! Y -->
      size_n=mxGetN(prhs(2))
      size=size_m*size_n
      plhs(2)=mxCreateFull (size_m, size_n, 0)
      p=mxGetPr(plhs(2))
      do j=0,size_n-1
        do i=0,size_m-1
          bffr(i+j*size_m) = y(i,j)
        enddo
      enddo
      call mxCopyReal8ToPtr (bffr, p, size)
      return
      end

!
!          *****   ******   ******   ******   **       **
!         **   **  **   **  **   **  **   **  **       **
!         **       **       **   **  **       **       **
!          *****   *****    **   **  *****    **       **
!              **  **       ******   **       **       **
!         **   **  **   **  **       **   **  **   **  **
!          *****   ******   **       ******   *******  **
!

      subroutine SEPELI (INTL,IORDER, A,B, M,MBDCND, BDA,ALPHA,
     &                      BDB,BETA, C,D, N,NBDCND, BDC,GAMA,
     &                  BDD,XNU,GRHS,USOL, IDMN, W, PERTRB, IERROR)
!
! dimension of    bda(n+1), bdb(n+1), bdc(m+1), bdd(m+1),
! arguments       usol(idmn,n+1),grhs(idmn,n+1), w(see argument list)
!
! latest      march 1985
! revision
!
! purpose     sepeli solves for either a second- or a fourth-order
!             finite difference approximation to a separable elliptic
!             equation
!
!                    d^2 u           du                     d^2 u
!             af(x)*------- + bf(x)*---- + cf(x)*u + df(y)*-------
!                     dx^2           dx                     dy^2
!
!                                          du
!                                 + ef(y)*---- + ff(y)*u = g(x,y)
!                                          dy
!
!
!             on a rectangle (x greater than or equal to a and less
!             than or equal to b; y greater than or equal to c and
!             less than or equal to d). Any combination of periodic
!             or mixed boundary conditions is allowed.
!
!             the possible boundary conditions are:
!             in the x-direction:
!                 (0) periodic, u(x+b-a,y)=u(x,y) for all y,x (1)
!                            u(a,y),u(b,y) are specified for all y
!                 (2) u(a,y), du(b,y)/dx+beta*u(b,y) are specified
!                                                       for all y
!                 (3) du(a,y)/dx+alpha*u(a,y),du(b,y)/dx+beta*u(b,y)
!                             are specified for all y
!                 (4) du(a,y)/dx+alpha*u(a,y),u(b,y) are specified
!                                             for all y
!
!             in the y-direction:
!                 (0) periodic, u(x,y+d-c)=u(x,y) for all x,y
!                 (1) u(x,c),u(x,d) are specified for all x
!                 (2) u(x,c),du(x,d)/dy+xnu*u(x,d) are specified
!                                                    for all x
!                 (3) du(x,c)/dy+gama*u(x,c),du(x,d)/dy+xnu*u(x,d)
!                                       are specified for all x
!                 (4) du(x,c)/dy+gama*u(x,c),u(x,d) are specified
!                                                       for all x
!
! usage       call sepeli (intl,iorder,a,b,m,mbdcnd,bda,
!                          alpha,bdb,beta,c,d,n,nbdcnd,bdc,
!                          gama,bdd,xnu,cofx,cofy,grhs,usol,
c                                     idmn,w,pertrb,ierror)
c
c arguments
c on input               intl
c                          = 0 on initial entry to sepeli or if any
c                              of the arguments c,d, n, nbdcnd, cofy
c                              are changed from a previous call
c                          = 1 if c, d, n, nbdcnd, cofy are unchanged
c                              from the previous call.
c
c                        iorder
c                          = 2 if a second-order approximation
c                              is sought
c                          = 4 if a fourth-order approximation
c                              is sought
c
c                        a,b
c                          the range of the x-independent variable,
c                          i.e., x is greater than or equal to a
c                          and less than or equal to b.  a must be
c                          less than b.
c
c                        m
c                          the number of panels into which the
c                          interval [a,b] is subdivided. hence,
c                          there will be m+1 grid points in the x-
c                          direction given by xi=a+(i-1)*dlx
c                          for i=1,2,...,m+1 where dlx=(b-a)/m is
c                          the panel width.  m must be less than
c                          idmn and greater than 5.
c
c                        mbdcnd
c                          indicates the type of boundary condition
c                          at x=a and x=b
c
c                          = 0 if the solution is periodic in x, i.e.,
c                              u(x+b-a,y)=u(x,y)  for all y,x
c                          = 1 if the solution is specified at x=a
c                              and x=b, i.e., u(a,y) and u(b,y) are
c                              specified for all y
c                          = 2 if the solution is specified at x=a and
c                              the boundary condition is mixed at x=b,
c                              i.e., u(a,y) and du(b,y)/dx+beta*u(b,y)
c                              are specified for all y
c                          = 3 if the boundary conditions at x=a and
c                              x=b are mixed, i.e.,
c                              du(a,y)/dx+alpha*u(a,y) and
c                              du(b,y)/dx+beta*u(b,y) are specified
c                              for all y
c                          = 4 if the boundary condition at x=a is
c                              mixed and the solution is specified
c                              at x=b, i.e., du(a,y)/dx+alpha*u(a,y)
c                              and u(b,y) are specified for all y
c
c                        bda
c                          a one-dimensional array of length n+1
c                          that specifies the values of
c                          du(a,y)/dx+ alpha*u(a,y) at x=a, when
c                          mbdcnd=3 or 4.
c                          bda(j) = du(a,yj)/dx+alpha*u(a,yj),
c                          j=1,2,...,n+1. when mbdcnd has any other
c                          other value, bda is a dummy parameter.
c
c                        alpha
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at x=a
c                          (see argument bda).  if mbdcnd is not
c                          equal to 3 or 4 then alpha is a dummy
c                          parameter.
c
c                        bdb
c                          a one-dimensional array of length n+1
c                          that specifies the values of
c                          du(b,y)/dx+ beta*u(b,y) at x=b.
c                          when mbdcnd=2 or 3
c                          bdb(j) = du(b,yj)/dx+beta*u(b,yj),
c                          j=1,2,...,n+1. when mbdcnd has any other
c                          other value, bdb is a dummy parameter.
c
c                        beta
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          x=b (see argument bdb).  if mbdcnd is
c                          not equal to 2 or 3 then beta is a dummy
c                          parameter.
c
c                        c,d
c                          the range of the y-independent variable,
c                          i.e., y is greater than or equal to c
c                          and less than or equal to d.  c must be
c                          less than d.
c
c                        n
c                          the number of panels into which the
c                          interval [c,d] is subdivided.
c                          hence, there will be n+1 grid points
c                          in the y-direction given by
c                          yj=c+(j-1)*dly for j=1,2,...,n+1 where
c                          dly=(d-c)/n is the panel width.
c                          in addition, n must be greater than 4.
c
c                        nbdcnd
c                          indicates the types of boundary conditions
c                          at y=c and y=d
c
c                          = 0 if the solution is periodic in y,
c                              i.e., u(x,y+d-c)=u(x,y)  for all x,y
c                          = 1 if the solution is specified at y=c
c                              and y = d, i.e., u(x,c) and u(x,d)
c                              are specified for all x
c                          = 2 if the solution is specified at y=c
c                              and the boundary condition is mixed
c                              at y=d, i.e., u(x,c) and
c                              du(x,d)/dy+xnu*u(x,d) are specified
c                              for all x
c                          = 3 if the boundary conditions are mixed
c                              at y=c and y=d, i.e.,
c                              du(x,d)/dy+gama*u(x,c) and
c                              du(x,d)/dy+xnu*u(x,d) are specified
c                              for all x
c                          = 4 if the boundary condition is mixed
c                              at y=c and the solution is specified
c                              at y=d, i.e. du(x,c)/dy+gama*u(x,c)
c                              and u(x,d) are specified for all x
c
c                        bdc
c                          a one-dimensional array of length m+1
c                          that specifies the value of
c                          du(x,c)/dy+gama*u(x,c) at y=c.
c                          when nbdcnd=3 or 4 bdc(i) = du(xi,c)/dy +
c                          gama*u(xi,c), i=1,2,...,m+1.
c                          when nbdcnd has any other value, bdc
c                          is a dummy parameter.
c
c                        gama
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          y=c (see argument bdc).  if nbdcnd is
c                          not equal to 3 or 4 then gama is a dummy
c                          parameter.
c
c                        bdd
c                          a one-dimensional array of length m+1
c                          that specifies the value of
c                          du(x,d)/dy + xnu*u(x,d) at y=c.
c                          when nbdcnd=2 or 3 bdd(i) = du(xi,d)/dy +
c                          xnu*u(xi,d), i=1,2,...,m+1.
c                          when nbdcnd has any other value, bdd
c                          is a dummy parameter.
c
c                        xnu
c                          the scalar multiplying the solution in
c                          case of a mixed boundary condition at
c                          y=d (see argument bdd).  if nbdcnd is
c                          not equal to 2 or 3 then xnu is a
c                          dummy parameter.
c
c                        cofx
c                          a user-supplied subprogram with
c                          parameters x, afun, bfun, cfun which
c                          returns the values of the x-dependent
c                          coefficients af(x), bf(x), cf(x) in the
c                          elliptic equation at x.
c
c                        cofy
c                          a user-supplied subprogram with parameters
c                          y, dfun, efun, ffun which returns the
c                          values of the y-dependent coefficients
c                          df(y), ef(y), ff(y) in the elliptic
c                          equation at y.
c
c                          note:  cofx and cofy must be declared
c                          external in the calling routine.
c                          the values returned in afun and dfun
c                          must satisfy afun*dfun greater than 0
c                          for a less than x less than b, c less
c                          than y less than d (see ierror=10).
c                          the coefficients provided may lead to a
c                          matrix equation which is not diagonally
c                          dominant in which case solution may fail
c                          (see ierror=4).
c
c                        grhs
c                          a two-dimensional array that specifies the
c                          values of the right-hand side of the
c                          elliptic equation, i.e.,
c                          grhs(i,j)=g(xi,yi), for i=2,...,m,
c                          j=2,...,n.  at the boundaries, grhs is
c                          defined by
c
c                          mbdcnd   grhs(1,j)   grhs(m+1,j)
c                          ------   ---------   -----------
c                            0      g(a,yj)     g(b,yj)
c                            1         *           *
c                            2         *        g(b,yj)  j=1,2,...,n+1
c                            3      g(a,yj)     g(b,yj)
c                            4      g(a,yj)        *
c
c                          nbdcnd   grhs(i,1)   grhs(i,n+1)
c                          ------   ---------   -----------
c                            0      g(xi,c)     g(xi,d)
c                            1         *           *
c                            2         *        g(xi,d)  i=1,2,...,m+1
c                            3      g(xi,c)     g(xi,d)
c                            4      g(xi,c)        *
c
c                          where * means these quantities are not used.
c                          grhs should be dimensioned idmn by at least
c                          n+1 in the calling routine.
c
c                        usol
c                          a two-dimensional array that specifies the
c                          values of the solution along the boundaries.
c                          at the boundaries, usol is defined by
c
c                          mbdcnd   usol(1,j)   usol(m+1,j)
c                          ------   ---------   -----------
c                            0         *           *
c                            1      u(a,yj)     u(b,yj)
c                            2      u(a,yj)        *     j=1,2,...,n+1
c                            3         *           *
c                            4         *        u(b,yj)
c
c                          nbdcnd   usol(i,1)   usol(i,n+1)
c                          ------   ---------   -----------
c                            0         *           *
c                            1      u(xi,c)     u(xi,d)
c                            2      u(xi,c)        *     i=1,2,...,m+1
c                            3         *           *
c                            4         *        u(xi,d)
c
c                          where * means the quantities are not used
c                          in the solution.
c
c                          if iorder=2, the user may equivalence grhs
c                          and usol to save space.  note that in this
c                          case the tables specifying the boundaries
c                          of the grhs and usol arrays determine the
c                          boundaries uniquely except at the corners.
c                          if the tables call for both g(x,y) and
c                          u(x,y) at a corner then the solution must
c                          be chosen.  for example, if mbdcnd=2 and
c                          nbdcnd=4, then u(a,c), u(a,d), u(b,d) must
c                          be chosen at the corners in addition
c                          to g(b,c).
c
c                          if iorder=4, then the two arrays, usol and
c                          grhs, must be distinct.
c
c                          usol should be dimensioned idmn by at least
c                          n+1 in the calling routine.
c
c                        idmn
c                          the row (or first) dimension of the arrays
c                          grhs and usol as it appears in the program
c                          calling sepeli.  this parameter is used
c                          to specify the variable dimension of grhs
c                          and usol.  idmn must be at least 7 and
c                          greater than or equal to m+1.
c
c                        w
c                          a one-dimensional array that must be
c                          provided by the user for work space.
c                          let k=int(log2(n+1))+1 and set l=2**(k+1).
c                          then (k-2)*l+k+10*n+12*m+27 will suffice
c                          as a length of w.  the actual length of w
c                          in the calling routine must be set in w(1)
c                          (see ierror=11).
c
c on output              usol
c                          contains the approximate solution to the
c                          elliptic equation.
c                          usol(i,j) is the approximation to u(xi,yj)
c                          for i=1,2...,m+1 and j=1,2,...,n+1.
c                          the approximation has error
c                          o(dlx**2+dly**2) if called with iorder=2
c                          and o(dlx**4+dly**4) if called with
c                          iorder=4.
c
c                        w
c                          contains intermediate values that must not
c                          be destroyed if sepeli is called again with
c                          intl=1.  in addition w(1) contains the
c                          exact minimal length (in floating point)
c                          required for the work space (see ierror=11).
c
c                        pertrb
c                          if a combination of periodic or derivative
c                          boundary conditions
c                          (i.e., alpha=beta=0 if mbdcnd=3;
c                          gama=xnu=0 if nbdcnd=3) is specified
c                          and if the coefficients of u(x,y) in the
c                          separable elliptic equation are zero
c                          (i.e., cf(x)=0 for x greater than or equal
c                          to a and less than or equal to b;
c                          ff(y)=0 for y greater than or equal to c
c                          and less than or equal to d) then a
c                          solution may not exist.  pertrb is a
c                          constant calculated and subtracted from
c                          the right-hand side of the matrix equations
c                          generated by sepeli which insures that a
c                          solution exists. sepeli then computes this
c                          solution which is a weighted minimal least
c                          squares solution to the original problem.
c
c                        ierror
c                          an error flag that indicates invalid input
c                          parameters or failure to find a solution
c                          = 0 no error
c                          = 1 if a greater than b or c greater than d
c                          = 2 if mbdcnd less than 0 or mbdcnd greater
c                              than 4
c                          = 3 if nbdcnd less than 0 or nbdcnd greater
c                              than 4
c                          = 4 if attempt to find a solution fails.
c                              (the linear system generated is not
c                              diagonally dominant.)
c                          = 5 if idmn is too small
c                              (see discussion of idmn)
c                          = 6 if m is too small or too large
c                              (see discussion of m)
c                          = 7 if n is too small (see discussion of n)
c                          = 8 if iorder is not 2 or 4
c                          = 9 if intl is not 0 or 1
c                          = 10 if afun*dfun less than or equal to 0
c                               for some interior mesh point (xi,yj)
c                          = 11 if the work space length input in w(1)
c                               is less than the exact minimal work
c                               space length required output in w(1).
c
c                          note (concerning ierror=4):  for the
c                          coefficients input through cofx, cofy,
c                          the discretization may lead to a block
c                          tridiagonal linear system which is not
c                          diagonally dominant (for example, this
c                          happens if cfun=0 and bfun/(2.*dlx) greater
c                          than afun/dlx**2).  in this case solution
c                          may fail.  this cannot happen in the limit
c                          as dlx, dly approach zero.  hence, the
c                          condition may be remedied by taking larger
c                          values for m or n.
c
c special conditions     see cofx, cofy argument descriptions above.
c
c i/o                    none
c
c precision              single
c
c required library       blktri, comf, q8qst4, and sepaux, which are
c files                  loaded by default on ncar's cray machines.
c
c language               fortran
c
c history                developed at ncar during 1975-76 by
c                        john c. adams of the scientific computing
c                        division.  released on ncar's public software
c                        libraries in january 1980.
c
c portability            fortran 66
c
c algorithm              sepeli automatically discretizes the
c                        separable elliptic equation which is then
c                        solved by a generalized cyclic reduction
c                        algorithm in the subroutine, blktri.  the
c                        fourth-order solution is obtained using
c                        'deferred corrections' which is described
c                        and referenced in sections, references and
c                        method.
c
c timing                 the operational count is proportional to
c                        m*n*log2(n).
c
c accuracy               the following accuracy results were obtained
c                        on a cdc 7600.  note that the fourth-order
c                        accuracy is not realized until the mesh is
c                        sufficiently refined.
c
c                                     second-order  fourth-order
c                            m    n     error         error
c
c                             6    6    6.8e-1        1.2e0
c                            14   14    1.4e-1        1.8e-1
c                            30   30    3.2e-2        9.7e-3
c                            62   62    7.5e-3        3.0e-4
c                           126  126    1.8e-3        3.5e-6
c
c portability            fortran 66
c
c references             keller, h.b., numerical methods for two-point
c                        boundary-value problems, blaisdel (1968),
c                        waltham, mass.
c
c                        swarztrauber, p., and r. sweet (1975):
c                        efficient fortran subprograms for the
c                        solution of elliptic partial differential
c                        equations.  ncar technical note
c                        ncar-tn/ia-109, pp. 135-137.
c***********************************************************************



      implicit none
      INT_TYPE INTL,IORDER, M,MBDCND, N,NBDCND, IDMN, IERROR
      REAL_TYPE  A,B, BDA(N+1),ALPHA, BDB(N+1),BETA, C,D,
     &           BDC(M+1),GAMA,  BDD(M+1),XNU,      
     &       GRHS(IDMN,N+1), USOL(IDMN,N+1), W(*), PERTRB

!-- internal -->

      INT_TYPE L, LOGB2N, LL, K, LENGTH, LINPUT, LOUTPT,
     &           I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13
#ifdef VERBOSE
      write(*,*) 'SEPELI: INTL,IORDER =', INTL,IORDER,
     &           '  A,B =', A,B,  ' M,MBDCND =',M,MBDCND,
     &           '  C,D =', C,D,  ' N,NBDCND=', N,NBDCND
#endif



!
! Check input parameters
!
      call CHKPRM(INTL,IORDER, A,B,M, MBDCND, C,D,N, NBDCND,
     &                                         IDMN, IERROR)
      if (IERROR.NE.0) return
!
! Compute minimum work space and check work space length input
!
      L = N + 1
      If (NBDCND.EQ.0) L = N
      LOGB2N = INT(ALOG(FLOAT(L)+0.5)/ALOG(2.0)) + 1
      LL = 2 ** (LOGB2N+1)
      K = M + 1
      L = N + 1
      LENGTH = (LOGB2N-2) * LL + LOGB2N + MAX0(2*L,6*K) + 5
      If (NBDCND.EQ.0) LENGTH = LENGTH + 2 * L
      IERROR = 11
      LINPUT = INT(W(1)+0.5)
      LOUTPT = LENGTH + 6 * (K+L) + 1
      W(1) = FLOAT(LOUTPT)
      If (LOUTPT.GT.LINPUT) return
      IERROR = 0
!
! Set work space indices
!
      I1 = LENGTH + 2
      I2 = I1 + L
      I3 = I2 + L
      I4 = I3 + L
      I5 = I4 + L
      I6 = I5 + L
      I7 = I6 + L
      I8 = I7 + K
      I9 = I8 + K
      I10 = I9 + K
      I11 = I10 + K
      I12 = I11 + K
      I13 = 2

      call SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N,
     &            NBDCND,BDC,GAMA,BDD,XNU,W(I1),W(I2),W(I3),W(I4),
     &            W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),W(I12),
     &            GRHS,USOL, IDMN,W(I13),PERTRB,IERROR)
      return
      end




      subroutine SPELIP (INTL,IORDER, A,B, M,MBDCND, BDA,ALPHA,
     &                   BDB,BETA, C,D, N,NBDCND, BDC,GAMA, BDD,XNU,
     &                   AN,BN,CN,DN,UN,ZN, AM,BM,CM,DM,UM,ZM,
     &                   GRHS,USOL, IDMN, W, PERTRB,IERROR)
!
! spelip sets up vectors and arrays for input to blktri and computes
! a second-order solution in usol.  a return jump to sepeli occurrs
! if iorder=2.  if iorder=4 a fourth order solution is generated in
! usol.
!

      implicit none
      INT_TYPE INTL,IORDER, M,MBDCND, N,NBDCND, IDMN, IERROR
      REAL_TYPE A,B, BDA(1),ALPHA, BDB(1),BETA,
     &          C,D, BDC(1),GAMA,  BDD(1),XNU,
     &          AN(1),BN(1),CN(1),DN(1),UN(1),ZN(1),
     &          AM(1),BM(1),CM(1),DM(1),UM(1),ZM(1),
     &          GRHS(IDMN,1), USOL(IDMN,1), W(1), PERTRB

!-- internal -->

      logical SINGLR
      INT_TYPE I,J, MP,NP, I1, IORD
      REAL_TYPE XI, AI,BI,CI, AXI,BXI,CXI, AX1,CXM,
     &          YJ, DJ,EJ,FJ, DYJ,EYJ,FYJ, DY1,FYN, PRTRB

#include "splp.h"

!
! set parameters internally
!
      KSWX = MBDCND + 1
      KSWY = NBDCND + 1
      K = M + 1
      L = N + 1
      AIT = A
      BIT = B
      CIT = C
      DIT = D
!
! Set right hand side values from grhs in usol on the interior
! and non-specified boundaries.
!
      do I=2,M
        do J=2,N
          USOL(I,J) = GRHS(I,J)
        enddo
      enddo
      if (KSWX.NE.2 .AND. KSWX.NE.3) then
        do J=2,N
          USOL(1,J) = GRHS(1,J)
        enddo
      endif
      if (KSWX.NE.2 .AND. KSWX.NE.5) then
        do J=2,N
          USOL(K,J) = GRHS(K,J)
        enddo
      endif
      if (KSWY.NE.2 .AND. KSWY.NE.3) then
        do I=2,M
          USOL(I,1) = GRHS(I,1)
        enddo
      endif
      if (KSWY.NE.2 .AND. KSWY.NE.5) then
        do I=2,M
          USOL(I,L) = GRHS(I,L)
        enddo
      endif
      If (KSWX.NE.2 .AND. KSWX.NE.3 .AND. KSWY.NE.2 .AND. KSWY.NE.3) 
     &    USOL(1,1) = GRHS(1,1)
      If (KSWX.NE.2 .AND. KSWX.NE.5 .AND. KSWY.NE.2 .AND. KSWY.NE.3) 
     &    USOL(K,1) = GRHS(K,1)
      If (KSWX.NE.2 .AND. KSWX.NE.3 .AND. KSWY.NE.2 .AND. KSWY.NE.5) 
     &    USOL(1,L) = GRHS(1,L)
      If (KSWX.NE.2 .AND. KSWX.NE.5 .AND. KSWY.NE.2 .AND. KSWY.NE.5) 
     &    USOL(K,L) = GRHS(K,L)
      I1 = 1
!
! Set switches for periodic or non-periodic boundaries
!
      MP = 1
      NP = 1
      If (KSWX.EQ.1) MP = 0
      If (KSWY.EQ.1) NP = 0
!
! Set dlx,dly and size of block tri-diagonal system generated
! in nint,mint
!
      DLX = (BIT-AIT) / FLOAT(M)
      MIT = K - 1
      If (KSWX.EQ.2) MIT = K - 2
      If (KSWX.EQ.4) MIT = K
      DLY = (DIT-CIT) / FLOAT(N)
      NIT = L - 1
      If (KSWY.EQ.2) NIT = L - 2
      If (KSWY.EQ.4) NIT = L
      TDLX3 = 2.0 * DLX ** 3
      DLX4 = DLX ** 4
      TDLY3 = 2.0 * DLY ** 3
      DLY4 = DLY ** 4
!
! Set subscript limits for portion of array to input to blktri
!
      IS = 1
      JS = 1
      If (KSWX.EQ.2.OR.KSWX.EQ.3) IS = 2
      If (KSWY.EQ.2.OR.KSWY.EQ.3) JS = 2
      NS = NIT + JS - 1
      MS = MIT + IS - 1
!
! Set x - direction
!
      do I=1,MIT
        XI = AIT + FLOAT(IS+I-2) * DLX
        call COFX(XI,AI,BI,CI)
        AXI = (AI/DLX-0.5*BI) / DLX
        BXI = -2. * AI / DLX ** 2 + CI
        CXI = (AI/DLX+0.5*BI) / DLX
        AM(I) = AXI
        BM(I) = BXI
        CM(I) = CXI
      enddo
!
! Set y direction
!
      do J=1,NIT
        YJ = CIT + FLOAT(JS+J-2) * DLY
        call COFY(YJ,DJ,EJ,FJ)
        DYJ = (DJ/DLY-0.5*EJ) / DLY
        EYJ = (-2.*DJ/DLY**2+FJ)
        FYJ = (DJ/DLY+0.5*EJ) / DLY
        AN(J) = DYJ
        BN(J) = EYJ
        CN(J) = FYJ
      enddo
!
!  adjust edges in x direction unless periodic
!
      AX1 = AM(1)
      CXM = CM(MIT)
      Go To (130,90,110,120,100), KSWX
c
c     dirichlet-dirichlet in x direction
c
   90 AM(1) = 0.0
      CM(MIT) = 0.0
      Go To 130
c
c     mixed-dirichlet in x direction
c
  100 AM(1) = 0.0
      BM(1) = BM(1) + 2. * ALPHA * DLX * AX1
      CM(1) = CM(1) + AX1
      CM(MIT) = 0.0
      Go To 130
c
c     dirichlet-mixed in x direction
c
  110 AM(1) = 0.0
      AM(MIT) = AM(MIT) + CXM
      BM(MIT) = BM(MIT) - 2. * BETA * DLX * CXM
      CM(MIT) = 0.0
      Go To 130
c
c     mixed - mixed in x direction
c
  120 Continue
      AM(1) = 0.0
      BM(1) = BM(1) + 2. * DLX * ALPHA * AX1
      CM(1) = CM(1) + AX1
      AM(MIT) = AM(MIT) + CXM
      BM(MIT) = BM(MIT) - 2. * DLX * BETA * CXM
      CM(MIT) = 0.0
  130 Continue
c
c     adjust in y direction unless periodic
c
      DY1 = AN(1)
      FYN = CN(NIT)
      Go To (180,140,160,170,150), KSWY
c
c     dirichlet-dirichlet in y direction
c
  140 Continue
      AN(1) = 0.0
      CN(NIT) = 0.0
      Go To 180
c
c     mixed-dirichlet in y direction
c
  150 Continue
      AN(1) = 0.0
      BN(1) = BN(1) + 2. * DLY * GAMA * DY1
      CN(1) = CN(1) + DY1
      CN(NIT) = 0.0
      Go To 180
c
c     dirichlet-mixed in y direction
c
  160 AN(1) = 0.0
      AN(NIT) = AN(NIT) + FYN
      BN(NIT) = BN(NIT) - 2. * DLY * XNU * FYN
      CN(NIT) = 0.0
      Go To 180
c
c     mixed - mixed direction in y direction
c
  170 Continue
      AN(1) = 0.0
      BN(1) = BN(1) + 2. * DLY * GAMA * DY1
      CN(1) = CN(1) + DY1
      AN(NIT) = AN(NIT) + FYN
      BN(NIT) = BN(NIT) - 2.0 * DLY * XNU * FYN
      CN(NIT) = 0.0
  180 continue
!
! adjust usol along x edge
!
      if (KSWX.NE.1) then
        do J=JS,NS
          if (KSWX.EQ.2 .OR. KSWX.EQ.3) then
            USOL(IS,J) = USOL(IS,J) - AX1 * USOL(1,J)
          else
            USOL(IS,J) = USOL(IS,J) + 2.0 * DLX * AX1 * BDA(J)
          endif
          if (KSWX.EQ.2 .OR. KSWX.EQ.5) then
            USOL(MS,J) = USOL(MS,J) - CXM * USOL(K,J)
          else
            USOL(MS,J) = USOL(MS,J) - 2.0 * DLX * CXM * BDB(J)
          endif
        enddo
      endif
!
! adjust usol along y edge
!
      if (KSWY.NE.1) then
        do I=IS,MS
          if (KSWY.EQ.2 .OR. KSWY.EQ.3) then
            USOL(I,JS) = USOL(I,JS) - DY1 * USOL(I,1)
          else
            USOL(I,JS) = USOL(I,JS) + 2.0 * DLY * DY1 * BDC(I)
          endif 
          if (KSWY.EQ.2 .OR. KSWY.EQ.5) then
            USOL(I,NS) = USOL(I,NS) - FYN * USOL(I,L)
          else
            USOL(I,NS) = USOL(I,NS) - 2.0 * DLY * FYN * BDD(I)
          endif
        enddo
      endif
!
! save adjusted edges in grhs if iorder=4
!
      if (IORDER.EQ.4) then
        do J = JS, NS
          GRHS(IS,J) = USOL(IS,J)
          GRHS(MS,J) = USOL(MS,J)
        enddo
        do I = IS, MS
          GRHS(I,JS) = USOL(I,JS)
          GRHS(I,NS) = USOL(I,NS)
        enddo
      endif
      IORD = IORDER
      PERTRB = 0.0
!
! check if operator is singular
!
      call CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,SINGLR)
!
! compute non-zero eigenvector in null space of transpose
! if singular
!
      if (SINGLR) call SEPTRI(MIT,AM,BM,CM,DM,UM,ZM)
      if (SINGLR) call SEPTRI(NIT,AN,BN,CN,DN,UN,ZN)
!
! make initialization call to blktri
!
      if (INTL.EQ.0) call BLKTRI(INTL,NP,NIT, AN,BN,CN, MP,MIT,
     &                    AM,BM,CM, IDMN,USOL(IS,JS), IERROR,W)
      if (IERROR.NE.0) return
!
! adjust right hand side if necessary
!
  230 Continue
        if (SINGLR) call SEPORT(USOL,IDMN,ZN,ZM,PERTRB)
!
! compute solution
!
        call BLKTRI (I1,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,IDMN,
     &                      USOL(IS,JS), IERROR,W)
        if (IERROR.NE.0) return
!
! set periodic boundaries if necessary
!
        if (KSWX.EQ.1) then
          do J=1,L
            USOL(K,J) = USOL(1,J)
          enddo
        endif
        if (KSWY.EQ.1) then
          do I=1,K
            USOL(I,L) = USOL(I,1)
          enddo
        endif
!
! minimize solution with respect to weighted least squares
! norm if operator is singular
!
        if (SINGLR) call SEPMIN(USOL,IDMN,ZN,ZM,PRTRB)
!
! return if deferred corrections and a fourth order solution are
! not flagged
!
        if (IORD.EQ.2) return
          IORD = 2
!
! compute new right hand side for fourth order solution
!
        call DEFER(IDMN,USOL,GRHS)
       goto 230
      end





      subroutine CHKPRM (INTL,IORDER, A,B, M,MBDCND,
     &                   C,D, N,NBDCND, IDMN,IERROR)
!
! this program checks the input parameters for errors
!
      implicit none
      INT_TYPE INTL,IORDER, M,MBDCND, N,NBDCND, IDMN,IERROR
      REAL_TYPE A,B, C,D
!-- internal -->
      INT_TYPE I,J
      REAL_TYPE DLX, XI, AI,BI,CI,   DLY, YJ, DJ,EJ,FJ

      IERROR=0
      if (A.GE.B .or. C.GE.D) then             ! check definition
         IERROR=1                              ! of solution region

      elseif (MBDCND.LT.0 .or. MBDCND.GT.4) then   ! check boundary
        IERROR=2                                   ! switches
      elseif (NBDCND.LT.0 .or. NBDCND.GT.4) then
        IERROR=3

      elseif (IDMN.LT.7) then               ! check first dimension
         IERROR=5                           ! in calling routine

      elseif (M.GT.(IDMN-1) .or. M.LT.6) then      ! check m
        IERROR=6

      elseif (N.LT.5) then                         ! check n
        IERROR=7

      elseif (IORDER.NE.2 .AND. IORDER.NE.4) then  ! check iorder
        IERROR=8

      elseif (INTL.NE.0 .AND. INTL.NE.1) then      ! check intl
        IERROR=9
      else
        DLX=(B-A)/FLOAT(M)                         ! check whether 
        DLY=(D-C)/FLOAT(N)                         ! the equation is
        do I = 2, M                                ! elliptic
          XI=A+FLOAT(I-1)*DLX
          call COFX(XI,AI,BI,CI)
          do J=2,N
            YJ = C + FLOAT(J-1) * DLY
            call COFY(YJ,DJ,EJ,FJ)
            if (AI*DJ.LE.0.0) then
              IERROR = 10
            endif 
          enddo
        enddo
      endif
      return
      end




      subroutine CHKSNG (MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,SINGLR)
!
! check if the pde sepeli must solve is a singular operator
!
      implicit none
      INT_TYPE MBDCND,NBDCND 
      REAL_TYPE ALPHA,BETA,GAMA,XNU
      logical SINGLR
#include "splp.h"
      INT_TYPE I,J
      REAL_TYPE XI,AI,BI,CI, YJ,DJ,EJ,FJ

      SINGLR = .FALSE.
!
!  check if the boundary conditions are entirely periodic and/or mixed
!
      if ((MBDCND.NE.0 .AND. MBDCND.NE.3) .OR.
     &    (NBDCND.NE.0 .AND. NBDCND.NE.3)) return
!
! check that mixed conditions are pure neuman
!
      if (MBDCND.EQ.3) then 
        if (ALPHA.NE.0.0 .or. BETA.NE.0.0) return
      endif 
      if (NBDCND.EQ.3) then
        If (GAMA.NE.0.0  .or. XNU.NE.0.0) return
      endif 

      do I = IS, MS                    ! Check that non-derivative
        XI = AIT + FLOAT(I-1) * DLX    ! coefficient functions are
        call COFX(XI,AI,BI,CI)         ! zero 
        if (CI.NE.0.0) return
      enddo
      do J = JS, NS
        YJ = CIT + FLOAT(J-1) * DLY
        call COFY(YJ,DJ,EJ,FJ)
        if (FJ.NE.0.0) return
      enddo
!
! the operator must be singular if this point is reached
!
      SINGLR = .TRUE.
      return
      end





      subroutine DEFER(IDMN,USOL,GRHS)
!
! this subroutine first approximates the truncation error given by
!
!              trun1(x,y)=dlx**2*tx+dly**2*ty
! where
!           tx=afun(x)*uxxxx/12.0+bfun(x)*uxxx/6.0
!
! on the interior and at the boundaries if periodic (here uxxx,uxxxx
! are the third and fourth partial derivatives of u with respect to
! x).  tx is of the form afun(x)/3.0*(uxxxx/4.0+uxxx/dlx) at x=a or
! x=b if the boundary condition there is mixed.  tx=0.0 along
! specified boundaries.  ty has symmetric form in y with x,afun(x),
! bfun(x) replaced by y,dfun(y),efun(y).  the second order solution
! in usol is used to approximate (via second order finite
! differencing) the truncation error and the result is added to the
! right hand side in grhs and then transferred to usol to be used as
! a new right hand side when calling blktri for a fourth order
! solution.
!

      implicit none
      INT_TYPE IDMN
      REAL_TYPE GRHS(IDMN,1), USOL(IDMN,1)
#include "splp.h"
      INT_TYPE I,J
      REAL_TYPE  XI,AI,BI,CI, TX, UXXX,UXXXX,
     &           YJ,DJ,EJ,FJ, TY, UYYY,UYYYY

!
! compute truncation error approximation over the entire mesh
!
      write(*,*) 'DEFER'

      do J = JS, NS
        YJ = CIT + FLOAT(J-1) * DLY
        call COFY(YJ,DJ,EJ,FJ)
        do I = IS, MS
          XI = AIT + FLOAT(I-1) * DLX
          call COFX(XI,AI,BI,CI)
!
! compute partial derivative approximations at (xi,yj)
!
          call SEPDX(USOL,IDMN,I,J,UXXX,UXXXX)
          call SEPDY(USOL,IDMN,I,J,UYYY,UYYYY)
          TX = AI * UXXXX / 12.0 + BI * UXXX / 6.0
          TY = DJ * UYYYY / 12.0 + EJ * UYYY / 6.0
!
! reset form of truncation if at boundary which is non-periodic
!
          if (.NOT.(KSWX.EQ.1.OR.(I.GT.1.AND.I.LT.K))) then
            TX = AI / 3.0 * (UXXXX/4.0+UXXX/DLX)
          endif 
          if (.NOT.(KSWY.EQ.1.OR.(J.GT.1.AND.J.LT.L))) then
            TY = DJ / 3.0 * (UYYYY/4.0+UYYY/DLY)
          endif 
          GRHS(I,J) = GRHS(I,J) + DLX ** 2 * TX + DLY ** 2 * TY
        enddo
      enddo
!
! reset the right hand side in usol
!
      do I = IS, MS
        do J = JS, NS
          USOL(I,J) = GRHS(I,J)
        enddo
      enddo
      return
      end








      subroutine BLKTRI (IFLG, NP,N, AN,BN,CN, MP,M, AM,BM,CM,
     &                      IDIMY,Y,IERROR, W)
!
! dimension of    an(n),bn(n),cn(n),am(m),bm(m),cm(m),y(idimy,n),
! arguments          w(see argument list)
!
! latest revision        january 1985
!
! purpose      blktri solves a system of linear equations of the form
!
!              an(j)*x(i,j-1) + am(i)*x(i-1,j) + (bn(j)+bm(i))*x(i,j)
!                    + cn(j)*x(i,j+1) + cm(i)*x(i+1,j) = y(i,j)
!
!              for i = 1,2,...,m  and  j = 1,2,...,n.
!
!                  i+1 and i-1 are evaluated modulo m
!              and j+1 and j-1 modulo n, i.e.,
!
!                     x(i,0) = x(i,n),  x(i,n+1) = x(i,1),
!                     x(0,j) = x(m,j),  x(m+1,j) = x(1,j).
!
!              these equations usually result from the discretization
!              of separable elliptic equations.  boundary conditions
!              may be Dirichlet, Neumann, or Periodic.
!
! arguments
!
! on input     iflg = 0  initialization only:  certain quantities
!                        that depend on np, n, an, bn, and cn are
!                        computed and stored in the work array w.
!
!                   = 1  the quantities that were computed in the
!                        initialization are used to obtain the
!                        solution x(i,j).
!
!                  note: a call with iflg=0 takes approximately one
!                        half the time as a call with iflg = 1.
!                        however, the initialization does not have
!                        to be repeated unless np, n, an, bn, or cn
!                        change.
!
!                np = 0  if an(1) and cn(n) are not zero, which
!                        corresponds to periodic bounary conditions.
!
!                   = 1  if an(1) and cn(n) are zero.
!
!              n   the number of unknowns in the j-direction.
!                  n must be greater than 4.  the operation count
!                  is proportional to mnlog2(n), hence n should
!                  be selected less than or equal to m.
!
!              an,bn,cn  one-dimensional arrays of length n that
!                        specify the coefficients in the linear
!                        equations given above.
!
!                mp = 0  if am(1) and cm(m) are not zero, which
!                        corresponds to periodic boundary conditions.
!
!                   = 1  if am(1) = cm(m) = 0  .
!
!              m   the number of unknowns in the i-direction.
!                  m must be greater than 4.
!
!              am,bm,cm  one-dimensional arrays of length m that
!                        specify the coefficients in the linear
!                        equations given above.
!
!              idimy   the row (or first) dimension of the two-
!                      dimensional array y as it appears in the
!                      program calling blktri.  this parameter is
!                      used to specify the variable dimension of y.
!                      idimy must be at least m.
!
!              y  a two-dimensional array that specifies the values
!                 of the right side of the linear system of equations
!                 given above.  y must be dimensioned at least m*n.
!
!              w  a one-dimensional array that must be provided by
!                 the user for work space.  if np=1 define
!                 k=int(log2(n))+1 and set l=2**(k+1) then w must
!                 have dimension (k-2)*l+k+5+max(2n,6m)
!
!                 if np=0 define k=int(log2(n-1))+1 and set
!                 l=2**(k+1) then w must have dimension
!                          (k-2)*l+k+5+2n+max(2n,6m)
!
!  **important**  for purposes of checking, the required dimension
!                 of w is computed by blktri and stored in w(1) in
!                 floating point format.
!
!
! on output    y   contains the solution x.
!
!              ierror   an error flag that indicates invalid
!                       input parameters.  except for number zero,
!                       a solution is not attempted.
!
!                    = 0  no error.
!                    = 1  m is less than 5
!                    = 2  n is less than 5
!                    = 3  idimy is less than m.
!                    = 4  blktri failed while computing results
!                         that depend on the coefficient arrays
!                         an, bn, cn.  check these arrays.
!                    = 5  an(j)*cn(j-1) is less than 0 for some j.
!
!                  possible reasons for this condition are
!                    1. the arrays an and cn are not correct
!                    2. too large a grid spacing was used in the
!                          discretization of the elliptic equation.
!                    3. the linear equations resulted from a partial
!                       differential equation which was not elliptic.
!
!              w   contains intermediate values that must not be
!                  destroyed if blktri will be called again with
!                  iflg=1. w(1) contains the number of locations
!                  required by w in floating point format.
!
!
! special      the algorithm may fail if abs(bm(i)+bn(j)) is less
! conditions   than abs(am(i))+abs(an(j)) + abs(cm(i))+abs(cn(j))
!              for some i and j. the algorithm will also fail if
!              an(j)*cn(j-1) is less than zero for some j.
!              see the description of the output parameter ierror.
!
!
! required         comf and q8qst4, which are automatically
! library files    loaded by ncar's cray machines.
!
! history      written by Paul Swarztrauber at ncar in the early
!              1970's.  rewritten and released in january, 1980.
!
! algorithm    generalized cyclic reduction
!
!
! references   Swarztrauber, P. and R. Sweet, 'Efficient fortran
!              subprograms for the solution of elliptic equations'
!              NCAR tn/ia-109, july, 1975, 138 pp.
!
!              Swarztrauber P. N., A direct method for the discrete
!              solution of separable elliptic equations, S.I.A.M.
!              j. numer. anal.,11(1974) pp. 1136-1150.
!
      implicit none
      INT_TYPE IFLG, NP,N,  MP,M, IDIMY, IERROR
      REAL_TYPE AN(N),BN(N),CN(N), AM(M),BM(M),CM(M), Y(IDIMY,M),W(*)
#include "cblkt.h"
      INT_TYPE NH,NL,IWAH, IW1,IWBH,IW2,IW3,IWD,IWW,IWU
      external PROD, PRODP, CPROD, CPRODP
#ifdef VERBOSE
      write(*,*) 'BLKTRI'
#endif
!
! test m and n for the proper form
!
      NM = N
      IERROR = 0
      If (M.LT.5) then
        IERROR = 1
      else
        If (NM.LT.3) then
          IERROR = 2
        else
          If (IDIMY.LT.M) then
            IERROR = 3
          else
            NH = N
            NPP = NP
            if (NPP.NE.0) NH=NH+1 
            IK = 2
            K = 1
   10        IK = IK + IK
             K = K + 1
            If (NH.GT.IK) goto 10
            NL = IK
            IK = IK + IK
            NL = NL - 1
            IWAH = (K-2) * IK + K + 6
            If (NPP.NE.0) then
!
! divide w into working sub arrays
!
              IW1 = IWAH
              IWBH = IW1 + NM
              W(1) = FLOAT(IW1-1+MAX0(2*NM,6*M))
            else
              IWBH = IWAH + NM + NM
              IW1 = IWBH
              W(1) = FLOAT(IW1-1+MAX0(2*NM,6*M))
              NM = NM - 1
            endif 
!
! subroutine comp b computes the roots of the b polynomials
!
            If (IERROR.EQ.0) then
              IW2 = IW1 + M
              IW3 = IW2 + M
              IWD = IW3 + M
              IWW = IWD + M
              IWU = IWW + M
              If (IFLG.EQ.0) then
                call COMPB(NL,IERROR,AN,BN,CN,W(2),W(IWAH),W(IWBH))
              else
                If (MP.NE.0) then
!
! subroutine blktr1 solves the linear system
!
                  call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
     &                   W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),
     &                                                PROD, CPROD)
                else
                  call BLKTR1(NL,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,W(2),
     &                   W(IW1),W(IW2),W(IW3),W(IWD),W(IWW),W(IWU),
     &                                              PRODP, CPRODP)
                endif 
              endif 
            endif 
          endif
        endif
      endif
      return
      end



      subroutine BLKTR1(N, AN,BN,CN, M, AM,BM,CM, IDIMY,Y,
     &                  B,W1,W2,W3,WD,WW,WU, PRDCT,CPRDCT)
!
! blktr1 solves the linear system
!
! b  contains the roots of all the b polynomials w1,w2,w3,wd,ww,wu
! are all working arrays  prdct  is either prodp or prod depending on
! whether the boundary conditions in the m direction are periodic or
! not  cprdct is either cprodp or cprod which are the complex
! versions of prodp and prod. these are called in the event that some
! of the roots of the b sub p polynomial are complex
!
      implicit none
      INT_TYPE N, M, IDIMY
      REAL_TYPE AN(1), BN(1), CN(1), AM(1), BM(1), CM(1), B(1),
     &     W1(1), W2(1), W3(1), WD(1), WW(1), WU(1), Y(IDIMY,1)
      external PRDCT,CPRDCT

#include "cblkt.h"

      REAL_TYPE DUM
      INT_TYPE KDO,L, IF,IFD, I,I1,I2,I3,I4, IZ,IZR, IDXC,IDXA,
     &         IR,IRM1, IM1,IM2,IM3,IMI1,IMI2, NA,NC,NZ, NM1,NM2,NM3, 
     &         IP,IP1,IP2,IP3,IPI1,IPI2,IPI3, NP,NP1,NP2,NP3, J, LL,
     &         izero
      parameter (izero=0)

#ifdef VERBOSE
      write(*,*) 'BLKTR1'
#endif

      KDO = K - 1       !--> begin reduction phase
      do L = 1, KDO
        IR = L - 1
        I2 = 2 ** IR
        I1 = I2 / 2
        I3 = I2 + I1
        I4 = I2 + I2
        IRM1 = IR - 1
        call INDXB(I2,IR,IM2,NM2)
        call INDXB(I1,IRM1,IM3,NM3)
        call INDXB(I3,IRM1,IM1,NM1)
        call PRDCT(NM2,B(IM2), NM3,B(IM3), NM1,B(IM1),izero,DUM,
     &                    Y(1,I2),W3,M, AM,BM,CM,WD,WW,WU)
        IF = 2 ** K
        do I = I4, IF, I4
          If (I.LE.NM) then
            IPI1 = I + I1
            IPI2 = I + I2
            IPI3 = I + I3
            call INDXC(I,IR,IDXC,NC)
            If (I.LT.IF) then
              call INDXA(I,IR,IDXA,NA)
              call INDXB(I-I1,IRM1,IM1,NM1)
              call INDXB(IPI2,IR,IP2,NP2)
              call INDXB(IPI1,IRM1,IP1,NP1)
              call INDXB(IPI3,IRM1,IP3,NP3)
              call PRDCT(NM1,B(IM1), izero,DUM, izero,DUM, NA,
     &                      AN(IDXA),W3,W1,M,AM, BM,CM,WD,WW,WU)
              if (IPI2.GT.NM) then
                do J = 1, M
                  W3(J) = 0.
                  W2(J) = 0.
                enddo
              else
                call PRDCT(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),izero,
     &                      DUM, Y(1,IPI2),W3,M,AM,BM,CM,WD,WW,WU)
                call PRDCT(NP1,B(IP1),izero,DUM,izero,DUM,NC,
     &                     CN(IDXC), W3,W2,M, AM,BM,CM,WD,WW,WU)
              endif 
              do J = 1, M
                Y(J,I) = W1(J) + W2(J) + Y(J,I)
              enddo
            endif
          endif
        enddo
      enddo
      if (NPP.EQ.0) then
!
! the periodic case is treated using the capacitance matrix method
!
        IF = 2 ** K
        I = IF / 2
        I1 = I / 2
        call INDXB(I-I1,K-2,IM1,NM1)
        call INDXB(I+I1,K-2,IP1,NP1)
        call INDXB(I,K-1,IZ,NZ)
        call PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),izero,DUM,Y(1,I),
     &                              W1,M,AM, BM,CM,WD,WW,WU)
        IZR = I
        do J = 1, M
          W2(J) = W1(J)
        enddo
        do LL = 2, K
          L = K - LL + 1
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I = I2
          call INDXC(I,IR,IDXC,NC)
          call INDXB(I,IR,IZ,NZ)
          call INDXB(I-I1,IR-1,IM1,NM1)
          call INDXB(I+I1,IR-1,IP1,NP1)
          call PRDCT (NP1,B(IP1),izero,DUM,izero,DUM,NC,CN(IDXC),
     &                               W1,W1,M,AM,BM,CM,WD,WW,WU)
          do J = 1, M
            W1(J) = Y(J,I) + W1(J)
          enddo
          call PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),izero,DUM,
     &                               W1,W1,M,AM,BM,CM,WD,WW,WU)
        enddo
        do LL = 2, K
          L = K - LL + 1
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I4 = I2 + I2
          IFD = IF - I2
          do I = I2, IFD, I4
            if (I-I2-IZR.EQ.0) then
              if (I.GT.NM) goto 100
              call INDXA(I,IR,IDXA,NA)
              call INDXB(I,IR,IZ,NZ)
              call INDXB(I-I1,IR-1,IM1,NM1)
              call INDXB(I+I1,IR-1,IP1,NP1)
              call PRDCT (NM1,B(IM1), izero,DUM, izero, DUM, NA,
     &                     AN(IDXA), W2,W2,M,AM, BM,CM,WD,WW,WU)
              do J = 1, M
                W2(J) = Y(J,I) + W2(J)
              enddo
              call PRDCT (NZ,B(IZ), NM1,B(IM1), NP1,B(IP1),izero,
     &                        DUM,W2,W2,M, AM,BM,CM,WD,WW,WU)
              IZR = I
              if (I.EQ.NM) goto 110
            endif
          enddo
  100    continue
        enddo
       
  110   do J = 1, M
          Y(J,NM+1) = Y(J,NM+1) - CN(NM+1) * W1(J) - AN(NM+1) * W2(J)
        enddo
        call INDXB(IF/2,K-1,IM1,NM1)
        call INDXB(IF,K-1,IP,NP)
        if (NCMPLX.NE.0) then
          call CPRDCT(NM+1,B(IP),NM1,B(IM1),izero,DUM,izero,DUM,
     &                  Y(1,NM+1),Y(1,NM+1),M,AM,BM,CM,W1,W3,WW)
        else
          call PRDCT(NM+1,B(IP),NM1,B(IM1),izero,DUM,izero,DUM,
     &                  Y(1,NM+1),Y(1,NM+1),M,AM,BM,CM,WD,WW,WU)
        endif
        do J = 1, M
          W1(J) = AN(1) * Y(J,NM+1)
          W2(J) = CN(NM) * Y(J,NM+1)
          Y(J,1) = Y(J,1) - W1(J)
          Y(J,NM) = Y(J,NM) - W2(J)
        enddo
        do L = 1, KDO
          IR = L - 1
          I2 = 2 ** IR
          I4 = I2 + I2
          I1 = I2 / 2
          I = I4
          call INDXA(I,IR,IDXA,NA)
          call INDXB(I-I2,IR,IM2,NM2)
          call INDXB(I-I2-I1,IR-1,IM3,NM3)
          call INDXB(I-I1,IR-1,IM1,NM1)
          call PRDCT (NM2,B(IM2),NM3,B(IM3),NM1,B(IM1),izero,DUM,
     &                                W1,W1,M,AM,BM,CM,WD,WW,WU)
          call PRDCT(NM1,B(IM1),izero,DUM,izero,DUM,NA,AN(IDXA),
     &                                 W1,W1,M,AM,BM,CM,WD,WW,WU)
          do J = 1, M
            Y(J,I) = Y(J,I) - W1(J)
          enddo
        enddo

        IZR = NM
        do L = 1, KDO
          IR = L - 1
          I2 = 2 ** IR
          I1 = I2 / 2
          I3 = I2 + I1
          I4 = I2 + I2
          IRM1 = IR - 1
          do I = I4, IF, I4
            IPI1 = I + I1
            IPI2 = I + I2
            IPI3 = I + I3
            if (IPI2.NE.IZR) then
              If (I-IZR) 170, 180, 170
            endif 
            call INDXC(I,IR,IDXC,NC)
            call INDXB(IPI2,IR,IP2,NP2)
            call INDXB(IPI1,IRM1,IP1,NP1)
            call INDXB(IPI3,IRM1,IP3,NP3)
            call PRDCT(NP2,B(IP2),NP1,B(IP1),NP3,B(IP3),izero,DUM,
     &                                 W2,W2,M, AM,BM,CM,WD,WW,WU)
            call PRDCT(NP1,B(IP1),izero,DUM,izero,DUM,NC,CN(IDXC),
     &                                  W2,W2,M,AM,BM,CM,WD,WW,WU)
            do J = 1, M
              Y(J,I) = Y(J,I) - W2(J)
            enddo
            IZR = I
            goto 180
  170      continue 
          enddo
  180    continue
        enddo
      endif 
!
! begin back substitution phase
!
      do LL = 1, K
        L = K - LL + 1
        IR = L - 1
        IRM1 = IR - 1
        I2 = 2 ** IR
        I1 = I2 / 2
        I4 = I2 + I2
        IFD = IF - I2
        do I = I2, IFD, I4
          if (I.LE.NM) then
            IMI1 = I - I1
            IMI2 = I - I2
            IPI1 = I + I1
            IPI2 = I + I2
            call INDXA(I,IR,IDXA,NA)
            call INDXC(I,IR,IDXC,NC)
            call INDXB(I,IR,IZ,NZ)
            call INDXB(IMI1,IRM1,IM1,NM1)
            call INDXB(IPI1,IRM1,IP1,NP1)
            if (I.LE.I2) then
              do J = 1, M
                W1(J) = 0.
              enddo
            else
              call PRDCT (NM1,B(IM1),izero,DUM,izero,DUM,NA,AN(IDXA),
     &                           Y(1,IMI2), W1,M,AM,BM,CM,WD,WW,WU)
            endif
            if (IPI2.GT.NM) then
              do J = 1, M
                W2(J) = 0.
              enddo
            else
              call PRDCT (NP1,B(IP1),izero,DUM,izero,DUM,NC,CN(IDXC),
     &                     Y(1,IPI2), W2,M,AM,BM,CM,WD,WW,WU)
            endif 
            do J = 1, M
              W1(J) = Y(J,I) + W1(J) + W2(J)
            enddo
            call PRDCT (NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),izero,DUM,
     &                            W1,Y(1,I),M, AM,BM,CM,WD,WW,WU)
          endif
        enddo
      enddo
      return
      end



      subroutine INDXB (I,IR,IDX,IDP)
!
! b(idx) is the location of the first root of the b(i,ir) polynomial
!
      implicit none
      INT_TYPE I,IR,IDX,IDP
#include "cblkt.h"
      INT_TYPE IZH,ID,IPL

      IDP = 0
      If (IR) 30, 10, 20
   10 If (I.GT.NM) goto 30
      IDX = I
      IDP = 1
      return
   20 IZH = 2 ** IR
      ID = I - IZH - IZH
      IDX = ID + ID + (IR-1) * IK + IR + (IK-I) / IZH + 4
      IPL = IZH - 1
      IDP = IZH + IZH - 1
      If (I-IPL-NM.GT.0) then
        IDP = 0
        return
      endif 
      if (I+IPL-NM.GT.0) then
        IDP = NM + IPL - I + 1
      endif 
   30 return
      end 

      subroutine INDXA (I,IR,IDXA,NA)
      implicit none
      INT_TYPE I,IR,IDXA,NA
#include "cblkt.h"
      NA = 2 ** IR
      IDXA = I - NA + 1
      if (I.gt.NM) NA = 0
      return
      end

      subroutine INDXC (I,IR,IDXC,NC)
      implicit none
      INT_TYPE I,IR,IDXC,NC
#include "cblkt.h"
      NC = 2 ** IR
      IDXC = I
      if (IDXC+NC-1-NM.GT.0) NC = 0
      return
      end




      subroutine PROD(ND,BD, NM1,BM1, NM2,BM2, NA,AA, X,Y,
     &                     M, A,B,C,D,W,U)
!
! prod applies a sequence of matrix operations to the vector x and
! stores the result in y
!
! bd,bm1,bm2 are arrays containing roots of certian b polynomials
! nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
! aa   array containing scalar multipliers of the vector x
! na is the length of the array aa
! x,y  the matrix operations are applied to x and the result is y
! a,b,c  are arrays which contain the tridiagonal matrix
! m  is the order of the matrix
! d,w,u are working arrays
! is  determines whether or not a change in sign is made
!
      implicit none
      INT_TYPE ND, NM1, NM2, NA, M
      REAL_TYPE BD(ND), BM1(NM1), BM2(NM2),  AA(NA),
     &     X(M),Y(M), A(M),B(M),C(M), D(M),W(M),U(M)

!-- internal -->

      REAL_TYPE DEN, RT
      INT_TYPE J,K, IA,ID,IBR, M1,M2,MM
#if VERBOSE > 2
      write(*,*) 'PROD'
#endif

      do J = 1, M
        W(J) = X(J)
        Y(J) = W(J)
      enddo 
      MM = M - 1
      ID = ND
      IBR = 0
      M1 = NM1
      M2 = NM2
      IA = NA
   20 if (IA.GT.0) then
        RT = AA(IA)
        If (ND.EQ.0) RT = -RT
        IA = IA - 1
        do J = 1, M              ! scalar multiplication
          Y(J) = RT * W(J)
        enddo
      endif 
      if (ID.GT.0) then
        RT = BD(ID)
        ID = ID - 1
        If (ID.EQ.0) IBR = 1
!
! begin solution to system
!
        D(M) = A(M) / (B(M)-RT)
        W(M) = Y(M) / (B(M)-RT)
        do J = 2, MM
          K = M - J
          DEN = B(K+1) - RT - C(K+1) * D(K+2)
          D(K+1) = A(K+1) / DEN
          W(K+1) = (Y(K+1)-C(K+1)*W(K+2)) / DEN
        enddo
        DEN = B(1) - RT - C(1) * D(2)
        W(1) = 1.
        If (DEN.NE.0) then
          W(1) = (Y(1)-C(1)*W(2)) / DEN
        endif 
        do J = 2, M
          W(J) = W(J) - D(J) * W(J-1)
        enddo
        If (NA) 80, 80, 20
   60   do J = 1, M
          Y(J) = W(J)
        enddo
        IBR = 1
        Go To 20
   80   If (M1.LE.0) then
          If (M2) 60, 60, 90
        endif 
        If (M2.GT.0) then
          If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90
        endif 
        If (IBR.LE.0) then
          If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60
        endif 
        RT = RT - BM1(M1)
        M1 = M1 - 1
        Go To 100
   90   If (IBR.LE.0) then
          If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60
        endif 
        RT = RT - BM2(M2)
        M2 = M2 - 1
  100   do J = 1, M
          Y(J) = Y(J) + RT * W(J)
        enddo
       goto 20
      endif
      return
      end




      subroutine PRODP (ND,BD, NM1,BM1, NM2,BM2, NA,AA,
     &                           X,Y, M, A,B,C,D, U,W)
!
! prodp applies a sequence of matrix operations to the vector x and
! stores the result in y        periodic boundary conditions
!
! bd,bm1,bm2 are arrays containing roots of certian b polynomials
! nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
! aa   array containing scalar multipliers of the vector x
! na is the length of the array aa
! x,y  the matrix operations are applied to x and the result is y
! a,b,c  are arrays which contain the tridiagonal matrix
! m  is the order of the matrix
! d,u,w are working arrays
!
      implicit none
      INT_TYPE ND, NM1, NM2, NA, M
      REAL_TYPE BD(ND), BM1(NM2), BM2(NM2), AA(NA),
     &    X(M),Y(M), A(M),B(M),C(M), D(M),U(M),W(M)

!-- internal -->

      REAL_TYPE RT,BH,YM,DEN,V,AM
      INT_TYPE J,K, IA,ID,IBR, M1,M2,MM,MM2

#if VERBOSE > 2
      write(*,*) 'PRODP'
#endif
      do J = 1, M
        Y(J) = X(J)
        W(J) = Y(J)
      enddo
      MM = M - 1
      MM2 = M - 2
      ID = ND
      IBR = 0
      M1 = NM1
      M2 = NM2
      IA = NA
   20 if (IA.GT.0) then
        RT = AA(IA)
        If (ND.EQ.0) RT = -RT
        IA = IA - 1
        do J = 1, M
          Y(J) = RT * W(J)
        enddo
      endif
      If (ID.GT.0) then
        RT = BD(ID)
        ID = ID - 1
        If (ID.EQ.0) IBR = 1
!
! begin solution to system
!
        BH = B(M) - RT
        YM = Y(M)
        DEN = B(1) - RT
        D(1) = C(1) / DEN
        U(1) = A(1) / DEN
        W(1) = Y(1) / DEN
        V = C(M)
        If (MM2.GE.2) then
          do J = 2, MM2
            DEN = B(J) - RT - A(J) * D(J-1)
            D(J) = C(J) / DEN
            U(J) = -A(J) * U(J-1) / DEN
            W(J) = (Y(J)-A(J)*W(J-1)) / DEN
            BH = BH - V * U(J-1)
            YM = YM - V * W(J-1)
            V = -V * D(J-1)
          enddo
        endif 
        DEN = B(M-1) - RT - A(M-1) * D(M-2)
        D(M-1) = (C(M-1)-A(M-1)*U(M-2)) / DEN
        W(M-1) = (Y(M-1)-A(M-1)*W(M-2)) / DEN
        AM = A(M) - V * D(M-2)
        BH = BH - V * U(M-2)
        YM = YM - V * W(M-2)
        DEN = BH - AM * D(M-1)
        if (DEN.NE.0.) then
          W(M) = (YM-AM*W(M-1)) / DEN
        else
          W(M) = 1.
        endif 
        W(M-1) = W(M-1) - D(M-1) * W(M)
        do J = 2, MM
          K = M - J
          W(K) = W(K) - D(K) * W(K+1) - U(K) * W(M)
        enddo
        If (NA) 80, 80, 20
   60   continue 
        do J = 1, M
          Y(J) = W(J)
        enddo
        IBR = 1
        Go To 20
   80   If (M1.LE.0) then
          If (M2) 60, 60, 90
        endif 
        If (M2.GT.0) then
          If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90
        endif 
        If (IBR.LE.0) then
          If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60
        endif 
        RT = RT - BM1(M1)
        M1 = M1 - 1
        goto 100
   90   If (IBR.LE.0) then
          If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60
        endif 
        RT = RT - BM2(M2)
        M2 = M2 - 1
  100   continue
        do J = 1, M
          Y(J) = Y(J) + RT * W(J)
        enddo
        goto 20
      endif 
      return
      end




      subroutine CPROD (ND,BD, NM1,BM1, NM2,BM2, NA,AA,
     &                       X,YY, M, A,B,C, D,W,Y)
!
! prod applies a sequence of matrix operations to the vector x and
! stores the result in yy           (complex case)
! aa   array containing scalar multipliers of the vector x
! nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
! bd,bm1,bm2 are arrays containing roots of certian b polynomials
! na is the length of the array aa
! x,yy the matrix operations are applied to x and the result is yy
! a,b,c  are arrays which contain the tridiagonal matrix
! m  is the order of the matrix
! d,w,y are working arrays
!
      implicit none
      INT_TYPE ND, NM1, NM2, NA, M
      COMPLEX_TYPE BD(ND), Y(M), D(M), W(M)
      REAL_TYPE BM1(NM1), BM2(NM2), AA(NA),
     &         X(M), YY(M), A(M), B(M), C(M)

!-- internal -->

      INT_TYPE J,K, IA,ID,M1,M2,MM, IFLG 
      COMPLEX_TYPE CRT, DEN, Y1, Y2
      REAL_TYPE RT, zero
      parameter (zero=0.)
#if VERBOSE > 2
      write(*,*) 'CPROD'
#endif
      do J = 1, M
        Y(J) = CMPLX(X(J),zero)
      enddo
      MM = M - 1
      ID = ND
      M1 = NM1
      M2 = NM2
      IA = NA
   20 IFLG = 0
      If (ID.GT.0) then
        CRT = BD(ID)
        ID = ID - 1
c
c begin solution to system
c
        D(M) = A(M) / (B(M)-CRT)
        W(M) = Y(M) / (B(M)-CRT)
        do J = 2, MM
          K = M - J
          DEN = B(K+1) - CRT - C(K+1) * D(K+2)
          D(K+1) = A(K+1) / DEN
          W(K+1) = (Y(K+1)-C(K+1)*W(K+2)) / DEN
        enddo
        DEN = B(1) - CRT - C(1) * D(2)
        If (CABS(DEN).NE.0) then
          Y(1) = (Y(1)-C(1)*W(2)) / DEN
        else
          Y(1) = (1.,0.)
        endif 
        do J = 2, M
          Y(J) = W(J) - D(J) * Y(J-1)
        enddo
      endif 
      if (M1.LE.0) then
        if (M2.LE.0) Go To 60
        RT = BM2(M2)
        M2 = M2 - 1
      else
        If (M2.LE.0) then
          RT = BM1(M1)
          M1 = M1 - 1
        else
          If (ABS(BM1(M1)).GT.ABS(BM2(M2))) then
            RT = BM1(M1)
            M1 = M1 - 1
          else
            RT = BM2(M2)
            M2 = M2 - 1
          endif 
        endif 
      endif 
      Y1 = (B(1)-RT) * Y(1) + C(1) * Y(2)
      if (MM.GE.2) then    !--> matrix multiplication
        do J = 2, MM
          Y2 = A(J) * Y(J-1) + (B(J)-RT) * Y(J) + C(J) * Y(J+1)
          Y(J-1) = Y1
          Y1 = Y2
        enddo
      endif 
      Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M)
      Y(M-1) = Y1
      IFLG = 1
      Go To 20
   60 if (IA.GT.0) then
        RT = AA(IA)
        IA = IA - 1
        IFLG = 1
        do J = 1, M       !--> scalar multiplication
          Y(J) = RT * Y(J)
        enddo
      endif 
      if (IFLG.GT.0) goto 20
      do J = 1, M
        YY(J) = REAL_PART(Y(J))
      enddo
      return
      end




      subroutine CPRODP (ND,BD, NM1,BM1, NM2,BM2, NA,AA,
     &                           X,YY, M, A,B,C,D,U,Y)
!
! prodp applies a sequence of matrix operations to the vector x and
! stores the result in yy       periodic boundary conditions
! and  complex  case
!
! bd,bm1,bm2 are arrays containing roots of certian b polynomials
! nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively
! aa   array containing scalar multipliers of the vector x
! na is the length of the array aa
! x,yy the matrix operations are applied to x and the result is yy
! a,b,c  are arrays which contain the tridiagonal matrix
! m  is the order of the matrix
! d,u,y are working arrays
!
      implicit none
      INT_TYPE ND, NM1, NM2, NA, M
      COMPLEX_TYPE BD(ND), D(M), U(M), Y(M)
      REAL_TYPE BM1(NM1), BM2(NM2), AA(NA),
     &          X(M), YY(M),  A(M), B(M), C(M)

!-- internal -->

      INT_TYPE J,K, IA,ID, M1,M2,MM,MM2, IFLG
      COMPLEX_TYPE V, DEN, BH, YM, AM, Y1, Y2, YH, CRT
      REAL_TYPE RT, zero
      parameter (zero=0.D0)

#if VERBOSE > 2
      write(*,*) 'CPRODP'
#endif

      do J = 1, M
        Y(J) = CMPLX(X(J),zero)
      enddo
      MM = M - 1
      MM2 = M - 2
      ID = ND
      M1 = NM1
      M2 = NM2
      IA = NA
   20 IFLG = 0
      If (ID.GT.0) then
        CRT = BD(ID)         !--> begin solution to system
        ID = ID - 1
        IFLG = 1
        BH = B(M) - CRT
        YM = Y(M)
        DEN = B(1) - CRT
        D(1) = C(1) / DEN
        U(1) = A(1) / DEN
        Y(1) = Y(1) / DEN
        V = CMPLX(C(M),0.)
        If (MM2.GE.2) then
          do J = 2, MM2
            DEN = B(J) - CRT - A(J) * D(J-1)
            D(J) = C(J) / DEN
            U(J) = -A(J) * U(J-1) / DEN
            Y(J) = (Y(J)-A(J)*Y(J-1)) / DEN
            BH = BH - V * U(J-1)
            YM = YM - V * Y(J-1)
            V = -V * D(J-1)
          enddo
        endif 
        DEN = B(M-1) - CRT - A(M-1) * D(M-2)
        D(M-1) = (C(M-1)-A(M-1)*U(M-2)) / DEN
        Y(M-1) = (Y(M-1)-A(M-1)*Y(M-2)) / DEN
        AM = A(M) - V * D(M-2)
        BH = BH - V * U(M-2)
        YM = YM - V * Y(M-2)
        DEN = BH - AM * D(M-1)
        If (CABS(DEN).NE.0) then
          Y(M) = (YM-AM*Y(M-1)) / DEN
        else
          Y(M) = (1.,0.)
        endif 
        Y(M-1) = Y(M-1) - D(M-1) * Y(M)
        do J = 2, MM
          K = M - J
          Y(K) = Y(K) - D(K) * Y(K+1) - U(K) * Y(M)
        enddo
      endif 
      If (M1.LE.0) then
        If (M2.LE.0) Go To 60
        RT = BM2(M2)
        M2 = M2 - 1
      else
        If (M2.LE.0) then
          RT = BM1(M1)
          M1 = M1 - 1
        else
          If (ABS(BM1(M1)).GT.ABS(BM2(M2))) then
            RT = BM1(M1)
            M1 = M1 - 1
          else
            RT = BM2(M2)
            M2 = M2 - 1
          endif 
        endif 
      endif 
!
! matrix multiplication
!
      YH = Y(1)
      Y1 = (B(1)-RT)*Y(1) + C(1)*Y(2) + A(1)*Y(M)
      If (MM.GE.2) then
        do J = 2, MM
          Y2 = A(J)*Y(J-1) + (B(J)-RT)*Y(J) + C(J)*Y(J+1)
          Y(J-1) = Y1
          Y1 = Y2
        enddo
      endif 
      Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M) + C(M) * YH
      Y(M-1) = Y1
      IFLG = 1
      Go To 20
   60 if (IA.GT.0) then
        RT = AA(IA)
        IA = IA - 1
        IFLG = 1
        do J = 1, M           ! scalar multiplication
          Y(J) = RT * Y(J)
        enddo
      endif 
      If (IFLG.GT.0) Go To 20
      do J = 1, M
        YY(J) = REAL_PART(Y(J))
      enddo
      return
      end




      subroutine PPADD (N,IERROR,A,C,CBP,BP,BH)
!
! ppadd computes the eigenvalues of the periodic tridiagonal matrix
!     with coefficients an,bn,cn
! n is the order of the bh and bp polynomials;
! on output bp contians the eigenvalues
! cbp is the same as bp except type complex
! bh is used to temporarily store the roots of the b hat polynomial
! which enters through bp
!
      implicit none
      INT_TYPE N,IERROR
      REAL_TYPE A(N), C(N), BP(N), BH(N)
      COMPLEX_TYPE CBP(N)
#include "cblkt.h"
   
      INT_TYPE IZ,IZM,IZM2, ICV, IS,IG,I2,I3,IF, IT,NT,MODIZ, J,NHALF
      COMPLEX_TYPE CX, FSG, HSG, DD, F, FP, FPP, CDIS, R1, R2, R3
      REAL_TYPE SCNV, XL,XR,XM, DB, SGN, PSG

      REAL_TYPE BSRH, PSGF, PPSGF, PPSPF  ! these are functions
      external PSGF, PPSGF, PPSPF
#ifdef VERBOSE
      write(*,*) 'PPADD'
#endif
      SCNV = SQRT(CNV)
      IZ = N
      IZM = IZ - 1
      IZM2 = IZ - 2
      If (BP(N)-BP(1)) 10, 240, 30
   10 do J = 1, N
        NT = N - J
        BH(J) = BP(NT+1)
      enddo
      goto 50
   30 do J = 1, N
        BH(J) = BP(J)
      enddo
   50 NCMPLX = 0
      MODIZ = MOD(IZ,2)
      IS = 1
      If (MODIZ.NE.0) then
        If (A(1)) 80, 240, 60
      endif 
   60 XL = BH(1)
      DB = BH(3) - BH(1)
   70 XL = XL - DB
      If (PSGF(XL,IZ,C,A,BH).LE.0) Go To 70
      SGN = -1.
      CBP(1) = CMPLX(BSRH(XL,BH(1),IZ,C,A,BH,PSGF,SGN),0.)
      IS = 2
   80 IF = IZ - 1
      If (MODIZ.NE.0) then
        If (A(1)) 90, 240, 110
      endif 
   90 XR = BH(IZ)
      DB = BH(IZ) - BH(IZ-2)
  100 XR = XR + DB
      If (PSGF(XR,IZ,C,A,BH).LT.0) Go To 100
      SGN = 1.
      CBP(IZ) = CMPLX(BSRH(BH(IZ),XR,IZ,C,A,BH,PSGF,SGN),0.)
      IF = IZ - 2
  110 Do 180 IG = IS, IF, 2
        XL = BH(IG)
        XR = BH(IG+1)
        SGN = -1.
        XM = BSRH(XL,XR,IZ,C,A,BH,PPSPF,SGN)
        PSG = PSGF(XM,IZ,C,A,BH)
        If (ABS(PSG).GT.EPS) then
          If (PSG*PPSGF(XM,IZ,C,A,BH)) 120, 130, 140
!
! case of a real zero
!
  120     SGN = 1.
          CBP(IG) = CMPLX(BSRH(BH(IG),XM,IZ,C,A,BH,PSGF,SGN),0.)
          SGN = -1.
          CBP(IG+1) = CMPLX(BSRH(XM,BH(IG+1),IZ,C,A,BH,PSGF,SGN),0.)
          Go To 180
        endif 
!
! case of a multiple zero
!
  130   CBP(IG) = CMPLX(XM,0.)
        CBP(IG+1) = CMPLX(XM,0.)
        goto 180
!
! case of a complex zero
!
  140   IT = 0
        ICV = 0
        CX = CMPLX(XM,0.)
  150   FSG = (1.,0.)
        HSG = (1.,0.)
        FP = (0.,0.)
        FPP = (0.,0.)
        do J = 1, IZ
          DD = 1. / (CX-BH(J))
          FSG = FSG * A(J) * DD
          HSG = HSG * C(J) * DD
          FP = FP + DD
          FPP = FPP - DD * DD
        enddo
        if (MODIZ.EQ.0) then
          F = (1.,0.) - FSG - HSG
        else
          F = (1.,0.) + FSG + HSG
        endif 
        I3 = 0
        If (CABS(FP).GT.0) then
          I3 = 1
          R3 = -F / FP
        endif 
        I2 = 0
        if (CABS(FPP).GT.0) then
          I2 = 1
          CDIS = CSQRT(FP**2-2.*F*FPP)
          R1 = CDIS - FP
          R2 = -FP - CDIS
          if (CABS(R1).GT.CABS(R2)) then
            R1 = R1 / FPP
          else
            R1 = R2 / FPP
          endif 
          R2 = 2. * F / FPP / R1
          If (CABS(R2).LT.CABS(R1)) R1 = R2
          If (I3.LE.0) Go To 170
          If (CABS(R3).LT.CABS(R1)) R1 = R3
        else
          R1 = R3
        endif 
  170   CX = CX + R1
        IT = IT + 1
        If (IT.GT.50) Go To 240
        If (CABS(R1).GT.SCNV) Go To 150
        If (ICV.LE.0) then
          ICV = 1
          Go To 150
        endif 
        CBP(IG) = CX
        CBP(IG+1) = CONJG(CX)
  180 Continue
      If (CABS(CBP(N))-CABS(CBP(1))) 190, 240, 210
  190 NHALF = N / 2
      do J = 1, NHALF
        NT = N - J
        CX = CBP(J)
        CBP(J) = CBP(NT+1)
        CBP(NT+1) = CX
      enddo
  210 NCMPLX = 1
      do J = 2, IZ
        if (IMAG_PART(CBP(J)).NE.0) Go To 250
      enddo
      NCMPLX = 0
      do J = 2, IZ
        BP(J) = REAL_PART(CBP(J))
      enddo
      goto 250
  240 IERROR = 4
  250 continue
      return
      end


      function BSRH (XLL,XRR,IZ,C,A,BH,F,SGN)
      implicit none
      INT_TYPE IZ 
      REAL_TYPE BSRH, XLL,XRR, C(1),A(1),BH(1), F,SGN, XL,XR,DX,X
      external F
#include "cblkt.h"
      XL = XLL
      XR = XRR
      DX = .5 * ABS(XR-XL)
   10 X = .5 * (XL+XR)
      If (SGN*F(X,IZ,C,A,BH)) 30, 50, 20
   20 XR = X
      Go To 40
   30 XL = X
   40 DX = .5 * DX
      If (DX.GT.CNV) Go To 10
   50 BSRH = .5 * (XL+XR)
      return
      end


      function PSGF(X,IZ,C,A,BH)
      implicit none
      INT_TYPE IZ, J
      REAL_TYPE PSGF, X, A(1), C(1), BH(1), FSG,HSG,DD
      FSG = 1.
      HSG = 1.
      do J = 1, IZ
        DD = 1. / (X-BH(J))
        FSG = FSG * A(J) * DD
        HSG = HSG * C(J) * DD
      enddo
      if (MOD(IZ,2).EQ.0) then
        PSGF = 1. - FSG - HSG
      else
        PSGF = 1. + FSG + HSG
      endif
      return
      end

      function PPSGF (X,IZ,C,A,BH)
      implicit none
      INT_TYPE IZ, J
      REAL_TYPE PPSGF, X, A(1), C(1), BH(1), SUM
      SUM = 0.
      do J = 1, IZ
        SUM = SUM - 1. / (X-BH(J)) ** 2
      enddo
      PPSGF = SUM
      return
      end


      function PPSPF (X,IZ,C,A,BH)
      implicit none
      INT_TYPE IZ, J
      REAL_TYPE PPSPF, X, A(1), C(1), BH(1), SUM
      SUM = 0.
      do J = 1, IZ
        SUM = SUM + 1. / (X-BH(J))
      enddo
      PPSPF = SUM
      return
      end




      subroutine COMPB (N,IERROR, AN,BN,CN,B,AH,BH)
!
! compb computes the roots of the b polynomials using subroutine
! tevls which is a modification the eispack program tqlrat.
! ierror is set to 4 if either tevls fails or if a(j+1)*c(j) is
! less than zero for some j.  ah,bh are temporary work arrays.
!
      implicit none
      INT_TYPE N,IERROR
      REAL_TYPE AN(1), BN(1), CN(1), B(1), AH(1), BH(1)
#include "cblkt.h"
      INT_TYPE I,I2,I4,IB,IF,IFD,IPL,IR, J,J1,J2,JF,JS,
     &                L,L1,L2,LH,LS, KDO, NB,NMP,N2M2
      REAL_TYPE DUM, BNORM, ARG, D1,D2,D3
      REAL_TYPE EPMACH

#ifdef VERBOSE
      write(*,*) 'COMPB'
#endif
      EPS = EPMACH(DUM)
      BNORM = ABS(BN(1))
      do J = 2, NM
        BNORM = max(BNORM,ABS(BN(J)))
        ARG = AN(J) * CN(J-1)
        If (ARG.LT.0) Go To 110
        B(J) = SIGN(SQRT(ARG),AN(J))
      enddo
      CNV = EPS * BNORM
      IF = 2 ** K
      KDO = K - 1
      do L = 1, KDO
        IR = L - 1
        I2 = 2 ** IR
        I4 = I2 + I2
        IPL = I4 - 1
        IFD = IF - I4
        do I = I4, IFD, I4
          call INDXB(I,L,IB,NB)
          If (NB.LE.0) goto 50
          JS = I - IPL
          JF = JS + NB - 1
          LS = 0
          do J = JS, JF
            LS = LS + 1
            BH(LS) = BN(J)
            AH(LS) = B(J)
          enddo
          call TEVLS(NB,BH,AH,IERROR)
          If (IERROR.NE.0) Go To 100
          LH = IB - 1
          do J = 1, NB
            LH = LH + 1
            B(LH) = -BH(J)
          enddo
        enddo
   50  continue
      enddo

      do J = 1, NM
        B(J) = -BN(J)
      enddo
      if (NPP.EQ.0) then
        NMP = NM + 1
        NB = NM + NMP
        do J = 1, NB
          L1 = MOD(J-1,NMP) + 1
          L2 = MOD(J+NM-1,NMP) + 1
          ARG = AN(L1) * CN(L2)
          If (ARG.LT.0) Go To 110
          BH(J) = SIGN(SQRT(ARG),-AN(L1))
          AH(J) = -BN(L1)
        enddo
        call TEVLS(NB,AH,BH,IERROR)
        If (IERROR.NE.0) Go To 100
        call INDXB(IF,K-1,J2,LH)
        call INDXB(IF/2,K-1,J1,LH)
        J2 = J2 + 1
        LH = J2
        N2M2 = J2 + NM + NM - 2
   80   D1 = ABS(B(J1)-B(J2-1))
        D2 = ABS(B(J1)-B(J2))
        D3 = ABS(B(J1)-B(J2+1))
        If (.NOT.((D2.LT.D1).AND.(D2.LT.D3))) then
          B(LH) = B(J2)
          J2 = J2 + 1
          LH = LH + 1
          If (J2-N2M2) 80, 80, 90
        endif 
        J2 = J2 + 1
        J1 = J1 + 1
        If (J2.LE.N2M2) Go To 80
   90   B(LH) = B(N2M2+1)
        call INDXB(IF,K-1,J1,J2)
        J2 = J1 + NMP + NMP
        call PPADD(NM+1,IERROR,AN,CN,B(J1),B(J1),B(J2))
      endif 
      return
  100 IERROR = 4
      return
  110 IERROR = 5
      return
      end






      subroutine TEVLS(N, D,E2, IERR)
      implicit none
      INT_TYPE N, IERR
      REAL_TYPE D(N), E2(N)
#include "cblkt.h"
      INT_TYPE I, J, L, M,  II, L1, MML, NTOP,NHALF
      REAL_TYPE B, C, F, G, H, P, R, S, DHOLD
!
! This subroutine is a modification of the eispack subroutine tqlrat
! algorithm 464, comm. acm 16, 689(1973) by reinsch. It finds the
! eigenvalues of a symmetric tridiagonal matrix by the rational ql
! method.
!
! input
!
!     n is the order of the matrix,
!
!     d contains the diagonal elements of the input matrix,
!
!     e2 contains the subdiagonal elements of the input matrix in
!        its last n-1 positions.  e2(1) is arbitrary.
!
! output
!
!     d contains the eigenvalues in ascending order.  if an error
!       exit is made, the eigenvalues are correct and ordered for
!       indices 1,2,...ierr-1, but may not be the smallest
!       eigenvalues,
!
!     e2 has been destroyed,
!
!     ierr is set to
!          zero       for normal return,
!          j          if the j-th eigenvalue has not been
!                     determined after 30 iterations.
!
! questions and comments should be directed to b. s. garbow,
! applied mathematics division, argonne national laboratory
!
!  ********** machep is a machine dependent parameter specifying
!             the relative precision of floating point arithmetic.
!

      IERR = 0
      If (N.NE.1) then

        do I = 2, N
          E2(I-1) = E2(I) * E2(I)
        enddo

        F = 0.0
        B = 0.0
        E2(N) = 0.0

        do L = 1, N
          J = 0
          H = EPS * (ABS(D(L))+SQRT(E2(L)))
          If (B.LE.H) then
            B = H
            C = B * B
          endif 
!
! ********** look for small squared sub-diagonal element **********
!
          do M = L, N
            If (E2(M).LE.C) Go To 30
!
! ********** e2(n) is always zero, so there is no exit
!             through the bottom of the loop **********
!
          enddo

   30     if (M.NE.L) then
   40       if (J.EQ.30) goto 110
            J = J + 1
!
! ********** form shift **********
!
            L1 = L + 1
            S = SQRT(E2(L))
            G = D(L)
            P = (D(L1)-G) / (2.0*S)
            R = SQRT(P*P+1.0)
            D(L) = S / (P+SIGN(R,P))
            H = G - D(L)

            do I = L1, N
              D(I) = D(I) - H
            enddo

            F = F + H
!
! ********** rational ql transformation **********
!
            G = D(M)
            If (G.EQ.0.0) G = B
            H = G
            S = 0.0
            MML = M - L
!
! ********** for i=m-1 step -1 until l do -- **********
!
            do II = 1, MML
              I = M - II
              P = G * H
              R = P + E2(I)
              E2(I+1) = S * R
              S = E2(I) / R
              D(I+1) = H + S * (H+D(I))
              G = D(I) - E2(I) / G
              If (G.EQ.0.0) G = B
              H = G * P / R
            enddo

            E2(L) = S * G
            D(L) = H
!
! ********** guard against underflowed h **********
!
            If (H.NE.0.0) then
              If (ABS(E2(L)).GT.ABS(C/H)) then
                E2(L) = H * E2(L)
                If (E2(L).NE.0.0) Go To 40
              endif 
            endif 
          endif 
          P = D(L) + F
!
! ********** order eigenvalues **********
!
          if (L.NE.1) then
!
! ********** for i=l step -1 until 2 do -- **********
!
            do II = 2, L
              I = L + 2 - II
              If (P.GE.D(I-1)) Go To 80
              D(I) = D(I-1)
            enddo
          endif 

          I = 1
   80     D(I) = P
        enddo

        If (ABS(D(N)).GE.ABS(D(1))) Go To 120
        NHALF = N / 2
        do I = 1, NHALF
          NTOP = N - I
          DHOLD = D(I)
          D(I) = D(NTOP+1)
          D(NTOP+1) = DHOLD
        enddo
        goto 120
!
! ********** set error -- no convergence to an
!            eigenvalue after 30 iterations **********
!
  110   IERR = L
      endif
  120 return
      end





!
! package comf   the entries in this package are lowlevel entries,
!                supporting ulib packages blktri and cblktri.  That
!                is, these routines are not called directly by users,
!                but rather by entries within blktri and cblktri.
!                Description of entries epmach and pimach follow
!                below.
!
! latest revision        january 1985
!
! function epmach (dum)
!
! purpose        to compute an approximate machine accuracy epsilon
!                according to the following definition: epsilon is
!                the smallest number such that (1.+epsilon).gt.1.)
!
! usage                  eps = epmach (dum)
!
! arguments
! on input               dum
!                          dummy value
!
! arguments
! on output              none
!
! history        the original version, written when the blktri
!                package was converted from the cdc 7600 to run on
!                the cray-1, calculated machine accuracy by
!                successive divisions by 10.  use of this constant
!                caused blktri to compute solutions on the cray-1
!                with four fewer places of accuracy than the version
!                on the 7600.  it was found that computing machine
!                accuracy by successive divisions of 2 produced a
!                machine accuracy 29% less than the value generated
!                by successive divisions by 10, and that use of this
!                machine constant in the blktri package recovered
!                the accuracy that appeared to be lost on conversion.
!
! algorithm      computes machine accuracy by successive divisions
!                of two.
!
! portability    this code will execute on machines other than the
!                cray1, but the returned value may be unsatisfactory.
!                see history above.


!
! function pimach (dum)
!
! purpose        to supply the value of the constant pi correct to
!                machine precision where
!                pi=3.141592653589793238462643383279502884197
!                             1693993751058209749446
!
! usage          pi = pimach (dum)
!
! arguments
! on input       dum
!                        dummy value
!
! arguments
! on output      none
!
! algorithm      the value of pi is set in a constant.
!
! portability    this entry is portable, but users should check
!                to see whether greater accuracy is required.

      function EPMACH(DUM)
      implicit none
      REAL_TYPE EPMACH, DUM,V, EPS
      common /VALUE/ V
      EPS = 1.
   10 EPS = EPS / 2.
      call STORE(EPS+1.)
      If (V.GT.1.) goto 10
      EPMACH = 100. * EPS
      return
      end

      subroutine STORE(X)
      implicit none
      REAL_TYPE X,V
      common /VALUE/ V
      V = X
      return
      end

      function PIMACH(DUM)
      implicit none
      REAL_TYPE PIMACH, DUM
! pi=3.1415926535897932384626433832795028841971693993751058209749446
      PIMACH = 3.14159265358979
      return
      end

!
! package sepaux   contains no user entry points.
!
! latest revision  march 1985
!
! purpose          this package contains auxiliary routines for ncar
!                  public software packages such as sepeli and sepx4.
!
! usage            since this package contains no user entries,
!                  no usage instructions or argument descriptions
!                  are given here.
!
! history          developed in the late 1970's by John C. Adams of
!                  NCAR's scienttific computing division.
!

      subroutine SEPORT(USOL,IDMN,ZN,ZM,PERTRB)
!
! this subroutine orthoganalizes the array usol with respect
! to the constant array in a weighted least squares norm
!
      implicit none
      INT_TYPE IDMN
      REAL_TYPE USOL(IDMN,1), ZN(1), ZM(1), PERTRB

#include "splp.h"

      REAL_TYPE UTE,ETE
      INT_TYPE ISTR,IFNL,JSTR,JFNL, I,II, J,JJ

#ifdef VERBOSE
      write(*,*) 'SEPORT'
#endif

      ISTR = IS
      IFNL = MS
      JSTR = JS
      JFNL = NS

      UTE = 0.0               !--> compute weighted inner products
      ETE = 0.0
      do I = IS, MS
        II = I - IS + 1
        do J = JS, NS
          JJ = J - JS + 1
          ETE = ETE + ZM(II) * ZN(JJ)
          UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ)
        enddo
      enddo

      PERTRB = UTE / ETE      !<-- set perturbation parameter

      do I = ISTR, IFNL       !--> subtract off constant pertrb
        do J = JSTR, JFNL
          USOL(I,J) = USOL(I,J) - PERTRB
        enddo
      enddo
      return
      end




      subroutine SEPMIN (USOL,IDMN,ZN,ZM,PERTB)
!
! this subroutine orhtogonalizes the array usol with respect
! to the constant array in a weighted least squares norm
!
      implicit none
      INT_TYPE IDMN
      REAL_TYPE USOL(IDMN,1), ZN(1), ZM(1), PERTB
#include "splp.h"
      REAL_TYPE UTE,ETE, PERTRB
      INT_TYPE ISTR,IFNL,JSTR,JFNL, I,II, J,JJ 
!
! entry at sepmin occurrs when the final solution is to be minimized
! with respect to the weighted least squares norm
!
#ifdef VERBOSE
      write(*,*) 'SEPMIN'
#endif

      ISTR = 1
      IFNL = K
      JSTR = 1
      JFNL = L
                              !--> Compute weighted inner products
      UTE = 0.0
      ETE = 0.0
      do I = IS, MS
        II = I - IS + 1
        do J = JS, NS
          JJ = J - JS + 1
          ETE = ETE + ZM(II) * ZN(JJ)
          UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ)
        enddo
      enddo

      PERTRB = UTE / ETE      !<-- set perturbation parameter

      do I = ISTR, IFNL       !--> subtract off constant pertrb
        do J = JSTR, JFNL
          USOL(I,J) = USOL(I,J) - PERTRB
        enddo
      enddo
      return
      end





      subroutine SEPTRI (N, A,B,C,D,U,Z)
!
! this subroutine solves for a non-zero eigenvector corresponding
! to the zero eigenvalue of the transpose of the rank deficient one
! matrix with subdiagonal a, diagonal b, and superdiagonal c , with
! a(1) in the (1,n) position, with c(n) in the (n,1) position, and
! all other elements zero.
!
      implicit none
      INT_TYPE N
      REAL_TYPE A(N), B(N), C(N), D(N), U(N), Z(N)

!-- internal -->

      INT_TYPE J,K, NM1,NM2
      REAL_TYPE AN, BN, V, DEN
#ifdef VERBOSE
      write(*,*) 'SEPTRI'
#endif
      BN = B(N)
      D(1) = A(2) / B(1)
      V = A(1)
      U(1) = C(N) / B(1)
      NM2 = N - 2
      do J = 2, NM2
        DEN = B(J) - C(J-1) * D(J-1)
        D(J) = A(J+1) / DEN
        U(J) = -C(J-1) * U(J-1) / DEN
        BN = BN - V * U(J-1)
        V = -V * D(J-1)
      enddo
      DEN = B(N-1) - C(N-2) * D(N-2)
      D(N-1) = (A(N)-C(N-2)*U(N-2)) / DEN
      AN = C(N-1) - V * D(N-2)
      BN = BN - V * U(N-2)
      DEN = BN - AN * D(N-1)

      Z(N) = 1.0           !<-- set the last component equal to one
      Z(N-1) = -D(N-1)
      NM1 = N - 1
      do J = 2, NM1
        K = N - J
        Z(K) = -D(K) * Z(K+1) - U(K) * Z(N)
      enddo
      return
      end



      subroutine SEPDX(U, IDMN,I,J, UXXX,UXXXX)
!
! Compute second-order finite difference approximations to the third
! and fourth x partial derivatives of u at the (i,j) mesh point
!
      implicit none
      INT_TYPE IDMN, I,J
      REAL_TYPE U(IDMN,1), UXXX,UXXXX
#include "splp.h"

      If (I.LE.2.OR.I.GE.(K-1)) then
        If (I.NE.1) then
          If (I.EQ.2) Go To 10
          If (I.EQ.K-1) Go To 20
          If (I.EQ.K) Go To 30
        endif 
!
! compute partial derivative approximations at x=a
!
        if (KSWX.NE.1) then
          UXXX = ( -5.0*U(1,J)+18.0*U(2,J)-24.0*U(3,J)+14.0*U(4,J)
     &                            -3.0*U(5,J)) / (TDLX3)
          UXXXX = ( 3.0*U(1,J)-14.0*U(2,J)+26.0*U(3,J)-24.0*U(4,J)
     &                           +11.0*U(5,J)-2.0*U(6,J)) / DLX4
          return
        endif 
!
! periodic at x=a
!
        UXXX = (-U(K-2,J)+2.0*U(K-1,J)-2.0*U(2,J)+U(3,J)) / (TDLX3)
        UXXXX = (U(K-2,J)-4.0*U(K-1,J)+6.0*U(1,J)-4.0*U(2,J)
     &                                             +U(3,J)) / DLX4
        return
!
! compute partial derivative approximations at x=a+dlx
!
   10   if (KSWX.NE.1) then
          UXXX = ( -3.0*U(1,J)+10.0*U(2,J)-12.0*U(3,J)+6.0*U(4,J)
     &                                          -U(5,J)) / TDLX3
          UXXXX = ( 2.0*U(1,J)-9.0*U(2,J)+16.0*U(3,J)-14.0*U(4,J)
     &                                 +6.0*U(5,J)-U(6,J)) / DLX4
          return
        endif 
!
! periodic at x=a+dlx
!
        UXXX = (-U(K-1,J)+2.0*U(1,J)-2.0*U(3,J)+U(4,J)) / (TDLX3)
        UXXXX = (U(K-1,J)-4.0*U(1,J)+6.0*U(2,J)-4.0*U(3,J)+U(4,J)) / 
     *      DLX4
        return
      endif
!
! compute partial derivative approximations on the interior
!
      UXXX = (-U(I-2,J)+2.0*U(I-1,J)-2.0*U(I+1,J)+U(I+2,J)) / TDLX3
      UXXXX = (U(I-2,J)-4.0*U(I-1,J)+6.0*U(I,J)-4.0*U(I+1,J)
     &                              +U(I+2,J)) / DLX4
      return
!
! compute partial derivative approximations at x=b-dlx
!
   20 if (KSWX.NE.1) then
        UXXX = (U(K-4,J)-6.0*U(K-3,J)+12.0*U(K-2,J)-10.0*U(K-1,J)
     &                              +3.0*U(K,J)) / TDLX3
        UXXXX = (-U(K-5,J)+6.0*U(K-4,J)-14.0*U(K-3,J)+16.0*U(K-2,J)
     &                             -9.0*U(K-1,J)+2.0*U(K,J)) / DLX4
        return
      endif 
!
! periodic at x=b-dlx
!
      UXXX = (-U(K-3,J)+2.0*U(K-2,J)-2.0*U(1,J)+U(2,J)) / TDLX3
      UXXXX = (U(K-3,J)-4.0*U(K-2,J)+6.0*U(K-1,J)-4.0*U(1,J)
     &                         +U(2,J)) / DLX4
      return
!
! compute partial derivative approximations at x=b
!
   30 UXXX = -( 3.0*U(K-4,J)-14.0*U(K-3,J)+24.0*U(K-2,J)
     &      -18.0*U(K-1,J)+5.0*U(K,J)) / TDLX3
      UXXXX = ( -2.0*U(K-5,J)+11.0*U(K-4,J)-24.0*U(K-3,J)
     &        +26.0*U(K-2,J)-14.0*U(K-1,J)+3.0*U(K,J)) / DLX4
      return
      end





      subroutine SEPDY (U,IDMN,I,J,UYYY,UYYYY)
!
! Compute second-order finite difference approximations to the third
! and fourth y partial derivatives of u at the (i,j) mesh point
!
      implicit none
      INT_TYPE IDMN, I,J
      REAL_TYPE U(IDMN,*), UYYY,UYYYY

#include "splp.h"

      if (J.LE.2.OR.J.GE.(L-1)) then
        if (J.NE.1) then
          if (J.EQ.  2) goto 10
          if (J.EQ.L-1) goto 20
          if (J.EQ.L  ) goto 30
        endif 
!
! compute partial derivative approximations at y=c
!
        If (KSWY.NE.1) then
          UYYY = (-5.0*U(I,1)+18.0*U(I,2)-24.0*U(I,3)+14.0*U(I,4)
     &                                       -3.0*U(I,5)) / TDLY3
          UYYYY = ( 3.0*U(I,1)-14.0*U(I,2)+26.0*U(I,3)-24.0*U(I,4)
     &                           +11.0*U(I,5)-2.0*U(I,6)) / DLY4
          return
        endif 
!
! periodic at x=a
!
        UYYY = (-U(I,L-2)+2.0*U(I,L-1)-2.0*U(I,2)+U(I,3)) / TDLY3
        UYYYY = (U(I,L-2)-4.0*U(I,L-1)+6.0*U(I,1)-4.0*U(I,2)+U(I,3))
     &      /DLY4
        return
!
! compute partial derivative approximations at y=c+dly
!
   10   If (KSWY.NE.1) then
          UYYY = (-3.0*U(I,1)+10.0*U(I,2)-12.0*U(I,3)+6.0*U(I,4)-
     *        U(I,5)) / TDLY3
          UYYYY = ( 2.0*U(I,1)-9.0*U(I,2)+16.0*U(I,3)-14.0*U(I,4)
     &                +6.0*U(I,5)-U(I,6) ) / DLY4
          return
        endif 
!
! periodic at y=c+dly
!
        UYYY = (-U(I,L-1)+2.0*U(I,1)-2.0*U(I,3)+U(I,4)) / TDLY3
        UYYYY = (U(I,L-1)-4.0*U(I,1)+6.0*U(I,2)-4.0*U(I,3)+U(I,4))
     &     / DLY4
        return
      endif 
!
! compute partial derivative approximations on the interior
!
      UYYY = (-U(I,J-2)+2.0*U(I,J-1)-2.0*U(I,J+1)+U(I,J+2)) / TDLY3
      UYYYY = (U(I,J-2)-4.0*U(I,J-1)+6.0*U(I,J)-4.0*U(I,J+1)
     &                                         +U(I,J+2)) / DLY4
      return
!
! compute partial derivative approximations at y=d-dly
!
   20 if (KSWY.NE.1) then
        UYYY = (U(I,L-4)-6.0*U(I,L-3)+12.0*U(I,L-2)-10.0*U(I,L-1)
     &                         +3.0*U(I,L)) / TDLY3
        UYYYY = (-U(I,L-5)+6.0*U(I,L-4)-14.0*U(I,L-3)+16.0*U(I,L-2)
     &                         -9.0*U(I,L-1)+2.0*U(I,L)) / DLY4
        return
      endif
!
! periodic at y=d-dly
!
      UYYY = (-U(I,L-3)+2.0*U(I,L-2)-2.0*U(I,1)+U(I,2)) / TDLY3
      UYYYY = (U(I,L-3)-4.0*U(I,L-2)+6.0*U(I,L-1)-4.0*U(I,1)
     &                                +U(I,2)) / DLY4
      return
!
! compute partial derivative approximations at y=d
!
   30 UYYY = -(3.0*U(I,L-4)-14.0*U(I,L-3)+24.0*U(I,L-2)
     &                 -18.0*U(I,L-1)+5.0*U(I,L)) / TDLY3
      UYYYY = (-2.0*U(I,L-5)+11.0*U(I,L-4)-24.0*U(I,L-3)
     &      +26.0*U(I,L-2)-14.0*U(I,L-1)+3.0*U(I,L)) / DLY4
      return
      end


!
! subroutines to compute the coefficients of the elliptic  equation
! solved to fill in the grid. They are passed (as arguments) to the
! elliptic solver sepeli
!
      subroutine COFX (XX,AFUN,BFUN,CFUN)
      implicit none
      REAL_TYPE XX,AFUN,BFUN,CFUN,  DXDI
      INT_TYPE I
#include "mexsepeli.h"
      I = NINT(XX)
      DXDI = 0.5 * (SXI(I+1)-SXI(I-1))
      AFUN = 1. / DXDI ** 2
      BFUN = (-SXI(I+1)+2.*SXI(I)-SXI(I-1)) / DXDI ** 3
      CFUN = 0.
      return
      end

      subroutine COFY (YY,DFUN,EFUN,FFUN)
      implicit none 
      REAL_TYPE YY,DFUN,EFUN,FFUN, DEDJ
      INT_TYPE J
#include "mexsepeli.h"
      J = NINT(YY)
      DEDJ = 0.5 * (SETA(J+1)-SETA(J-1))
      DFUN = 1. / DEDJ ** 2
      EFUN = (-SETA(J+1)+2.*SETA(J)-SETA(J-1)) / DEDJ ** 3
      FFUN = 0.
      return
      end




