      SUBROUTINE SEPELI (INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,
     1                   D,N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,GRHS,
     2                   USOL,IDMN,W,PERTRB,IERROR)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                        F I S H P A K                          *
C     *                                                               *
C     *                                                               *
C     *     A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF      *
C     *                                                               *
C     *      SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS        *
C     *                                                               *
C     *                  (VERSION 3.1 , OCTOBER 1980)                  *
C     *                                                               *
C     *                             BY                                *
C     *                                                               *
C     *        JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET         *
C     *                                                               *
C     *                             OF                                *
C     *                                                               *
C     *         THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH          *
C     *                                                               *
C     *                BOULDER, COLORADO  (80307)  U.S.A.             *
C     *                                                               *
C     *                   WHICH IS SPONSORED BY                       *
C     *                                                               *
C     *              THE NATIONAL SCIENCE FOUNDATION                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C
C DIMENSION OF           BDA(N+1), BDB(N+1), BDC(M+1), BDD(M+1),
C ARGUMENTS              USOL(IDMN,N+1), GRHS(IDMN,N+1),
C                        W (SEE ARGUMENT LIST)
C
C LATEST REVISION        SEPTEMBER 1996
C
C PURPOSE                SEPELI SOLVES FOR EITHER THE SECOND-ORDER
C                        FINITE DIFFERENCE APPROXIMATION OR A
C                        FOURTH-ORDER APPROXIMATION TO A SEPARABLE
C                        ELLIPTIC EQUATION
C
C                                    2    2
C                             AF(X)*D U/DX + BF(X)*DU/DX  + CF(X)*U +
C                                    2    2
C                             DF(Y)*D U/DY  + EF(Y)*DU/DY + FF(Y)*U
C
C                             = G(X,Y)
C
C                        ON A RECTANGLE (X GREATER THAN OR EQUAL TO A
C                        AND LESS THAN OR EQUAL TO B; Y GREATER THAN
C                        OR EQUAL TO C AND LESS THAN OR EQUAL TO D).
C                        ANY COMBINATION OF PERIODIC OR MIXED BOUNDARY
C                        CONDITIONS IS ALLOWED.
C
C PURPOSE                THE POSSIBLE BOUNDARY CONDITIONS ARE:
C                        IN THE X-DIRECTION:
C                         (0) PERIODIC, U(X+B-A,Y)=U(X,Y) FOR ALL Y,X
C                         (1) U(A,Y), U(B,Y) ARE SPECIFIED FOR ALL Y
C                         (2) U(A,Y), DU(B,Y)/DX+BETA*U(B,Y) ARE
C                             SPECIFIED FOR ALL Y
C                         (3) DU(A,Y)/DX+ALPHA*U(A,Y),DU(B,Y)/DX+
C                             BETA*U(B,Y) ARE SPECIFIED FOR ALL Y
C                         (4) DU(A,Y)/DX+ALPHA*U(A,Y),U(B,Y) ARE
C                             SPECIFIED FOR ALL Y
C
C                        IN THE Y-DIRECTION:
C                         (0) PERIODIC, U(X,Y+D-C)=U(X,Y) FOR ALL X,Y
C                         (1) U(X,C),U(X,D) ARE SPECIFIED FOR ALL X
C                         (2) U(X,C),DU(X,D)/DY+XNU*U(X,D) ARE SPECIFIED
C                             FOR ALL X
C                         (3) DU(X,C)/DY+GAMA*U(X,C),DU(X,D)/DY+
C                             XNU*U(X,D) ARE SPECIFIED FOR ALL X
C                         (4) DU(X,C)/DY+GAMA*U(X,C),U(X,D) ARE
C                             SPECIFIED FOR ALL X
C
C ARGUMENTS
C
C ON INPUT               INTL
C                          = 0 ON INITIAL ENTRY TO SEPELI OR IF ANY OF
C                              THE ARGUMENTS C, D, N, NBDCND, COFY ARE
C                              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 IS SOUGHT
C                          = 4 IF A FOURTH-ORDER APPROXIMATION 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 AND
C                          LESS THAN OR EQUAL TO B.  A MUST BE LESS THAN
C                          B.
C
C                        M
C                          THE NUMBER OF PANELS INTO WHICH THE INTERVAL
C                          [A,B] IS SUBDIVIDED.  HENCE, THERE WILL BE
C                          M+1 GRID POINTS IN THE X-DIRECTION GIVEN BY
C                          XI=A+(I-1)*DLX FOR I=1,2,...,M+1 WHERE
C                          DLX=(B-A)/M IS THE PANEL WIDTH.  M MUST BE
C                          LESS THAN IDMN AND GREATER THAN 5.
C
C                        MBDCND
C                          INDICATES THE TYPE OF BOUNDARY CONDITION AT
C                          X=A AND X=B
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 AND
C                              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 X=B
C                              ARE MIXED; I.E., DU(A,Y)/DX+ALPHA*U(A,Y)
C                              AND DU(B,Y)/DX+BETA*U(B,Y) ARE SPECIFIED
C                              FOR ALL Y
C                          = 4 IF THE BOUNDARY CONDITION AT X=A IS MIXED
C                              AND THE SOLUTION IS SPECIFIED AT X=B;
C                              I.E., DU(A,Y)/DX+ALPHA*U(A,Y) AND U(B,Y)
C                              ARE SPECIFIED FOR ALL Y
C
C                        BDA
C                          A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT
C                          SPECIFIES THE VALUES OF DU(A,Y)/DX+
C                          ALPHA*U(A,Y) AT X=A, WHEN MBDCND=3 OR 4.
C                               BDA(J) = DU(A,YJ)/DX+ALPHA*U(A,YJ);
C                               J=1,2,...,N+1
C                          WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A
C                          DUMMY PARAMETER.
C
C ON INPUT               ALPHA
C                          THE SCALAR MULTIPLYING THE SOLUTION IN CASE
C                          OF A MIXED BOUNDARY CONDITION AT X=A (SEE
C                          ARGUMENT BDA).  IF MBDCND ' 3,4 THEN ALPHA IS
C                          A DUMMY PARAMETER.
C
C                        BDB
C                          A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT
C                          SPECIFIES THE VALUES OF DU(B,Y)/DX+
C                          BETA*U(B,Y) AT X=B.  WHEN MBDCND=2 OR 3
C                               BDB(J) = DU(B,YJ)/DX+BETA*U(B,YJ);
C                               J=1,2,...,N+1
C                          WHEN MBDCND HAS ANY OTHER VALUE, BDB IS A
C                          DUMMY PARAMETER.
C
C                        BETA
C                          THE SCALAR MULTIPLYING THE SOLUTION IN CASE
C                          OF A MIXED BOUNDARY CONDITION AT X=B (SEE
C                          ARGUMENT BDB).  IF MBDCND'2,3 THEN BETA IS A
C                          DUMMY PARAMETER.
C
C                        C,D
C                          THE RANGE OF THE Y-INDEPENDENT VARIABLE;
C                          I.E., Y IS GREATER THAN OR EQUAL TO C AND
C                          LESS THAN OR EQUAL TO D.  C MUST BE LESS THAN
C                          D.
C
C                        N
C                          THE NUMBER OF PANELS INTO WHICH THE INTERVAL
C                          [C,D] IS SUBDIVIDED.  HENCE, THERE WILL BE
C                          N+1 GRID POINTS 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.  IN ADDITION,
C                          N MUST BE GREATER THAN 4.
C
C                        NBDCND
C                          INDICATES THE TYPES OF BOUNDARY CONDITIONS AT
C                          Y=C AND Y=D
C                          = 0 IF THE SOLUTION IS PERIODIC IN Y; I.E.,
C                              U(X,Y+D-C)=U(X,Y) FOR ALL X,Y
C                          = 1 IF THE SOLUTION IS SPECIFIED AT Y=C AND
C                              Y = D, I.E., U(X,C) AND U(X,D) ARE
C                              SPECIFIED FOR ALL X
C                          = 2 IF THE SOLUTION IS SPECIFIED AT Y=C AND
C                              THE BOUNDARY CONDITION IS MIXED AT Y=D;
C                              I.E., U(X,C) AND DU(X,D)/DY+XNU*U(X,D)
C                              ARE SPECIFIED FOR ALL X
C                          = 3 IF THE BOUNDARY CONDITIONS ARE MIXED AT
C                              Y=C AND Y=D; I.E., DU(X,C)/DY+GAMA*U(X,C)
C                              AND DU(X,D)/DY+XNU*U(X,D) ARE SPECIFIED
C                              FOR ALL X
C                          = 4 IF THE BOUNDARY CONDITION IS MIXED AT Y=C
C                              AND THE SOLUTION IS SPECIFIED AT Y=D;
C                              I.E. DU(X,C)/DY+GAMA*U(X,C) AND U(X,D)
C                              ARE SPECIFIED FOR ALL X
C
C                        BDC
C                          A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT
C                          SPECIFIES THE VALUE OF DU(X,C)/DY+GAMA*U(X,C)
C                          AT Y=C.  WHEN NBDCND=3 OR 4
C                             BDC(I) = DU(XI,C)/DY + GAMA*U(XI,C);
C                             I=1,2,...,M+1.
C                          WHEN NBDCND HAS ANY OTHER VALUE, BDC IS A
C                          DUMMY PARAMETER.
C
C                        GAMA
C                          THE SCALAR MULTIPLYING THE SOLUTION IN CASE
C                          OF A MIXED BOUNDARY CONDITION AT Y=C (SEE
C                          ARGUMENT BDC).  IF NBDCND'3,4 THEN GAMA IS A
C                          DUMMY PARAMETER.
C
C                        BDD
C                          A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT
C                          SPECIFIES THE VALUE OF DU(X,D)/DY +
C                          XNU*U(X,D) AT Y=C.  WHEN NBDCND=2 OR 3
C                            BDD(I) = DU(XI,D)/DY + XNU*U(XI,D);
C                            I=1,2,...,M+1.
C                          WHEN NBDCND HAS ANY OTHER VALUE, BDD IS A
C                          DUMMY PARAMETER.
C
C                        XNU
C                          THE SCALAR MULTIPLYING THE SOLUTION IN CASE
C                          OF A MIXED BOUNDARY CONDITION AT Y=D (SEE
C                          ARGUMENT BDD).  IF NBDCND'2 OR 3 THEN XNU IS
C                          A 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
C                          THE ELLIPTIC EQUATION AT X.
C
C                        COFY
C                          A USER-SUPPLIED SUBPROGRAM WITH
C                          PARAMETERS Y, DFUN, EFUN, FFUN WHICH
C                          RETURNS THE VALUES OF THE Y-DEPENDENT
C                          COEFFICIENTS DF(Y), EF(Y), FF(Y) IN
C                          THE ELLIPTIC EQUATION AT Y.
C
C                        NOTE:  COFX AND COFY MUST BE DECLARED EXTERNAL
C                        IN THE CALLING ROUTINE.  THE VALUES RETURNED IN
C                        AFUN AND DFUN MUST SATISFY AFUN*DFUN GREATER
C                        THAN 0 FOR A LESS THAN X LESS THAN B,
C                        C LESS THAN Y LESS THAN D (SEE IERROR=10).
C                        THE COEFFICIENTS PROVIDED MAY LEAD TO A MATRIX
C                        EQUATION WHICH IS NOT DIAGONALLY DOMINANT IN
C                        WHICH CASE SOLUTION MAY FAIL (SEE IERROR=4).
C
C                        GRHS
C                          A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE
C                          VALUES OF THE RIGHT-HAND SIDE OF THE ELLIPTIC
C                          EQUATION; I.E., GRHS(I,J)=G(XI,YI), FOR
C                          I=2,...,M; J=2,...,N.  AT THE BOUNDARIES,
C                          GRHS IS 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 QUANTITES 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 QUANTITES ARE NOT USED IN
C                          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 OF
C                          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 BE
C                          CHOSEN.  FOR EXAMPLE, IF MBDCND=2 AND
C                          NBDCND=4, THEN U(A,C), U(A,D), U(B,D) MUST BE
C                          CHOSEN AT THE CORNERS IN ADDITION 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 TO
C                          SPECIFY THE VARIABLE DIMENSION OF GRHS AND
C                          USOL.  IDMN MUST BE AT LEAST 7 AND GREATER
C                          THAN OR EQUAL TO M+1.
C
C                        W
C                          A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED
C                          BY THE USER FOR WORK SPACE.  LET
C                          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 IN
C                          THE CALLING ROUTINE MUST BE SET IN W(1) (SEE
C                          IERROR=11).
C
C ON OUTPUT              USOL
C                          CONTAINS THE APPROXIMATE SOLUTION TO THE
C                          ELLIPTIC EQUATION.  USOL(I,J) IS THE
C                          APPROXIMATION TO U(XI,YJ) FOR I=1,2...,M+1
C                          AND J=1,2,...,N+1.  THE APPROXIMATION HAS
C                          ERROR O(DLX**2+DLY**2) IF CALLED WITH
C                          IORDER=2 AND O(DLX**4+DLY**4) IF CALLED WITH
C                          IORDER=4.
C
C                        W
C                          CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE
C                          DESTROYED IF SEPELI IS CALLED AGAIN WITH
C                          INTL=1.  IN ADDITION W(1) CONTAINS THE EXACT
C                          MINIMAL LENGTH (IN FLOATING POINT) REQUIRED
C                          FOR THE WORK SPACE (SEE IERROR=11).
C
C                        PERTRB
C                          IF A COMBINATION OF PERIODIC OR DERIVATIVE
C                          BOUNDARY CONDITIONS (I.E., ALPHA=BETA=0 IF
C                          MBDCND=3; GAMA=XNU=0 IF NBDCND=3) IS
C                          SPECIFIED AND IF THE COEFFICIENTS OF U(X,Y)
C                          IN THE SEPARABLE ELLIPTIC EQUATION ARE ZERO
C                          (I.E., CF(X)=0 FOR X GREATER THAN OR EQUAL TO
C                          A AND LESS THAN OR EQUAL TO B; FF(Y)=0 FOR
C                          Y GREATER THAN OR EQUAL TO C AND LESS THAN
C                          OR EQUAL TO D) THEN A SOLUTION MAY NOT EXIST.
C                          PERTRB IS A CONSTANT CALCULATED AND
C                          SUBTRACTED FROM THE RIGHT-HAND SIDE OF THE
C                          MATRIX EQUATIONS GENERATED BY SEPELI WHICH
C                          INSURES THAT A SOLUTION EXISTS.  SEPELI THEN
C                          COMPUTES THIS SOLUTION WHICH IS A WEIGHTED
C                          MINIMAL LEAST SQUARES SOLUTION TO THE
C                          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 (SEE DISCUSSION OF
C                              IDMN)
C                          = 6 IF M IS TOO SMALL OR TOO LARGE (SEE
C                              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 FOR
C                               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, THE
C                          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 MAY
C                          FAIL.  THIS CANNOT HAPPEN IN THE LIMIT AS
C                          DLX, DLY APPROACH ZERO.  HENCE, THE CONDITION
C                          MAY BE REMEDIED BY TAKING LARGER VALUES FOR M
C                          OR N.
C
C ENTRY POINTS           SEPELI, SPELIP, CHKPRM, CHKSNG, ORTHOG, MINSOL,
C                        TRISP, DEFER, DX, DY, BLKTRI, BLKTR1, INDXB,
C                        INDXA, INDXC, PROD, PRODP, CPROD, CPRODP,
C                        PPADD, PSGF, BSRH, PPSGF, PPSPF, COMPB,
C                        TRUN1, STOR1, TQLRAT
C
C SPECIAL CONDITIONS     NONE
C
C COMMON BLOCKS          SPLP, CBLKT, VALU1
C
C I/O                    NONE
C
C PRECISION              SINGLE
C
C SPECIALIST             JOHN C. ADAMS, NCAR, BOULDER, COLORADO  80307
C
C LANGUAGE               FORTRAN
C
C HISTORY                DEVELOPED AT NCAR DURING 1975-76.
C
C ALGORITHM              SEPELI AUTOMATICALLY DISCRETIZES THE SEPARABLE
C                        ELLIPTIC EQUATION WHICH IS THEN SOLVED BY A
C                        GENERALIZED CYCLIC REDUCTION ALGORITHM IN THE
C                        SUBROUTINE, BLKTRI.  THE FOURTH-ORDER SOLUTION
C                        IS OBTAINED USING 'DEFERRED CORRECTIONS' WHICH
C                        IS DESCRIBED AND REFERENCED IN SECTIONS,
C                        REFERENCES AND METHOD.
C
C SPACE REQUIRED         14654 (OCTAL) = 6572 (DECIMAL)
C
C ACCURACY AND TIMING    THE FOLLOWING COMPUTATIONAL RESULTS WERE
C                        OBTAINED BY SOLVING THE SAMPLE PROBLEM AT THE
C                        END OF THIS WRITE-UP ON THE CONTROL DATA 7600.
C                        THE OP COUNT IS PROPORTIONAL TO M*N*LOG2(N).
C                        IN CONTRAST TO THE OTHER ROUTINES IN THIS
C                        CHAPTER, ACCURACY IS TESTED BY COMPUTING AND
C                        TABULATING SECOND- AND FOURTH-ORDER
C                        DISCRETIZATION ERRORS.  BELOW IS A TABLE
C                        CONTAINING COMPUTATIONAL RESULTS.  THE TIMES
C                        GIVEN DO NOT INCLUDE INITIALIZATION (I.E.,
C                        TIMES ARE FOR INTL=1).  NOTE THAT THE
C                        FOURTH-ORDER ACCURACY IS NOT REALIZED UNTIL THE
C                        MESH IS SUFFICIENTLY REFINED.
C
C              SECOND-ORDER    FOURTH-ORDER   SECOND-ORDER  FOURTH-ORDER
C    M    N   EXECUTION TIME  EXECUTION TIME    ERROR         ERROR
C               (M SEC)         (M SEC)
C     6    6         6              14          6.8E-1        1.2E0
C    14   14        23              58          1.4E-1        1.8E-1
C    30   30       100             247          3.2E-2        9.7E-3
C    62   62       445           1,091          7.5E-3        3.0E-4
C   126  126     2,002           4,772          1.8E-3        3.5E-6
C
C PORTABILITY            THERE ARE NO MAHCINE-DEPENDENT CONSTANTS.
C
C REQUIRED RESIDENT      SQRTF, ABS, CABS, LOGF
C ROUTINES
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
C
C
C
      DIMENSION       GRHS(IDMN,1)           ,USOL(IDMN,1)
      DIMENSION       BDA(1)     ,BDB(1)     ,BDC(1)     ,BDD(1)     ,
     1                W(1)
      EXTERNAL        COFX       ,COFY
C
C     CHECK INPUT PARAMETERS
C
      CALL CHKPRM (INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY,
     1             IDMN,IERROR)
      IF (IERROR .NE. 0) RETURN
C
C     COMPUTE MINIMUM WORK SPACE AND CHECK WORK SPACE LENGTH INPUT
C
      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
C
C     SET WORK SPACE INDICES
C
      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,
     1             NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,W(I1),W(I2),W(I3),
     2             W(I4),W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),
     3             W(I12),GRHS,USOL,IDMN,W(I13),PERTRB,IERROR)
      RETURN
      END
