C     ALGORITHM 637 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 4,
C     DEC., 1985, P. 413-415.
C//////////////////////////////////////////////////////////////
C//////////////////  ALGORITHM  GENCOL ///////////////////////
C/////////////////////////////////////////////////////////////
C>>>>>>>>>>>>>>>>>> LOGICAL FILE 1 : DRIVER <<<<<<<<<<<<<<<<<<
C/////////////////////////////////////////////////////////////
C     PROGRAM DRIVE1
C
C
C  ---- T E S T   D R I V E R   F O R   G E N C O L -----
C
C
C MODIFIABLE PARAMETERS FOR DIMENSION SETTING.
C  THE FOLLOWING DIMENSIONS LIMIT THE SIZE OF THE PROBLEM THAT
C  CAN BE SOLVED.
C
C  NGXMAX - MAXIMUM NUMBER OF X GRID LINES
C  NGYMAX - MAXIMUM NUMBER OF Y GRID LINES
C  NTBXMX - MAXIMUM NUMBER OF X GRID LINES IN TABLED OUTPUT ON A
C           GRID OTHER THAN THE DISCRETIZATION GRID
C  NTBYMX - MAXIMUM NUMBER OF Y GRID LINES IN TABLED OUTPUT ON A
C           GRID OTHER THAN THE DISCRETIZATION GRID
C  NOUTMX - MAXIMUM NUMBER OF OUTPUT REQUESTS
C  NBPTMX - MAXIMUM NUMBER OF BOUNDARY POINTS (CORNERS OF DOMAIN
C           AND INTERSECTIONS OF DOMAIN BOUNARY WITH GRID)
C  NPMAX  - MAXIMUM NUMBER OF PIECES OF DOMAIN BOUNDARY
C  NEQMAX - MAXIMUM NUMBER OF EQUATIONS (= 4*NGXMAX*NGYMAX) )
C  NCOEMX - MAXIMUM NUMBER OF COEFFICIENTS
C  NWKMAX - MAXIMUM AMOUNT OF WORKSPACE (= 2+(12*NGYMAX + 23)*NEQMAX ) )
C
C     TABX        NTBXMX             BPTYPE      NBPTMX
C     TABY        NTBYMX             BNEIGH      NBPTMX
C     IDCOEF      NEQMAX,NCOEMX      BGRID       NBPTMX
C     OUTFNC      NOUTMX             GRIDX       NGXMAX
C     OUTTYP      NOUTMX             GRIDY       NGYMAX
C     GTYPE       NGXMAX,NGYMAX      BRANGE      2,NPMAX
C     XBOUND      NBPTMX             NODELM      NEQMAX
C
C
      REAL COEF(324,17),UNKVCT(324),WRKSP(42446),TABX(10),TABY(10)
      INTEGER IDCOEF(324,17),OUTFNC(4),OUTTYP(4),GTYPE(9,9),GIVOPT,
     .        GIVOPZ
      LOGICAL DOMFLG,PLOTIT,USECRN,PLOTIZ,USECRZ
C
C
C  MODIFICATIONS TO CHANGE PROBLEM
C
C  THE DRIVER IS SET TO RUN EXAMPLE 1
C
C  TO RUN EXAMPLE 2, CHANGE NOUT TO 4 AND MAKE THE MODIFICATIONS
C  TO BRANGE AND BY AS INDICATED BELOW
C
C
C********************* NON-STANDARD COMMON BLOCK USE ****************
C       MANY LABELED COMMON BLOCKS ARE DIMENSIONED WITH LENGTH 1    *
C       IN THIS ALGORITHM.  THE CORRECT LENGTHS MUST BE SET IN THE  *
C       DRIVER(MAIN) PROGRAM.  THIS IS NOT STANDARD FORTRAN 77 BUT  *
C       IT WORKS ON ALL SYSTEMS FORTRAN KNOWN TO US.  IT ALLOWS ONE *
C       TO USE THE PROGRAMS REPEATEDLY WITHOUT RECOMPILING THEM.    *
C    EXAMPLE.                                                       *
C      COMMON / GRIDXZ / GRIDX(1)                                   *
C    IS USED INSTEAD OF                                             *
C      COMMON / GRIDXZ / GRIDX(9)                                   *
C       IF YOUR COMPILER ENFORCES THE STANDARD THEN YOU MUST CHANGE *
C       THESE COMMON BLOCK LENGTHS AND RECOMPILE EACH TIME          *
C********************* NON-STANDARD COMMON BLOCK USE ****************
C
      COMMON /PROBR/AX,BX,AY,BY
      COMMON /XBOUZZ/XBOUND(50)
      COMMON /YBOUZZ/YBOUND(50)
      COMMON /BPARZZ/BPARAM(50)
      COMMON /PIECZZ/PIECE(50)
      COMMON /BPTYZZ/BPTYPE(50)
      COMMON /BNEIZZ/BNEIGH(50)
      COMMON /BGRIZZ/BGRID(50)
      COMMON /GRIDXZ/GRIDX(9)
      COMMON /GRIDYZ/GRIDY(9)
      COMMON /BRANZZ/BRANGE(2,4)
      COMMON /COLOPT/BCP1Z,BCP2Z,DSCARZ,PTSIZZ,GIVOPZ,PLOTIZ,USECRZ
      COMMON /COLNUM/NODELM(324)
C
C -------  MODIFIABLE DATA STATEMENT -------
C
C IGRID - INTEGER NOT 0 OR 1
C         IF POSITIVE, THEN UNIFORM IGRID X IGRID
C         IF NEGATIVE, THEN VALUES LISTED IN SUBROUTINE SETGRD
C
C OTHER PARAMETERS ARE DESCRIBED IN SUBROUTINE GENCOL
C
      DATA IGRID,BCP1,BCP2,DSCARE,GIVOPT,PLOTIT,USECRN,LEVEL/5,0.,0.,
     .     .05,1,.FALSE.,.FALSE.,1/
C
      CALL KILLIT(144)
      DOMFLG = .TRUE.
      PTSIZE = 7.
C
      MXNEQ = 324
      MXNCOE = 17
      NBDIM = 50
      MOUTPT = 6
C
      BCP1Z = BCP1
      BCP2Z = BCP2
      DSCARZ = DSCARE
      PTSIZZ = PTSIZE
      GIVOPZ = GIVOPT
      PLOTIZ = PLOTIT
      USECRZ = USECRN
C
C SET BOUNDARY DATA STRUCTURES
C
C   ----THIS CODE IS FOR EXAMPLE 1, THE UNIT SQUARE
C
      NBOUND = 4
      DO 10 I = 1,4
         BRANGE(1,I) = 0.
         BRANGE(2,I) = 1.
   10 CONTINUE
      AX = 0.
      BX = 1.
      AY = 0.
      BY = 1.
C  --------------------------------------------------
C
C    FOR EXAMPLE 2'S NONRECTANGULAR DOMAIN, INCLUDE
C       PI = 4.*ATAN(1.)
C       BRANGE(2,4) = PI/2.
C       BY = .5
C
C  --------------------------------------------------
CCCC    PI = 4.*ATAN(1.)
CCCC    BRANGE(2,4) = PI/2.
CCCC    BY = .5
      NPDIM = NBOUND
C
C SET GRID
C
      CALL SETGRD(GRIDX,GRIDY,NGRIDX,NGRIDY,AX,BX,AY,BY,IGRID)
C
C SET OUTPUT TYPES
C  SHOULD BE 2 FOR EXAMPLE 1, 4 FOR EXAMPLE 2 , 3
C
      NOUT = 4
C TABLE APPROXIMATE SOLUTION
      OUTFNC(1) = 1
      OUTTYP(1) = 3
C TABLE RESIDUAL ON 10 X 10 GRID
      OUTFNC(2) = 3
      OUTTYP(2) = 4
      NTABX = 10
      NTABY = 10
      DO 20 I = 1,10
         AIM1 = I - 1
         TABX(I) = AX + (BX-AX)*AIM1/9.
         TABY(I) = AY + (BY-AY)*AIM1/9.
   20 CONTINUE
C MAXIMUM ERROR (MEANINGLESS FOR EXAMPLE 1 WHERE TRUE IS UNKNOWN)
      OUTFNC(3) = 2
      OUTTYP(3) = 1
C MAXIMUM ERROR ON 10 X 10 GRID (MEANINGLESS FOR EXAMPLE 1)
      OUTFNC(4) = 2
      OUTTYP(4) = 2
C
      CALL GENCOL(GTYPE,NGRIDX,NGRIDY,NBOUND,NBNDPT,NPDIM,NBDIM,DOMFLG,
     .            AX,BX,AY,BY,NGRIDX,NGRIDY,COEF,IDCOEF,UNKVCT,MXNCOE,
     .            MXNEQ,OUTFNC,OUTTYP,NOUT,TABX,TABY,NTABX,NTABY,LEVEL,
     .            MOUTPT,WRKSP)
      STOP
      END
      SUBROUTINE KILLIT(N)
C     CDC-NOS DEPENDENT ROUTINE TO SUPPRESS EXP() UNDERFLOW MESSAGES.
      DIMENSION IS(6)
      DATA (IS(K),K=1,6)/-0,-0,-0,0,-0,-0/
      CALL SYSTEMC(N,IS)
      RETURN
      END
      SUBROUTINE REGION(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,BCOORD,
     .                  SCLOCK,SARC,SLEVEL,GTYPE,XBOUND,YBOUND,PIECE,
     .                  BPTYPE,BNEIGH,BGRID,BPARAM,NBDIM,NBPTS,FAIL)
C
C    DOMAIN PROCESSOR.  WRITTEN BY JOHN R. RICE, PURDUE UNIVERSITY.
C
C    THE PURPOSE OF THIS ALGORITHM IS TO RELATE A GENERAL 2 DIMENSIONAL
C    DOMAIN TO A RECTANGULAR GRID LAID OVER IT.  THE DOMAIN BOUNDARY IS
C    GIVEN IN PARAMETRIC FORM AND THE PROGRAM PRODUCES ARRAYS WHICH
C    DESCRIBE THE RELATIONSHIP BETWEEN THE TWO OBJECTS.  THESE ARRAYS
C    CONTAIN INFORMATION POINTING IN BOTH DIRECTIONS AND WHICH SHOULD
C    FACILITATE APPLICATIONS( E.G. NUMERICAL QUADRATURE, SURFACE
C    FITTING, DISCRETIZING PDE'S ) WHICH INVOLVE SUCH DOMAINS.
C
C
C  ------------------- INPUT INFORMATION FOR DOMAIN PROCESSOR --------
C
C 1** OUTPUT LEVEL CONTROL *******
C         SLEVEL = LEVEL = OUTPUT CONTROL SETTING
C                          DETAILS GIVEN IN SUBROUTINE DOMAIN
C
C 2** GRID DEFINITION ************
C
C         NGDIMX,NGDIMY  = NUMBER OF GRID LINES IN X AND Y COORDINATES
C                          NGRIDX = NGDIMX, NGRIDY = NGDIMY
C         XGRID(IX),YGRID(JY)  FOR  IX = 1 TO NGRIDX, JY = 1 TO NGRIDY
C                        = ARRAYS CONTAINING THE GRID LINE COORDINATES
C
C 3** BOUNDARY DEFINITION ********
C
C         NPDIM  = NBOUND = NUMBER OF BOUNDARY PIECES
C         NBDIM  = ARRAY DIMENSION FOR BOUNDARY POINTS
C         BCOORD = PARAMETERIZED DEFINITION OF THE BOUNDARY.
C                  BCOORD(P,X,Y,IPIECE) GIVES THE X,Y VALUES OF THE
C                  POINT ON PIECE IPIECE WITH PARAMETER VALUE = P.
C                  SEE NOTE BELOW FOR DISCRETE BOUNDARY DEFINITION
C
C         BRANGE(2,I) = FIRST AND LAST VALUES OF PARAMETERS DEFINING
C                       THE I-TH BOUNDARY PIECE
C         SCLOCK = CLOCK = SWITCH TO SPECIFY BOUNDARY ORIENTATION
C                  .TRUE.  MEANS BOUNDARY IS CLOCKWISE
C                  .FALSE. MEANS BOUNDARY IS COUNTER-CLOCKWISE
C         SARC   = ARC = .TRUE. MEANS DOMAIN IS AN ARC WITH NO INTERIOR
C
C  ---------------------------------------------------------------------
C  ---------  OUTPUT INFORMATION OF THE DOMAIN PROCESSOR IS IN TWO PARTS
C               ( THERE IS ALSO A FAILURE FLAG = FAIL )
C
C  PART 1 ********** GRID POINT TYPES ************
C
C    GTYPE(IX,JY) FOR  IX = 1 TO NGRIDX, JY = 1 TO NGRIDY
C      THE VALUES IN THIS ARRAY GIVE THE TYPE OF THE GRID POINTS AND
C      INFORMATION ABOUT THEIR RELATION TO THE BOUNDARY. DETAILS AND
C      EXAMPLES ARE GIVEN IN THE SUBROUTINE DOMAIN.
C      POSSIBLE VALUE ARE( ASSUMING PACKING FACTOR IPACKB = 1000)
C
C    = INTEGER OVER 1000
C      GRID POINT IS NEXT TO THE BOUNDARY AND THE GTYPE VALUE IS
C      GTYPE = K + 1000*J   WHERE
C            K IS THE INDEX OF THE LOWEST NUMBERED BOUNDARY NEIGHBOR
C            J = FOUR BITS TO NOTE LOCATION OF BOUNDARY POINTS
C                 0001 - BOUNDARY NEIGHBOR TO NORTH  (NOON)
C                 0010 - BOUNDARY NEIGHBOR TO EAST   (3 O'CLOCK)
C                 0100 - BOUNDARY NEIGHBOR TO SOUTH  (6 O'CLOCK)
C                 1000 - BOUNDARY NEIGHBOR TO WEST   (9 O'CLOCK)
C
C    = 999
C      GRID POINT IS INTERIOR TO REGION AND NOT CLOSE TO BOUNDARY
C
C    = INTEGER LESS THAN 999
C      GRID POINT IS ALSO BOUNDARY PT., GTYPE IS ITS INDEX
C
C    = 0
C      GRID POINT IS EXTERIOR FAR FROM THE BOUNDARY
C
C    = NEGATIVE INTEGER
C      GRID POINT IS EXTERIOR NEXT TO THE BOUNDARY, ITS LOCATION
C      RELATIVE TO THE BOUNDARY IS ENCODED AS FOR INTERIOR POINTS
C      NEAR THE BOUNDARY
C
C
C  PART 2 ********* BOUNDARY/GRID INTERSECTIONS **********
C
C         NBNDPT = NUMBER OF BOUNDARY POINTS ACTUALLY FOUND
C                  BOUNDARY POINT INDEX RANGE = 1 TO NBNDPT+1
C                  THE FIRST BOUNDARY POINT IS ADDED TO THE LIST TO
C                  MAKE IT CIRCULAR, HENCE THE +1 IN THE INDEX RANGE
C         XBOUND(I),YBOUND(I) = COORDINATES OF I-TH BOUNDARY POINT
C         BPARAM(I) = PARAMETER VALUE OF I-TH BOUNDARY POINT
C         PIECE(I)  = INDEX OF BOUNDARY PIECE TO WHICH PT. BELONGS
C                        SMALLEST NUMBER FOR CORNER POINTS
C         BPTYPE(I) = TYPE OF BOUNDARY POINT
C                   = HORZ  IF POINT IN ON A Y GRID LINE
C                   = VERT                   X
C                   = BOTH  IF POINT IS ALSO A GRID POINT
C                   = INTE  IF POINT IS NOT ON A GRID LINE
C                           HAPPENS ONLY FOR CORNERS OF THE BOUNDARY
C                   = JUMP  IF POINT PRECEEDS A HOLE
C         BNEIGH(I) = POINTER TO THE INTERIOR POINTS FROM THE I-TH
C                     BOUNDARY POINT. SAME SCHEME IS USED TO ENCODE
C                     DIRECTIONS AS FOR THE J PART OF GTYPE ABOVE
C         BGRID(I) = IX + IPACKB*JY WHEN PT. I IS IN GRID SQUARE IX,JY
C
C  *********************************************************************
C  **************  DOMAIN PROCESSOR ALGORITHM STRUCTURE  ***************
C
C     THERE ARE 2 BASIC PARTS TO THE ALGORITHM, PROCESSING THE BOUNDARY
C     AND TYPING THE GRID POINTS. AN OUTLINE FOLLOWS
C
C    MARK ALL GRID POINTS AS EXTERIOR
C
C*******   PART 1 *** PROCESS BOUNDARY
C
C    LOCATE FIRST BOUNDARY POINT
C
C    DO LOOP OVER BOUNDARY PIECES
C     1 DO WHILE NOT AT END OF PIECE
C     1  1 LOCATE NEXT INTERSECTION WITH GRID LINE( SUBROUTINE BWALK )
C     1  1 DETERMINE TYPE FOR INTERSECTION POINT FOUND
C     1  1 CHECK FOR CHANGING BOUNDARY PIECE(SUBROUTINE CHANGE DOES IT)
C     1  1 CHECK CONTINUITY OF BOUNDARY
C     1  +---
C     +------
C
C    CHECK CLOSING OF BOUNDARY
C    FINISH UP FIRST POINT
C
C******* PART 2 *** TYPE GRID POINTS AND MARK INTERIOR POINTS
C
C    PASS OVER GRID MARKING THE GRID POINT TYPES(SUBROUTINE GVALUS )
C    MARK INTERIOR POINTS( SUBROUTINE FILL )
C          LOCATE FIRST INTERIOR POINT
C          MARK THE INTERIOR POINTS
C    SET POINTERS FROM BOUNDARY TO INTERIOR( SUBROUTINE NEIGH )
C
C
C   **** NOTE. IF THE BOUNDARY IS DESCRIBED DISCRETELY BY THE SET OF
C   BOUNDARY-GRID INTERSECTIONS, THEN THE FOLLOWING FUNCTION PROGRAM
C   ALLOWS THE DOMAIN PROCESSOR TO BE APPLIED.  IT ESSENTIALLY PROVIDES
C   LINEAR INTERPOLATION BETWEEN THE DISCRETE POINTS
C
C    *     FUNCTION BCOORD(P,X,Y,IPIECE)
C    *     REAL XBOUND(***),YBOUND(***)
C    *     LOGICAL STARTD
C    *     DATA STARTD/.FALSE./
C    *     IF( STARTD )             GO TO 20
C    *           OBTAIN DATA FOR DISCRETE POINTS THROUGH A COMMON BLOCK
C    *           OR BY READING FROM A FILE.  THE READ CASE IS SHOWN.
C    *         READ(INFILE,10) NBNDPT,(XBOUND(I),YBOUND(I),I=1,NBNDPT-1)
C    *   10    FORMAT( ****** )
C    *         XBOUND(NBNDPT) = XBOUND(1)
C    *         YBOUND(NBNDPT) = YBOUND(1)
C    *
C    *            OBTAIN COORDINATES FROM ARRAYS OF DISCRETE DATA
C    *   20 IP = P
C    *      X = (P-IP)*XBOUND(IP+1) + (1.+P-IP)*XBOUND(IP)
C    *      Y = (P-IP)*YBOUND(IP+1) + (1.+P-IP)*YBOUND(IP)
C    *      RETURN
C    *      END
C
C  *****************************************************************
C  ***** REGION IS AN INTERFACE SUBROUTINE FOR THE DOMAIN PROCESSING
C            MAIN SUBPROGRAM = DOMAIN.
C      REGION MUST BE THE FIRST ROUTINE CALLED FOR PROCESSING
C            MULTIPLE REGIONS( E.G. WITH HOLES)
C      THE FUNCTIONS OF THIS INTERFACE PROGRAM ARE AS FOLLOWS.
C
C      1.  INITIALIZE OR COMPUTE VARIOUS CONSTANTS AND PLACE THEM IN
C          COMMON BLOCKS.  THESE ARE
C          A. I/O CONSTANT = MOUTPT
C          B. NUMERICAL CONTROL CONSTANTS = EPSGRD,EPSTAN
C          C. 5 CHARACTER STRING CONSTANTS
C          D. COUNTERS = N1BND,N1BDPT = 1
C          E. INTERIOR MARKER = NSIDE
C          F. PACKING FACTOR = IPACKB( SHOULD BE NSIDE+1 )
C
C      2.  CHECK THAT THE GRID IS PROPERLY DEFINED.
C
C      3.  MOVE VARIABLES FROM SUBROUTINE ARGUMENT LIST INTO COMMON.
C          A. DIMENSION SPECIFICATIONS = NGDIMX,NGDIMY,NBOUND
C          B. CONTROL VARIABLES = LEVEL, CLOCKW, ARC
C
C      4.  CALL DOMAIN.
C
C      5.  SET FAILURE FLAG FAIL, FINAL NUMBER NBPTS OF BOUNDARY POINTS
C
C
      REAL XGRID(NGDIMX),YGRID(NGDIMY),BRANGE(2,NPDIM)
      LOGICAL SCLOCK,SARC,FAIL
      INTEGER SLEVEL
      DIMENSION GTYPE(NGDIMX,NGDIMY),XBOUND(NBDIM),YBOUND(NBDIM),
     .          BPTYPE(NBDIM),BNEIGH(NBDIM),BGRID(NBDIM),BPARAM(NBDIM),
     .          PIECE(NBDIM)
C
C     ******** IF HOLE IS ALSO USED THEN THE CALLING PROGRAM  *****
C     ******** MUST HAVE THE FOLLOWING COMMON BLOCKS IN ORDER *****
C     ******** TO CONFORM TO FORTRAN STANDARDS                *****
C
      EXTERNAL BCOORD

      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN,EPSGRQ
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /NUMCON/EPSGRQ
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
C  CONSTANT CHARACTER STRINGS FOR NICE PRINT OUT
      CHARACTER *4 IQHORZ,IQVERT,IQBOTH,IQINTE,IQJUMP
      DATA IQHORZ,IQVERT,IQBOTH,IQINTE,IQJUMP/'HORZ','VERT','BOTH',
     .     'INTE','JUMP'/
C
C  INITIALIZE COMMON BLOCKS FROM SUBROUTINE ARGUMENTS
      NGRIDX = NGDIMX
      NGRIDY = NGDIMY
      NBOUND = NPDIM
      LEVEL = SLEVEL
      CLOCKW = SCLOCK
      ARC = SARC
C
C  INITIALIZE CONSTANT CHARACTER STRINGS FROM DATA STATEMENTS
      HORZ = IQHORZ
      VERT = IQVERT
      BOTH = IQBOTH
      INTER = IQINTE
      JUMP = IQJUMP
C                 SET UNIT NUMBER FOR OUTPUT
      MOUTPT = 6
C  PRINT GREETINGS
C
      IF (LEVEL.GT.0) WRITE (MOUTPT,9001)
C
C                 SET MARKER FOR INTERIOR POINTS OF REGION
C                 MUST EXCEED MAX NUMBER OF BOUNDARY POINTS
      NSIDE = 999
      IPACKB = NSIDE + 1
C                 SET LIMIT ON THE QUEUE LENGTH IN ROUTINE FILL
C                 QLIMIT = 250 SHOULD ALLOW WELL OVER 10,000
C                 INTERIOR POINTS IN THE DOMAIN
      QLIMIT = 250
C                 MARK THAT WE ARE NOT IN HOLE OR ARC
      ARC = .FALSE.
      INHOLE = .FALSE.
C                 INITIALIZE BOUNDARY LIST COUNTERS
      N1BND = 1
      N1BDPT = 1
C
C                 SET GEOMETRY TOLERANCE PARAMETER EPSGRD
C     ****************************************************************
C     ******                                                    ******
C     ******  THIS IS A MACHINE AND PROBLEM DEPENDENT CONSTANT  ******
C     ******                                                    ******
C     ****************************************************************
C
C     EPSGRD SHOULD BE LARGE ENOUGH TO INSULATE THE COMPUTATIONS FROM
C     MACHINE ROUND-OFF.  ALL POINTS AND LINES WITHIN EPSGRD OF ONE
C     ANOTHER ARE ASSUMED TO BE EQUAL.  IT IS PROBABLY SAFE TO TAKE
C     EPSGRD AS 50 UNITS IN THE LAST PLACE, IT MUST BE AT LEAST 20
C     UNITS IN THE LAST PLACE. THE CONVERGENCE TEST IN SECANT IS
C     .2*EPSGRD
C
C     ******  THIS PARAMETER IS RELATIVE TO THE DOMAIN SIZE  ******
C
C     EPSGRD SHOULD BE SMALL ENOUGH SO THAT THE ACCURACY IN THE PROBLEM
C     SOLUTION IS NOT AFFECTED BY AN UNCERTAINTY IN THE GEOMETRY OF
C     EPSGRD.
C
C        EPSGRD = 1.E-8 IS APPROPRIATE FOR MOST LONG WORD LENGTH
C                       MACHINES AND PROBLEMS
C                 2.E-5 IS APPROPRIATE( BUT NOT FAIL-SAFE ) FOR
C                       32 BIT MACHINES.
      XWIDTH = XGRID(NGRIDX) - XGRID(1)
      YWIDTH = YGRID(NGRIDY) - YGRID(1)
      EPSGRD = 2.E-5*AMAX1(XWIDTH,YWIDTH)
C     PUT EPSGRD IN COMMON /NUMCON/ FOR DOMAIN PROCESSING.
      EPSGRQ = EPSGRD
C
C         CHECK THAT THE GRIDS ARE PROPERLY DEFINED
      IF (NGRIDX.GE.2) GO TO 20
C               FATAL ERROR, GRID TOO SMALL
   10 CONTINUE
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9011)
      FATAL = .TRUE.
      RETURN

   20 IF (NGRIDY.LE.1) GO TO 10
C
C          CHECK THAT THE GRID INCREASES IN BOTH DIRECTIONS
C          FIND THE MINIMUM GRID WIDTHS
      XGMIN = XWIDTH
      DO 30 I = 2,NGRIDX
         XGMIN = AMIN1(XGMIN,XGRID(I)-XGRID(I-1))
   30 CONTINUE
      YGMIN = YWIDTH
      DO 40 I = 2,NGRIDY
         YGMIN = AMIN1(YGMIN,YGRID(I)-YGRID(I-1))
   40 CONTINUE
      IF (AMIN1(XGMIN,YGMIN).GT.EPSGRD) GO TO 50
C
C        FATAL ERROR, HAVE ZERO OR NEGATIVE GRID WIDTH
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9021)
      FATAL = .TRUE.
      RETURN
C
C           SET PARAMETER FOR TANGENCY TEST
C           SEE SUBROUTINE CHKTAN FOR DETAILS
   50 EPSTAN = .1E0*AMIN1(XGMIN,YGMIN,15.E0*SQRT(EPSGRD))
C
C  CALL DOMAIN PROCESSOR
      CALL DOMAIN(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,BCOORD,GTYPE,
     .            XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH,BGRID,BPARAM,NBDIM)
C
C           SET FAILURE FLAG AND FINAL NUMBER OF BOUNDARY POINTS
      FAIL = FATAL
C           THESE NUMBERS WILL BE USED BY HOLE IF IT IS CALLED
      NBPTS = NBNDPT
      N1BND = NBOUND + 1
C           CHECK FOR INTERIOR POINTS ON BOUNDARY OF GTYPE ARRAY
C           IF FOUND IT MEANS THERE IS SOMETHING WRONG
      DO 60 I = 1,NGDIMX
         IF (GTYPE(I,1).EQ.NSIDE .OR. GTYPE(I,NGDIMY).EQ.
     .       NSIDE) GO TO 80
   60 CONTINUE
      DO 70 J = 1,NGDIMY
         IF (GTYPE(1,J).EQ.NSIDE .OR. GTYPE(NGDIMX,J).EQ.
     .       NSIDE) GO TO 80
   70 CONTINUE
C                       NO PROBLEM DETECTED
      RETURN
C                       SERIOUS PROBLEM DETECTED
   80 CONTINUE
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9031)
      FATAL = .TRUE.
      RETURN

 9001 FORMAT (///1H ,19 (1H-)/1H ,18H DOMAIN PROCESSOR /1H ,
     .                      19 (1H-)//5X,
     .                      33H D O M A I N   P R O C E S S O R )
 9011 FORMAT (/5 (3H **),39H FATAL ERROR, MUST HAVE AT LEAST 2 GRID,
     .                 18H LINES IN X AND Y )
 9021 FORMAT (/5 (3H **),32H FATAL ERROR, X AND Y GRID LINES/9X,
     .                 19H MUST BE INCREASING)
 9031 FORMAT (/5 (3H **),33H FATAL ERROR, FOUND INTERIOR PTS /4X,
     .                 43HON THE EDGE OF THE GRID.  PROBABLE CAUSES  /
     .                 4X,40H   BOUNDARY ORIENTATION IS WRONG        /
     .                 4X,
     .                47H   BOUNDARY STARTS WHERE INTERIOR IS TOO THIN
     .                 /4X,
     .                 46H   BOUNDARY OSCILLATES TOO RAPIDLY SOME PLACE
     .                 )
      END
      SUBROUTINE SETUP(COEF,IDCOEF,MXNEQ,MXNCOE,NUMBEQ,NUMCOE,ABD,LDA,
     .                 UNKNWN,NBANDU,NBANDL)
C
C      INTERFACE TO SOLUTION MODULE
C
      REAL COEF(MXNEQ,*),UNKNWN(*),ABD(LDA,*)
      INTEGER IDCOEF(MXNEQ,*)
C
C  FIND BANDWIDTHS AND SET MATRIX SIZE
C
      IF (NBANDU*NBANDL.EQ.0) CALL BNDWTH(IDCOEF,MXNEQ,NUMBEQ,NUMCOE,
     .                             NBANDU,NBANDL)
      NROW = 2*NBANDL + NBANDU + 1
      NCOL = NUMBEQ
      KWORK = NROW*NCOL + NCOL
      WRITE (6,9001) NUMBEQ,NBANDL,NBANDU,KWORK
C
C  ZERO OUT ABD
C
      DO 20 J = 1,NCOL
         DO 10 I = 1,NROW
            ABD(I,J) = 0.0
   10    CONTINUE
   20 CONTINUE
C
C  LOAD ABD AND RIGHT SIDE
C
      M = NBANDL + NBANDU + 1
      DO 40 I = 1,NUMBEQ
         UNKNWN(I) = COEF(I,MXNCOE)
         DO 30 JJ = 1,NUMCOE
            J = IDCOEF(I,JJ)
            IF (J.EQ.0) GO TO 30
            K = I - J + M
            ABD(K,J) = COEF(I,JJ)
   30    CONTINUE
   40 CONTINUE
C
C
      RETURN

 9001 FORMAT (/1H ,38H   INTERFACE TO LINEAR EQUATION SOLVER/1H ,
     .  22H   NUMBER OF EQUATIONS,I9/1H ,18H   LOWER BANDWIDTH,I13/1H ,
     .  18H   UPPER BANDWIDTH,I13/1H ,21H   REQUIRED WORKSPACE,I10)
      END
      SUBROUTINE SOLVE(ABD,LDA,NUMBEQ,UNKNWN,PIVOTS,NBANDU,NBANDL)
C
C  WAYNE R. DYKSEN, JANUARY 1982
C
      REAL ABD(LDA,*),UNKNWN(*),PIVOTS(*)
C
C  FACTOR THE MATRIX
C
      WRITE (6,9001)
      CALL FACTR(ABD,LDA,NUMBEQ,NBANDL,NBANDU,PIVOTS,INFO)
      IF (INFO.NE.0) GO TO 10
C
C  SOLVE THE SYSTEM
C
      CALL BACKSU(ABD,LDA,NUMBEQ,NBANDL,NBANDU,PIVOTS,UNKNWN,0)
C
      WRITE (6,9011)
      RETURN
C
C  ERROR EXIT  --  DIVISION BY ZERO IN BACKSU
C
   10 CONTINUE
      WRITE (6,9021) INFO,INFO
      STOP
C
C
 9001 FORMAT (/1H ,25H   BAND GAUSS ELIMINATION)
 9011 FORMAT (1H ,23H   EXECUTION SUCCESSFUL)
 9021 FORMAT (1H ,36H    THE BAND FACTOR ROUTINE EXECUTED/1H ,
     .  35H    SUCCESSFULLY, BUT THE BAND BACK/1H ,
     .  40H    SOLVE ROUTINE WILL DIVIDE BY ZERO IF/1H ,
     .  36H    CALLED.  THE DIAGONAL ELEMENT IN/1H ,8H    ROW ,I10,
     .  20H OF THE UPPER FACTOR/1H ,12H    IS ZERO.)
      END
      SUBROUTINE OUTPUT(OUTFNC,OUTTYP,TABX,NTABX,TABY,NTABY,WKSPAC,
     .                  GTYPE,NGRDXD,NGRDYD,UNKNWN,RECTAN)
C
C
C  PURPOSE
C
C       CONTROLS THE TYPE OF OUTPUT OF THE COLLOCATION APPROX.
C
C  PARAMETERS
C
C       OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT
C               =1 APPROXIMATE SOLUTION
C               =2 ERROR (TRUE MUST BE SUPPLIED)
C               =3 RESIDUAL
C
C       OUTTYP - INDICATES WHAT INFORMATION ABOUT THE
C                OUTPUT FUNCTION IS TO BE PRINTED
C               =1 PRINT MAX, L1, L2 NORMS OF FUNCTION
C                  BASED ON DISCRETIZATION GRID
C               =2 PRINT MAX, L1, L2 NORMS OF FUNCTION
C                  BASED ON GIVEN GRID (TABX,TABY)
C               =3 PRINT TABLE OF FUNCTION ON DISC. GRID
C               =4 PRINT TABLE OF FUNCTION ON GIVEN GRID
C
C       TABX, TABY - X AND Y COORDINATES FOR OUTTYP 2 AND 4
C
C       NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY
C
C  ALL OTHER PARAMETERS ARE PASSED TO SUBROUTINES
C
C
C
      REAL TABX(NTABX),TABY(NTABY),WKSPAC(NTABX),UNKNWN(*)
      INTEGER OUTFNC,OUTTYP,GTYPE(NGRDXD,NGRDYD)
      LOGICAL RECTAN
C
      COMMON /GRIDXZ/GRIDX(1)
      COMMON /GRIDYZ/GRIDY(1)
      COMMON /PROBI/NGRIDX,NGRIDY
C
      GO TO (10,20,30,40),OUTTYP
C
C PRINT NORMS OF THE FUNCTION ON DISCRETIZATION GRID
C
   10 CALL FNCMAX(OUTFNC,GRIDX,NGRIDX,GRIDY,NGRIDY,GTYPE,NGRDXD,NGRDYD,
     .            UNKNWN,RECTAN)
      RETURN
C
C PRINT NORMS OF FUNCTION ON GIVEN GRID
C
   20 CALL FNCMAX(OUTFNC,TABX,NTABX,TABY,NTABY,GTYPE,NGRDXD,NGRDYD,
     .            UNKNWN,RECTAN)
      RETURN
C
C PRINT TABLE OF FUNCTION ON DISCRETIZATION GRID
C
   30 CALL TABLER(OUTFNC,GRIDX,NGRIDX,GRIDY,NGRIDY,WKSPAC,GTYPE,NGRDXD,
     .            NGRDYD,UNKNWN,RECTAN)
      RETURN
C
C PRINT TABLE OF FUNCTION ON GIVEN GRID
C
   40 CALL TABLER(OUTFNC,TABX,NTABX,TABY,NTABY,WKSPAC,GTYPE,NGRDXD,
     .            NGRDYD,UNKNWN,RECTAN)
      RETURN
      END
C//////////////////////////////////////////////////////////////////////
C///////////////// END OF LOGICAL FILE 1 ////////////////////////////////
C//////////////////////////////////////////////////////////////////////
C//////////////////////////////////////////////////////////////////////
C//////////////////// ALGORITHM  GENCOL //////////////////////////////
C/////////////////////////////////////////////////////////////////////
C>>>>>>>>>>>> LOGICAL FILE 2 : RELATED ROUTINES TO SOLVE S/R<<<<<<<<<<<
C/////////////////////////////////////////////////////////////////////
      SUBROUTINE BNDWTH(IDCOEF,MXNEQ,NUMBEQ,NUMCOE,NBANDU,NBANDL)
C
C  BNDWTH COMPUTES THE BANDWIDTH OF THE LINEAR SYSTEM WHOSE IDS
C  ARE STORED IN THE ARRAY IDCOEF
C
C  NBANDU = NUMBER OF BANDS ABOVE DIAGONAL
C  NBANDL = NUMBER OF BANDS BELOW DIAGONAL
C
      INTEGER IDCOEF(MXNEQ,*)
C
      NBANDU = 0
      NBANDL = 0
C
      DO 20 I = 1,NUMBEQ
         DO 10 JJ = 1,NUMCOE
            J = IDCOEF(I,JJ)
            IF (J.EQ.0) GO TO 10
            JIDIFF = J - I
            NBANDU = MAX0(NBANDU,JIDIFF)
            NBANDL = MIN0(NBANDL,JIDIFF)
   10    CONTINUE
   20 CONTINUE
      NBANDL = IABS(NBANDL)
C
      RETURN
      END
      SUBROUTINE BACKSU(ABD,LDA,N,ML,MU,PIVOTS,B,JOB)
C
C  P U R P O S E
C
C    BACKSU SOLVES THE REAL BAND SYSTEM A * X = B OR TRANS(A) * X = B
C    USING THE FACTORS COMPUTED BY FACTOR.
C
C  D E S C R I P T I O N
C
C    BACKSU IS A MODIFED VERSION OF THE LINPACK GENERAL BAND SOLVE
C    ROUTINE SGBSL WHICH USES AN INTEGER RATHER THAN A REAL PIVOT
C    VECTOR.
C
C  R E F E R E N C E S
C
C    J. J. DONGARRA, C. B. MOLER, J. R. BUNCH, AND G. W. STEWART,
C    LINPACK USERS' GUIDE. SIAM, PHILADELPHIA, 1979.
C
C  A U T H O R
C
C    WAYNE R. DYKSEN
C    DEPARTMENT OF COMPUTER SCIENCE
C    PURDUE UNIVERSITY
C    WEST LAFAYETTE, INDIANA  47097
C    317-494-6001
C
C  V E R S I O N
C
C    JANUARY 1982
C
C---------------------------------------------------------------------
C
      INTEGER LDA,N,ML,MU,JOB
      REAL ABD(LDA,*),B(*),PIVOTS(*)
C
C     BACKSU SOLVES THE REAL BAND SYSTEM
C     A * X = B  OR  TRANS(A) * X = B
C     USING THE FACTORS COMPUTED BY FACTOR.
C
C     ON ENTRY
C
C        ABD     REAL(LDA, N)
C                THE OUTPUT FROM FACTOR
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  ABD .
C
C        N       INTEGER
C                THE ORDER OF THE ORIGINAL MATRIX.
C
C        ML      INTEGER
C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
C
C        MU      INTEGER
C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C
C        PIVOTS  REAL(N)
C                THE PIVOT VECTOR FROM FACTOR.
C
C        B       REAL(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C                = NONZERO   TO SOLVE  TRANS(A)*X = B , WHERE
C                            TRANS(A)  IS THE TRANSPOSE.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
C        CALLED CORRECTLY AND IF FACTOR HAS SET INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL FACTOR(ABD,LDA,N,ML,MU,PIVOTS,INFO)
C           IF (INFO .NE. 0) GO TO ...
C           DO 10 J = 1, P
C              CALL BACKSU(ABD,LDA,N,ML,MU,PIVOTS,C(1,J),0)
C        10 CONTINUE
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS SAXPY,SDOT
C     FORTRAN MIN0
C
C     INTERNAL VARIABLES
C
      REAL SDOT,T
      INTEGER K,KB,L,LA,LB,LM,M,NM1
C
      M = MU + ML + 1
      NM1 = N - 1
      IF (JOB.NE.0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE L*Y = B
C
      IF (ML.EQ.0) GO TO 30
      IF (NM1.LT.1) GO TO 30
      DO 20 K = 1,NM1
         LM = MIN0(ML,N-K)
         L = PIVOTS(K)
         T = B(L)
         IF (L.EQ.K) GO TO 10
         B(L) = B(K)
         B(K) = T
   10    CONTINUE
         CALL SAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
   20 CONTINUE
   30 CONTINUE
C
C        NOW SOLVE  U*X = Y
C
      DO 40 KB = 1,N
         K = N + 1 - KB
         B(K) = B(K)/ABD(M,K)
         LM = MIN0(K,M) - 1
         LA = M - LM
         LB = K - LM
         T = -B(K)
         CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1)
   40 CONTINUE
      GO TO 100

   50 CONTINUE
C
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
C        FIRST SOLVE  TRANS(U)*Y = B
C
      DO 60 K = 1,N
         LM = MIN0(K,M) - 1
         LA = M - LM
         LB = K - LM
         T = SDOT(LM,ABD(LA,K),1,B(LB),1)
         B(K) = (B(K)-T)/ABD(M,K)
   60 CONTINUE
C
C        NOW SOLVE TRANS(L)*X = Y
C
      IF (ML.EQ.0) GO TO 90
      IF (NM1.LT.1) GO TO 90
      DO 80 KB = 1,NM1
         K = N - KB
         LM = MIN0(ML,N-K)
         B(K) = B(K) + SDOT(LM,ABD(M+1,K),1,B(K+1),1)
         L = PIVOTS(K)
         IF (L.EQ.K) GO TO 70
         T = B(L)
         B(L) = B(K)
         B(K) = T
   70    CONTINUE
   80 CONTINUE
   90 CONTINUE
  100 CONTINUE
      RETURN
      END
      SUBROUTINE FACTR(ABD,LDA,N,ML,MU,PIVOTS,INFO)
C
C  P U R P O S E
C
C    FACTOR FACTORS A REAL BAND MATRIX BY GAUSS ELIMINATION.
C
C  D E S C R I P T I O N
C
C    FACTOR FACTORS A REAL BAND MATRIX BY GAUSS ELIMINATION WITH
C    PARTIAL PIVOTING AND ROW EQUILIBRATION. FACTOR IS A MODIFIED
C    VERSION OF THE LINPACK GENERAL BAND FACTOR ROUTINE SGBFA WHICH
C    DOES NOT DO ROW EQUILIBRATION.
C
C  R E F E R E N C E S
C
C    J. J. DONGARRA, C. B. MOLER, J. R. BUNCH, AND G. W. STEWART,
C    LINPACK USERS' GUIDE. SIAM, PHILADELPHIA, 1979.
C
C  A U T H O R
C
C    WAYNE R. DYKSEN
C    DEPARTMENT OF COMPUTER SCIENCE
C    PURDUE UNIVERSITY
C    WEST LAFAYETTE, INDIANA  47097
C    317-494-6001
C
C  V E R S I O N
C
C    JANUARY 1982
C
C-----------------------------------------------------------------------
C
C
      INTEGER LDA,N,ML,MU,INFO
      REAL ABD(LDA,*),PIVOTS(*)
C
C     ON ENTRY
C
C        ABD     REAL(LDA, N)
C                CONTAINS THE MATRIX IN BAND STORAGE.  THE COLUMNS
C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  ABD  AND
C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
C                ML+1 THROUGH 2*ML+MU+1 OF  ABD .
C                SEE THE COMMENTS BELOW FOR DETAILS.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  ABD .
C                LDA MUST BE .GE. 2*ML + MU + 1 .
C
C        N       INTEGER
C                THE ORDER OF THE ORIGINAL MATRIX.
C
C        ML      INTEGER
C                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
C                0 .LE. ML .LT. N .
C
C        MU      INTEGER
C                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C                0 .LE. MU .LT. N .
C                MORE EFFICIENT IF  ML .LE. MU .
C     ON RETURN
C
C        ABD     AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        PIVOTS  REAL(N)
C                A REAL VECTOR OF PIVOT INDICES.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C                     INDICATE THAT EGBSL WILL DIVIDE BY ZERO IF
C                     CALLED.
C
C     BAND STORAGE
C
C           IF  A  IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
C           WILL SET UP THE INPUT.
C
C                   ML = (BAND WIDTH BELOW THE DIAGONAL)
C                   MU = (BAND WIDTH ABOVE THE DIAGONAL)
C                   M = ML + MU + 1
C                   DO 20 J = 1, N
C                      I1 = MAX0(1, J-MU)
C                      I2 = MIN0(N, J+ML)
C                      DO 10 I = I1, I2
C                         K = I - J + M
C                         ABD(K,J) = A(I,J)
C                10    CONTINUE
C                20 CONTINUE
C
C           THIS USES ROWS  ML+1  THROUGH  2*ML+MU+1  OF  ABD .
C           IN ADDITION, THE FIRST  ML  ROWS IN  ABD  ARE USED FOR
C           ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
C           THE TOTAL NUMBER OF ROWS NEEDED IN  ABD  IS  2*ML+MU+1 .
C           THE  ML+MU BY ML+MU  UPPER LEFT TRIANGLE AND THE
C           ML BY ML  LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
C
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS SAXPY,SSCAL
C     MODIFIED BLAS ISWMAX
C     FORTRAN MAX0,MIN0
C
C     INTERNAL VARIABLES
C
      REAL T
      INTEGER I,ISAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
      INTEGER KBEG,KEND,ISWMAX,IPVTK
C
C
      M = ML + MU + 1
      INFO = 0
C
C     FIND RECIPROCAL OF LARGEST ELEMENT IN EACH ROW
C     FOR ROW EQUILIBRATION
C
      DO 10 I = 1,N
         PIVOTS(I) = 0.
   10 CONTINUE
C
      DO 30 J = 1,N
         KBEG = M - MIN0(J-1,MU)
         KEND = M + MIN0(N-J,ML)
         DO 20 K = KBEG,KEND
            I = K + J - M
            PIVOTS(I) = AMAX1(PIVOTS(I),ABS(ABD(K,J)))
   20    CONTINUE
   30 CONTINUE
C
      DO 40 I = 1,N
         IF (PIVOTS(I).NE.0.) PIVOTS(I) = 1./PIVOTS(I)
   40 CONTINUE
C
C     ZERO INITIAL FILL-IN COLUMNS
C
      J0 = MU + 2
      J1 = MIN0(N,M) - 1
      IF (J1.LT.J0) GO TO 70
      DO 60 JZ = J0,J1
         I0 = M + 1 - JZ
         DO 50 I = I0,ML
            ABD(I,JZ) = 0.0E0
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      JZ = J1
      JU = 0
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      NM1 = N - 1
      IF (NM1.LT.1) GO TO 170
      DO 160 K = 1,NM1
         KP1 = K + 1
C
C        ZERO NEXT FILL-IN COLUMN
C
         JZ = JZ + 1
         IF (JZ.GT.N) GO TO 90
         IF (ML.LT.1) GO TO 90
         DO 80 I = 1,ML
            ABD(I,JZ) = 0.0E0
   80    CONTINUE
   90    CONTINUE
C
C        FIND L = PIVOT INDEX
C
         LM = MIN0(ML,N-K)
         ISAMAX = ISWMAX(LM+1,ABD(M,K),PIVOTS(K),1)
         L = ISAMAX + M - 1
         K1 = K + ISAMAX - 1
         PIVOTS(K1) = PIVOTS(K)
         PIVOTS(K) = L + K - M
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (ABD(L,K).EQ.0.0E0) GO TO 140
C
C           INTERCHANGE IF NECESSARY
C
         IF (L.EQ.M) GO TO 100
         T = ABD(L,K)
         ABD(L,K) = ABD(M,K)
         ABD(M,K) = T
  100    CONTINUE
C
C           COMPUTE MULTIPLIERS
C
         T = -1.0E0/ABD(M,K)
         CALL SSCAL(LM,T,ABD(M+1,K),1)
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C
         IPVTK = PIVOTS(K)
         JU = MIN0(MAX0(JU,MU+IPVTK),N)
         MM = M
         IF (JU.LT.KP1) GO TO 130
         DO 120 J = KP1,JU
            L = L - 1
            MM = MM - 1
            T = ABD(L,J)
            IF (L.EQ.MM) GO TO 110
            ABD(L,J) = ABD(MM,J)
            ABD(MM,J) = T
  110       CONTINUE
            CALL SAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
  120    CONTINUE
  130    CONTINUE
         GO TO 150

  140    CONTINUE
         INFO = K
  150    CONTINUE
  160 CONTINUE
  170 CONTINUE
      PIVOTS(N) = N
      IF (ABD(M,N).EQ.0.0E0) INFO = N
      RETURN
      END
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
C
C     OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY.
C
      REAL SX(*),SY(*),SA
      IF (N.LE.0 .OR. SA.EQ.0.E0) RETURN
      IF (INCX.EQ.INCY) IF (INCX-1) 10,30,70
   10 CONTINUE
C
C        CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 20 I = 1,N
         SY(IY) = SY(IY) + SA*SX(IX)
         IX = IX + INCX
         IY = IY + INCY
   20 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4.
C
   30 M = N - (N/4)*4
      IF (M.EQ.0) GO TO 50
      DO 40 I = 1,M
         SY(I) = SY(I) + SA*SX(I)
   40 CONTINUE
      IF (N.LT.4) RETURN
   50 MP1 = M + 1
      DO 60 I = MP1,N,4
         SY(I) = SY(I) + SA*SX(I)
         SY(I+1) = SY(I+1) + SA*SX(I+1)
         SY(I+2) = SY(I+2) + SA*SX(I+2)
         SY(I+3) = SY(I+3) + SA*SX(I+3)
   60 CONTINUE
      RETURN
C
C        CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
C
   70 CONTINUE
      NS = N*INCX
      DO 80 I = 1,NS,INCX
         SY(I) = SA*SX(I) + SY(I)
   80 CONTINUE
      RETURN
      END
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
C     RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY.
C
      REAL SX(1),SY(1)
      SDOT = 0.0E0
      IF (N.LE.0) RETURN
      IF (INCX.EQ.INCY) IF (INCX-1) 10,30,70
   10 CONTINUE
C
C        CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS.
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 20 I = 1,N
         SDOT = SDOT + SX(IX)*SY(IY)
         IX = IX + INCX
         IY = IY + INCY
   20 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5.
C
   30 M = N - (N/5)*5
      IF (M.EQ.0) GO TO 50
      DO 40 I = 1,M
         SDOT = SDOT + SX(I)*SY(I)
   40 CONTINUE
      IF (N.LT.5) RETURN
   50 MP1 = M + 1
      DO 60 I = MP1,N,5
         SDOT = SDOT + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) +
     .          SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4)
   60 CONTINUE
      RETURN
C
C        CODE FOR POSITIVE EQUAL INCREMENTS .NE.1.
C
   70 CONTINUE
      NS = N*INCX
      DO 80 I = 1,NS,INCX
         SDOT = SDOT + SX(I)*SY(I)
   80 CONTINUE
      RETURN
      END
      SUBROUTINE SSCAL(N,SA,SX,INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SA,SX(*)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1) GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I = 1,NINCX,INCX
         SX(I) = SA*SX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1,M
         SX(I) = SA*SX(I)
   30 CONTINUE
      IF (N.LT.5) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
         SX(I) = SA*SX(I)
         SX(I+1) = SA*SX(I+1)
         SX(I+2) = SA*SX(I+2)
         SX(I+3) = SA*SX(I+3)
         SX(I+4) = SA*SX(I+4)
   50 CONTINUE
      RETURN
      END
      INTEGER FUNCTION ISWMAX(N,SX,WEIGHT,INCX)
C
C     FINDS THE INDEX OF ELEMENT HAVING MAX. WEIGHTED ABSOLUTE VALUE.
C     WAYNE DYKSEN, JANUARY 1982
C
      REAL SX(1),WEIGHT(1),SMAX
      INTEGER I,INCX,IX,N
C
      ISWMAX = 0
      IF (N.LT.1) RETURN
      ISWMAX = 1
      IF (N.EQ.1) RETURN
      IF (INCX.EQ.1) GO TO 30
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      IX = 1
      SMAX = ABS(SX(1)*WEIGHT(1))
      IX = IX + INCX
      DO 20 I = 2,N
         IF (ABS(SX(IX)*WEIGHT(I)).LE.SMAX) GO TO 10
         ISWMAX = I
         SMAX = ABS(SX(IX)*WEIGHT(I))
   10    IX = IX + INCX
   20 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
   30 SMAX = ABS(SX(1)*WEIGHT(1))
      DO 40 I = 2,N
         IF (ABS(SX(I)*WEIGHT(I)).LE.SMAX) GO TO 40
         ISWMAX = I
         SMAX = ABS(SX(I)*WEIGHT(I))
   40 CONTINUE
      RETURN
      END
C///////////////////////////////////////////////////////////////////////
C//////////////////// END OF LOGICAL FILE 2 ////////////////////////////
C///////////////////////////////////////////////////////////////////////
C////////////////////////////////////////////////////////////////////////
C//////////////////// ALGORITHM  GENCOL ////////////////////////////////
C///////////////////////////////////////////////////////////////////////
C>>>>>>>>>> LOGICAL FILE 3 : RELATED SUBPROGRAMS TO S/R OUTPUT<<<<<<<<<<
C///////////////////////////////////////////////////////////////////////
      SUBROUTINE SETGRD(GRIDX,GRIDY,NGRIDX,NGRIDY,AX,BX,AY,BY,IGRID)
C
      REAL GRIDX(*),GRIDY(*)
C
      IF (IGRID.LT.0) GO TO 20
C
      DX = BX - AX
      DY = BY - AY
      NGRIDX = IGRID
      NGRIDY = IGRID
C
      AIGRM1 = IGRID - 1
      DO 10 I = 1,IGRID
         AIM1 = I - 1
         GRIDX(I) = AX + (AIM1)*DX/ (AIGRM1)
         GRIDY(I) = AY + (AIM1)*DY/ (AIGRM1)
   10 CONTINUE
C
      RETURN
C
   20 CONTINUE
      RETURN
      END
      SUBROUTINE FNDINT(X1A,Y1A,X1B,Y1B,X2A,Y2A,X2B,Y2B,X0,Y0,NOINT)
C
C ...............................................................
C
C     E L L P A C K    7 8    O U T P U T    M O D U L E
C
C     SUBROUTINE  FNDINT
C
C     PURPOSE
C         TO FIND THE INTERSECTION POINT (X0,Y0) OF LINE L1
C         DEFINED BY (X1A,Y1A) AND (X1B,Y1B) AND LINE L2
C         DEFINED BY (X2A,Y2A) AND (X2B,Y2B).
C
C     METHOD
C         THE SLOPE INTERCEPT FORM ( Y = M*X + B ) IS USED.
C
C ...............................................................
C
      REAL X1A,Y1A,X1B,Y1B,X2A,Y2A,X2B,Y2B,X0,Y0,M1,B1,M2,B2,DX,DM,
     .     EPSGRD
      COMMON /NUMCON/EPSGRD
      NOINT = 0
      NCASE = 1
C
C         IS L1 A VERTICAL LINE
C
      DX = X1B - X1A
      IF (ABS(DX).GT.EPSGRD) GO TO 10
      NCASE = 2
      GO TO 20

   10 CONTINUE
      M1 = (Y1B-Y1A)/DX
      B1 = Y1A - M1*X1A
   20 CONTINUE
C
C         IS L2 A VERTICAL LINE
C
      DX = X2B - X2A
      IF (ABS(DX).GT.EPSGRD) GO TO 30
      NCASE = NCASE + 2
      GO TO 40

   30 CONTINUE
      M2 = (Y2B-Y2A)/DX
      B2 = Y2A - M2*X2A
C
C         BRANCH ON NCASE
C
C         NCASE = 1 - BOTH LINES ARE NOT VERTICAL
C         NCASE = 2 - L1 IS VERTICAL, L2 IS NOT
C         NCASE = 3 - L2 IS VERTICAL, L1 IS NOT
C         NCASE = 4 - BOTH LINES ARE VERTICAL
C
   40 CONTINUE
      GO TO (50,60,70,80),NCASE
C
C         BOTH LINES ARE NOT VERTICAL
C
   50 CONTINUE
      DM = M2 - M1
      IF (ABS(DM).LE.EPSGRD) GO TO 80
      X0 = (B1-B2)/DM
      Y0 = M1*X0 + B1
      GO TO 90
C
C         L1 IS VERTICAL, L2 IS NOT
C
   60 CONTINUE
      X0 = X1A
      Y0 = M2*X0 + B2
      GO TO 90
C
C         L2 IS VERTICAL, L1 IS NOT
C
   70 CONTINUE
      X0 = X2A
      Y0 = M1*X0 + B1
      GO TO 90
C
C        L1 AND L2 ARE PARALLEL
C
   80 CONTINUE
      NOINT = 1
C
   90 CONTINUE
      RETURN
      END
      SUBROUTINE INSQAR(X,Y,GT,GX,GY,XBOUND,YBOUND,NBNDPT,IPNTYP,IERR)
C
C ...............................................................
C
C     E L L P A C K    7 8    O U T P U T    M O D U L E
C
C     SUBROUTINE  INSQAR
C
C     PURPOSE
C         TO DETERMINE IF THE POINT (X,Y) IS INSIDE THE REGION,
C         ON THE BOUNDARY, OR OUTSIDE THE REGION WHEN (X,Y) IS
C         INSIDE A GRID SQUARE.
C
C     PARAMETERS
C         X      - THE X COORDINATE OF THE POINT
C         Y      - THE Y COORDINATE OF THE POINT
C         GT     - THE GTYPES OF THE FOUR CORNER POINTS OF THE
C                  SQUARE BEGINNING WITH THE LOWER LEFT HAND
C                  POINT AND PROCEEDING COUNTERCLOCKWISE.
C         GX     - GX(1) IS THE X COORDINATE OF THE LEFT GRID
C                  LINE BOUNDING THE GRID SQUARE WHILE GX(2)
C                  IS THE X COORDINATE OF THE RIGHT BOUNDARY.
C         GY     - GY(1) IS THE Y COORDINATE OF THE BOTTOM OF
C                  THE GRID SQUARE WHILE GY(2) IS THE Y COORDI-
C                  NATE OF THE TOP.
C         IPNTYP - 1, IF (X,Y) IS INSIDE THE REGION
C                  2, IF (X,Y) IS ON THE BOUNDARY
C                  3, IF (X,Y) IS OUTSIDE THE REGION
C         IERR   - ERROR RETURN CODE.
C
C ...............................................................
C
      REAL X,Y,GX(4),GY(4),XBOUND(NBNDPT),YBOUND(NBNDPT),XBK,YBK,XBKP1,
     .     YBKP1,PLEN,BLEN,GXL,GYL,XINT,YINT
      INTEGER GT(4),IBOUND(8)
      COMMON /NUMCON/EPSGRD
C
C         NINSID - THE NUMBER OF NEIGHBORING GRID POINTS
C                  INSIDE THE DOMAIN
C         NOUTSD - THE NUMBER OF NEIGHBORING GRID POINTS
C                  OUTSIDE THE DOMAIN
C
      NINSID = 0
      NOUTSD = 0
C
      DO 20 K = 1,4
         IF (GT(K).LT.999) GO TO 10
         NINSID = NINSID + 1
         GO TO 20

   10    CONTINUE
         IF (GT(K).LE.0) NOUTSD = NOUTSD + 1
   20 CONTINUE
C
      NCASE = 5*NOUTSD + NINSID + 1
C
C         BRANCH ON NCASE
C         -----------------------------------
C          NOUTSD   NINSID   NCASE    (X,Y)
C         -----------------------------------
C
C             0        0       1     OUTSIDE
C             0        1       2     INSIDE
C             0        2       3     INSIDE
C             0        3       4     INSIDE
C             0        4       5     INSIDE
C             1        0       6     OUTSIDE
C             1        1       7     UNKNOWN
C             1        2       8     UNKNOWN
C             1        3       9     UNKNOWN
C             1        4      10     IMPOSS.
C             2        0      11     OUTSIDE
C             2        1      12     UNKNOWN
C             2        2      13     UNKNOWN
C             2        3      14     IMPOSS.
C             2        4      15     IMPOSS.
C             3        0      16     OUTSIDE
C             3        1      17     UNKNOWN
C             3        2      18     IMPOSS.
C             3        3      19     IMPOSS.
C             3        4      20     IMPOSS.
C             4        0      21     OUTSIDE
C
      GO TO (120,100,100,100,100,120,30,30,30,
     .       130,120,30,30,130,130,120,30,130,
     .       130,130,120),NCASE
C
C         THIS IS THE CASE WHERE AT LEAST ONE OF THE
C         NEIGHBORING GRID POINTS IS INSIDE AND AT
C         LEAST ONE IS OUTSIDE.
C
   30 CONTINUE
      KB = 0
C
      DO 40 K = 1,NBNDPT
         IF (XBOUND(K).LT.GX(1)-EPSGRD) GO TO 40
         IF (YBOUND(K).LT.GY(1)-EPSGRD) GO TO 40
         IF (XBOUND(K).GT.GX(3)+EPSGRD) GO TO 40
         IF (YBOUND(K).GT.GY(3)+EPSGRD) GO TO 40
         KB = KB + 1
         IBOUND(KB) = K
   40 CONTINUE
C
      IF (KB.GE.2) GO TO 50
      IERR = 3
      GO TO 120
C
   50 CONTINUE
      KSTOP = KB - 1
C
      DO 90 K = 1,KSTOP
         IBK = IBOUND(K)
         IBKP1 = IBOUND(K+1)
         IF (IBKP1-IBK.EQ.1) GO TO 60
         IF (IBK.NE.1) GO TO 90
         IF (IBKP1.NE.NBNDPT) GO TO 90
   60    CONTINUE
         XBK = XBOUND(IBK)
         YBK = YBOUND(IBK)
         XBKP1 = XBOUND(IBKP1)
         YBKP1 = YBOUND(IBKP1)
         IF (ABS(XBK-XBKP1).LT.EPSGRD) GO TO 90
         IF (ABS(YBK-YBKP1).LT.EPSGRD) GO TO 90
C
         DO 80 L = 1,4
            IF ((GT(L).GT.0) .AND. (GT(L).LT.999)) GO TO 80
            GXL = GX(L)
            GYL = GY(L)
            CALL FNDINT(X,Y,GXL,GYL,XBK,YBK,XBKP1,YBKP1,XINT,YINT,NOINT)
            IF (NOINT.NE.0) GO TO 80
            IF (XINT.GT.AMAX1(XBK,XBKP1)+EPSGRD) GO TO 80
            IF (XINT.LT.AMIN1(XBK,XBKP1)-EPSGRD) GO TO 80
            IF (YINT.GT.AMAX1(YBK,YBKP1)+EPSGRD) GO TO 80
            IF (YINT.LT.AMIN1(YBK,YBKP1)-EPSGRD) GO TO 80
C
C                 PLEN - THE DISTANCE FROM THE GRID POINT TO (X,Y)
C
C                 BLEN - THE DISTANCE FROM THE GRID POINT TO THE
C                        BOUNDARY INTERSECTION POINT
C
            PLEN = (GXL-X)**2 + (GYL-Y)**2
            BLEN = (GXL-XINT)**2 + (GYL-YINT)**2
C
C                 NOW DETERMINE THE RELATIONSHIP
C                 OF (X,Y) TO THE BOUNDARY.
C
            IF (ABS(PLEN-BLEN).LE.EPSGRD) GO TO 110
            IF (GT(L).LE.0) GO TO 70
            IF (PLEN.LT.BLEN) GO TO 100
            GO TO 120

   70       CONTINUE
            IF (PLEN.GT.BLEN) GO TO 100
            GO TO 120

   80    CONTINUE
   90 CONTINUE
      GO TO 120
C
C         (X,Y) IS INSIDE
C
  100 CONTINUE
      IPNTYP = 1
      GO TO 130
C
C         (X,Y) IS ON THE BOUNDARY
C
  110 CONTINUE
      IPNTYP = 2
      GO TO 130
C
C         (X,Y) IS OUTSIDE
C
  120 CONTINUE
      IPNTYP = 3
C
  130 CONTINUE
      RETURN
      END
      SUBROUTINE OUTSID(X,Y,GTYPE,GRIDX,GRIDY,XBOUND,YBOUND,NGRDXD,
     .                  NGRDYD,NBNDPT,IPNTYP,IERR)
C
C ...............................................................
C
C     E L L P A C K    7 8    O U T P U T    M O D U L E
C
C     SUBROUTINE  OUTSID
C
C     PURPOSE
C         TO DETERMINE IF THE POINT (X,Y) IS INSIDE THE REGION,
C         ON THE BOUNDARY, OR OUTSIDE THE REGION.
C
C     PARAMETERS
C         X      - THE X COORDINATE OF THE POINT
C         Y      - THE Y COORDINATE OF THE POINT
C         IPNTYP - 1, IF (X,Y) IS INSIDE THE REGION
C                  2, IF (X,Y) IS ON THE BOUNDARY
C                  3, IF (X,Y) IS OUTSIDE THE REGION
C         IERR   - ERROR RETURN CODE.
C                  1, IF NO NEIGHBORING VERTICAL GRID
C                     LINES WERE FOUND.
C                  2, IF NO NEIGHBORING HORIZONTAL GRID
C                     LINES WERE FOUND.
C                  3, IF AN ERROR OCCURED IN SUBROUTINE INSQAR.
C                  4, IF AN ERROR OCCURED IN SUBROUTINE ONLINE.
C
C ...............................................................
C
      REAL X,Y,GRIDX(NGRDXD),GRIDY(NGRDYD),XBOUND(NBNDPT),
     .     YBOUND(NBNDPT),GX(4),GY(4),EPSGRD,AX,BX,AY,BY
      INTEGER GTYPE(NGRDXD,NGRDYD),GT(4),IPNTYP,IPTCOM
      COMMON /PROBR/AX,BX,AY,BY
      COMMON /NUMCON/EPSGRD
      COMMON /PNTYPE/IPTCOM
C
C         CHECK FOR TRIVIAL CASES
C
      IERR = 0
      IF (X.LT.AX-EPSGRD) GO TO 130
      IF (Y.LT.AY-EPSGRD) GO TO 130
      IF (X.GT.BX+EPSGRD) GO TO 130
      IF (Y.GT.BY+EPSGRD) GO TO 130
      NCASE = 1
C
C         FIND VERTICAL GRID LINES SUCH THAT
C
C         GRIDX(I-1) .LT. X .LT. GRIDX(I)
C
C
      DO 20 I = 1,NGRDXD
         IF (ABS(GRIDX(I)-X).GT.EPSGRD) GO TO 10
         NCASE = 2
         GO TO 30

   10    CONTINUE
         IF (X.LE.GRIDX(I)) GO TO 30
   20 CONTINUE
C
      IERR = 1
      GO TO 130
C
C         FIND HORIZONTAL GRID LINES SUCH THAT
C
C         GRIDY(J-1) .LT. Y .LT. GRIDY(J)
C
   30 CONTINUE
C
      DO 50 J = 1,NGRDYD
         IF (ABS(GRIDY(J)-Y).GT.EPSGRD) GO TO 40
         NCASE = NCASE + 2
         GO TO 60

   40    CONTINUE
         IF (Y.LE.GRIDY(J)) GO TO 60
   50 CONTINUE
C
      IERR = 2
      GO TO 130
C
C         BRANCH ON NCASE
C
C         NCASE = 1 - (X,Y) IS STRICTLY ON THE INTERIOR
C                           OF A GRID SQUARE
C         NCASE = 2 - (X,Y) IS ON A VERTICAL GRID LINE
C         NCASE = 3 - (X,Y) IS ON A HORIZONTAL GRID LINE
C         NCASE = 4 - (X,Y) IS A GRID POINT
C
   60 CONTINUE
      GO TO (70,80,90,100),NCASE
C
C         (X,Y) IS AN INTERIOR POINT
C
   70 CONTINUE
      GX(1) = GRIDX(I-1)
      GX(2) = GRIDX(I)
      GX(3) = GX(2)
      GX(4) = GX(1)
      GY(1) = GRIDY(J-1)
      GY(2) = GY(1)
      GY(3) = GRIDY(J)
      GY(4) = GY(3)
      GT(1) = GTYPE(I-1,J-1)
      GT(2) = GTYPE(I,J-1)
      GT(3) = GTYPE(I,J)
      GT(4) = GTYPE(I-1,J)
      CALL INSQAR(X,Y,GT,GX,GY,XBOUND,YBOUND,NBNDPT,IPNTYP,IERR)
      GO TO 140
C
C         (X,Y) IS ON A VERTICAL GRID LINE
C
   80 CONTINUE
      GY(1) = GRIDY(J-1)
      GY(2) = GRIDY(J)
      GT(1) = GTYPE(I,J-1)
      GT(2) = GTYPE(I,J)
      CALL ONLINE(X,Y,GT,GY,XBOUND,YBOUND,NBNDPT,IPNTYP,IERR)
      GO TO 140
C
C         (X,Y) IS ON A HORIZONTAL GRID LINE
C
   90 CONTINUE
      GX(1) = GRIDX(I-1)
      GX(2) = GRIDX(I)
      GT(1) = GTYPE(I-1,J)
      GT(2) = GTYPE(I,J)
      CALL ONLINE(Y,X,GT,GX,YBOUND,XBOUND,NBNDPT,IPNTYP,IERR)
      GO TO 140
C
C         (X,Y) IS A GRID POINT
C
  100 CONTINUE
      IF (GTYPE(I,J).GT.999) GO TO 110
      IF (GTYPE(I,J).GT.0) GO TO 120
      GO TO 130
C
C         (X,Y) IS INSIDE
C
  110 CONTINUE
      IPNTYP = 1
      GO TO 140
C
C         (X,Y) IS ON THE BOUNDARY
C
  120 CONTINUE
      IPNTYP = 2
      GO TO 140
C
C         (X,Y) IS OUTSIDE
C
  130 CONTINUE
      IPNTYP = 3
C
  140 CONTINUE
      IPTCOM = IPNTYP
      RETURN
      END
      SUBROUTINE ONLINE(X,Y,GT,GY,XBOUND,YBOUND,NBNDPT,IPNTYP,IERR)
C
C ...............................................................
C
C     E L L P A C K    7 8    O U T P U T    M O D U L E
C
C     SUBROUTINE  ONLINE
C
C     PURPOSE
C         TO DETERMINE IF THE POINT (X,Y) IS INSIDE THE REGION,
C         ON THE BOUNDARY, OR OUTSIDE THE REGION WHEN (X,Y) IS
C         ON A VERTICAL GRID LINE.
C
C     PARAMETERS
C         X      - THE X COORDINATE OF THE POINT
C         Y      - THE Y COORDINATE OF THE POINT
C         GT     - GT(1) IS THE GTYPE OF THE LOWER GRID POINT
C                  (X,GY(1)) WHILE GT(2) IS THE GTYPE OF THE
C                  THE UPPER GRID POINT (X,GY(2))
C         GY     - GY(1) IS THE Y COORDINATE OF THE LOWER GRID
C                  POINT WHILE GY(2) IS THE Y COORDINATE OF THE
C                  UPPER GRID POINT.
C         IPNTYP - 1, IF (X,Y) IS INSIDE THE REGION
C                  2, IF (X,Y) IS ON THE BOUNDARY
C                  3, IF (X,Y) IS OUTSIDE THE REGION
C         IERR   - ERROR RETURN CODE.
C
C ...............................................................
C
      REAL X,Y,GY(2),XBOUND(NBNDPT),YBOUND(NBNDPT),PLEN,BLEN
      INTEGER GT(2),IGT(2)
      COMMON /NUMCON/EPSGRD
C
C         NINSID - THE NUMBER OF NEIGHBORING GRID POINTS
C                  INSIDE THE DOMAIN
C         NOUTSD - THE NUMBER OF NEIGHBORING GRID POINTS
C                  OUTSIDE THE DOMAIN
C
      NINSID = 0
      NOUTSD = 0
C
      DO 20 K = 1,2
         IF (GT(K).LT.999) GO TO 10
         NINSID = NINSID + 1
         GO TO 20

   10    CONTINUE
         IF (GT(K).LE.0) NOUTSD = NOUTSD + 1
   20 CONTINUE
C
      NCASE = 3*NOUTSD + NINSID + 1
C
C         BRANCH ON NCASE
C         -----------------------------------
C          NOUTSD   NINSID   NCASE    (X,Y)
C         -----------------------------------
C
C             0        0       1     BOUNDARY
C             0        1       2     INSIDE
C             0        2       3     INSIDE
C             1        0       4     OUTSIDE
C             1        1       5     UNKNOWN
C             1        2       6     IMPOSS.
C             2        0       7     OUTSIDE
C
      GO TO (100,90,90,110,30,120,110),NCASE
C
C         THIS IS THE CASE WHERE ONE OF THE NEIGHBORING GRID
C         POINTS IS INSIDE WHILE THE OTHER IS OUTSIDE.
C
C         FIRST, TRY TO USE GTYPE INFO SAVED IN GT TO FIND THE
C         BOUNDARY - GRID LINE INTERSECTION POINT.
C
   30 CONTINUE
C
      DO 40 K = 1,2
         KB = MOD(IABS(GT(K)),1000)
         IF (ABS(XBOUND(KB)-X).GT.EPSGRD) GO TO 40
         IF (YBOUND(KB).LT.GY(1)) GO TO 40
         IF (YBOUND(KB).GT.GY(2)) GO TO 40
         GO TO 60

   40 CONTINUE
C
C         USING GTYPE HAS FAILED.
C         NOW TRY A BRUTE FORCE SEARCH.
C
      DO 50 KB = 1,NBNDPT
         IF (ABS(XBOUND(KB)-X).GT.EPSGRD) GO TO 50
         IF (YBOUND(KB).LT.GY(1)) GO TO 50
         IF (YBOUND(KB).GT.GY(2)) GO TO 50
         GO TO 60

   50 CONTINUE
C
      IERR = 4
      GO TO 120
C
C         PLEN - THE DISTANCE FROM THE INTERIOR POINT TO (X,Y)
C
C         BLEN - THE DISTANCE FROM THE INTERIOR POINT TO THE
C                BOUNDARY - GRID LINE INTERSECTION POINT
C
   60 CONTINUE
      IF (GT(2).GE.999) GO TO 70
      PLEN = ABS(GY(1)-Y)
      BLEN = ABS(GY(1)-YBOUND(KB))
      GO TO 80

   70 CONTINUE
      PLEN = ABS(GY(2)-Y)
      BLEN = ABS(GY(2)-YBOUND(KB))
C
C         NOW DETERMINE THE RELATIONSHIP
C         OF (X,Y) TO THE BOUNDARY.
C
   80 CONTINUE
      IF (ABS(PLEN-BLEN).LE.EPSGRD) GO TO 100
      IF (PLEN.LT.BLEN) GO TO 90
      GO TO 110
C
C         (X,Y) IS INSIDE
C
   90 CONTINUE
      IPNTYP = 1
      GO TO 120
C
C         (X,Y) IS ON THE BOUNDARY
C
  100 CONTINUE
      IPNTYP = 2
      GO TO 120
C
C         (X,Y) IS OUTSIDE
C
  110 CONTINUE
      IPNTYP = 3
C
  120 CONTINUE
      RETURN
      END
      SUBROUTINE TABLER(OUTFNC,TABX,NTABX,TABY,NTABY,TPRINT,GTYPE,
     .                  NGRDXD,NGRDYD,UNKNWN,RECTAN)
C
C
C     PURPOSE
C
C         PRINTS THE VALUES OF THE OUTPUT FUNCTION ON THE GRID
C         DEFINED BY THE ARRAYS TABX AND TABY IF A PAR-
C         TICULAR GRID POINT IS INSIDE THE REGION.
C
C     PARAMETERS
C
C       OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT
C               =1 APPROXIMATE SOLUTION
C               =2 ERROR (TRUE MUST BE SUPPLIED)
C               =3 RESIDUAL
C
C       TABX, TABY - X AND Y COORDINATES FOR OUTTYP 2 AND 4
C
C       NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY
C
C       RECTAN - =TRUE - RECTANGULAR DOMAIN
C                =FALSE- NONRECTANGULAR DOMAIN
C
C       TPRINT - WORKSPACE OF DIMENSION >= NTABX
C
C  ALL OTHER PARAMETERS ARE PASSED TO SUBROUTINES
C
C
      REAL TABX(NTABX),TABY(NTABY),DERVSL(6),UNKNWN(1),TPRINT(NTABX)
C
      INTEGER GTYPE(NGRDXD,NGRDYD),OUTFNC
C
      LOGICAL RECTAN
C
      COMMON /GRIDXZ/GRIDX(1)
      COMMON /GRIDYZ/GRIDY(1)
      COMMON /XBOUZZ/XBOUND(1)
      COMMON /YBOUZZ/YBOUND(1)
      COMMON /BNDRY/IPIECE,NBOUND,NBNDPT
      COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW
C
C
C       PRINT HEADING
C
      IF (OUTFNC.EQ.1) WRITE (MOUTPT,9001) NTABX,NTABY
      IF (OUTFNC.EQ.2) WRITE (MOUTPT,9011) NTABX,NTABY
      IF (OUTFNC.EQ.3) WRITE (MOUTPT,9021) NTABX,NTABY
      WRITE (MOUTPT,9031)
      WRITE (MOUTPT,9051) (TABX(I),I=1,NTABX)
C
C
      DO 20 JUP = 1,NTABY
         J = NTABY - JUP + 1
         Y = TABY(J)
         WRITE (MOUTPT,9041) Y
C
         DO 10 I = 1,NTABX
            TPRINT(I) = 0.0
            X = TABX(I)
            IPNTYP = 0
            IF ( .NOT. RECTAN) CALL OUTSID(X,Y,GTYPE,GRIDX,GRIDY,XBOUND,
     .                              YBOUND,NGRDXD,NGRDYD,NBNDPT,IPNTYP,
     .                              IERR)
            IF (IPNTYP.EQ.3) GO TO 10
            IF (OUTFNC.NE.3) TPRINT(I) = COLAPR(X,Y,UNKNWN,DERVSL,
     .                                   RECTAN)
            IF (OUTFNC.EQ.2) TPRINT(I) = TPRINT(I) - TRUE(X,Y)
            IF (OUTFNC.EQ.3) TPRINT(I) = RESID(X,Y,IPNTYP,UNKNWN,RECTAN)
   10    CONTINUE
C
         WRITE (MOUTPT,9051) (TPRINT(I),I=1,NTABX)
   20 CONTINUE
C
      RETURN

 9001 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X,
     .                          9HTABLE OF ,8HSOLUTION,3H ON,I4,2H X,I4,
     .                          6H  GRID,4X,1H+)
 9011 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X,
     .                          9HTABLE OF ,8HERROR   ,3H ON,I4,2H X,I4,
     .                          6H  GRID,4X,1H+)
 9021 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X,
     .                         9HTABLE OF ,9HRESIDUALS,3H ON,I4,2H X,I4,
     .                          6H  GRID,4X,1H+)
 9031 FORMAT (1H ,10X,1H+,44X,1H+/1H ,10X,46 (1H+)///1H ,4X,
     .  15HX-ABSCISSAE ARE/1H ,4X,15 (1H-))
 9041 FORMAT (/1H ,7X,3HY =,E13.6/1H ,7X,16 (1H-))
 9051 FORMAT ((1H ,3X,E13.6,3X,E13.6,3X,E13.6,3X,E13.6))
      END
      SUBROUTINE FNCMAX(OUTFNC,TABX,NTABX,TABY,NTABY,GTYPE,NGRDXD,
     .                  NGRDYD,UNKNWN,RECTAN)
C
C
C  PURPOSE
C
C    COMPUTE AND PRINT MAX, L1 AND L2 NORMS OF THE OUTPUT
C    FUNCTION BASED ON THE GRID DEFINED BY TABX AND TABY
C
C  PARAMETERS
C
C       OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT
C               =1 APPROXIMATE SOLUTION
C               =2 ERROR (TRUE MUST BE SUPPLIED)
C               =3 RESIDUAL
C
C       TABX, TABY - GRID ON WHICH TO MEASURE NORMS
C
C       NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY
C
C       RECTAN - =TRUE - RECTANGULAR DOMAIN
C                =FALSE- NONRECTANGULAR DOMAIN
C
C
      REAL TABX(NTABX),TABY(NTABY),DERVSL(6),UNKNWN(*)
C
      INTEGER GTYPE(NGRDXD,NGRDYD),OUTFNC
C
      LOGICAL RECTAN
C
      COMMON /GRIDXZ/GRIDX(1)
      COMMON /GRIDYZ/GRIDY(1)
      COMMON /XBOUZZ/XBOUND(1)
      COMMON /YBOUZZ/YBOUND(1)
      COMMON /BNDRY/IPIECE,NBOUND,NBNDPT
      COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW
C
C
      R1NRMI = 0.0
      R1NRM1 = 0.0
      R1NRM2 = 0.0
      GPTS = 0
C
C
      DO 20 I = 1,NTABX
         X = TABX(I)
         DO 10 J = 1,NTABY
            Y = TABY(J)
            IPTYPE = 0
            IF ( .NOT. RECTAN) CALL OUTSID(X,Y,GTYPE,GRIDX,GRIDY,XBOUND,
     .                              YBOUND,NGRDXD,NGRDYD,NBNDPT,IPTYPE,
     .                              IERR)
            IF (IPTYPE.EQ.3) GO TO 10
            IF (OUTFNC.NE.3) FVALUE = COLAPR(X,Y,UNKNWN,DERVSL,RECTAN)
            IF (OUTFNC.EQ.2) FVALUE = ABS(FVALUE-TRUE(X,Y))
            IF (OUTFNC.EQ.3) FVALUE = RESID(X,Y,IPTYPE,UNKNWN,RECTAN)
            GPTS = GPTS + 1.
            R1NRMI = AMAX1(R1NRMI,FVALUE)
            R1NRM1 = R1NRM1 + FVALUE
            R1NRM2 = R1NRM2 + FVALUE**2
   10    CONTINUE
   20 CONTINUE
      R1NRM1 = R1NRM1/GPTS
      R1NRM2 = SQRT(R1NRM2/GPTS)
      IF (OUTFNC.EQ.1) WRITE (MOUTPT,9001) NTABX,NTABY,R1NRMI,NTABX,
     .    NTABY,R1NRM1,NTABX,NTABY,R1NRM2
      IF (OUTFNC.EQ.2) WRITE (MOUTPT,9011) NTABX,NTABY,R1NRMI,NTABX,
     .    NTABY,R1NRM1,NTABX,NTABY,R1NRM2
      IF (OUTFNC.EQ.3) WRITE (MOUTPT,9021) NTABX,NTABY,R1NRMI,NTABX,
     .    NTABY,R1NRM1,NTABX,NTABY,R1NRM2
      RETURN
C
C
 9001 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+  MAX( ABS(,
     .                    8HSOLUTION,7H) ) ON ,I3,3H X ,I3,7H GRID =,
     .                    1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L1 NORM( ,8HSOLUTION,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L2 NORM( ,8HSOLUTION,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    60 (1H+))
 9011 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+  MAX( ABS(,
     .                    8HERROR   ,7H) ) ON ,I3,3H X ,I3,7H GRID =,
     .                    1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L1 NORM( ,8HERROR   ,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L2 NORM( ,8HERROR   ,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    60 (1H+))
 9021 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+  MAX( ABS(,
     .                    8HRESIDUAL,7H) ) ON ,I3,3H X ,I3,7H GRID =,
     .                    1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L1 NORM( ,8HRESIDUAL,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    12H+  L2 NORM( ,8HRESIDUAL,7H)   ON ,I3,3H X ,
     .                    I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X,
     .                    60 (1H+))
      END
      REAL FUNCTION RESID(X,Y,IPTCOM,UNKNWN,RECTAN)
C
C
C     PURPOSE
C         TO DETERMINE THE VALUE OF THE RESIDUAL AT
C         THE POINT (X,Y).
C
C     PARAMETERS
C
C         X,Y - THE POINT AT WHICH TO COMPUTE THE RESIDUAL
C         IPTCOM - RESULT FROM SUBROUTINE OUTSID
C               =1 POINT IS INTERIOR TO DOMAIN
C               =2 POINT IS ON BOUNDARY OF DOMAIN
C
C
      REAL UNKNWN(*),BVALUS(4),DERVSL(6),CVALUS(7)
C
C
      LOGICAL RECTAN
C
      COMMON /BRANZZ/BRANGE(2,1)
      COMMON /BNDRY/IPIECE,NBOUND,NBNDPT
      COMMON /GRIDXZ/GRIDX(1)
      COMMON /GRIDYZ/GRIDY(1)
      COMMON /NUMCON/EPSGRD
      COMMON /PROBR/AX,BX,AY,BY
C
      RESID = 0.
C
C         CHECK TO SEE IF (X,Y) IS ON THE BOUNDARY.
C         IN THIS CASE THE RESIDUAL OF THE BOUNDARY
C         OPERATOR MUST BE COMPUTED.
C
C         FOR RECTANGULAR DOMAINS, NEED TO SET IPTCOM
C
      IF ( .NOT. RECTAN) GO TO 10
      IPTCOM = 1
      IF (X.EQ.AX .OR. X.EQ.BX .OR. Y.EQ.AY .OR. Y.EQ.BY) IPTCOM = 2
C
   10 IF (IPTCOM.EQ.2) GO TO 20
C
C         EVALUATE THE COEFFICIENTS OF THE PDE OPERATOR.
C
      CALL PDE(X,Y,CVALUS)
      GO TO 60
C
C         EVALUATE COEFFICIENTS OF BOUNDARY OPERATOR
C
   20 CONTINUE
C
C         FIND THE BOUNDARY PIECE THE POINT IS ON
C
      IF (RECTAN) GO TO 40
C
C            NONRECTANGULAR CASE
C
      DO 30 IP = 1,NBOUND
         IPIECE = IP
         P = BRANGE(1,IP)
         CALL BCOORD(P,X1,Y1,IPIECE)
         P = BRANGE(2,IP)
         CALL BCOORD(P,X2,Y2,IPIECE)
         IF (X.LT.AMIN1(X1,X2)-EPSGRD) GO TO 30
         IF (Y.LT.AMIN1(Y1,Y2)-EPSGRD) GO TO 30
         IF (X.GT.AMAX1(X1,X2)+EPSGRD) GO TO 30
         IF (Y.GT.AMAX1(Y1,Y2)+EPSGRD) GO TO 30
         GO TO 50

   30 CONTINUE
      GO TO 80
C
C        RECTANGULAR CASE
C
   40 IF (X.EQ.BX) IP = 1
      IF (Y.EQ.AY) IP = 2
      IF (X.EQ.AX) IP = 3
      IF (Y.EQ.BY) IP = 4
C
   50 CONTINUE
C
C       EVALUATE BOUNDARY COEFFICIENTS
C
      TEMP = BCOND(IP,X,Y,BVALUS)
C
C         EVALUATE COMPUTED SOLUTION AND ITS DERIVATIVES AT (X,Y)
C
   60 CONTINUE
      TEMP = COLAPR(X,Y,UNKNWN,DERVSL,RECTAN)
C
C         DERIVATIVES OF SOLUTION ARE NOW IN THE ARRAY DERVSL
C
C         COMPUTE THE RESIDUAL FOR ALL CASES
C
      IF (IPTCOM.EQ.2) GO TO 70
C
C         RESIDUAL OF PDE OPERATOR IS
C
      RESID = CVALUS(1)*DERVSL(1) + CVALUS(2)*DERVSL(2) +
     .        CVALUS(3)*DERVSL(3) + CVALUS(4)*DERVSL(4) +
     .        CVALUS(5)*DERVSL(5) + CVALUS(6)*DERVSL(6) - PDERHS(X,Y)
      GO TO 80
C
C         THE RESIDUAL OF THE BOUNDARY OPERATOR IS
C
   70 CONTINUE
      RESID = BVALUS(1)*DERVSL(6) + BVALUS(2)*DERVSL(4) +
     .        BVALUS(3)*DERVSL(5) - BVALUS(4)
C
   80 CONTINUE
      RETURN
      END
C/////////////////////////////////////////////////////////////////////
C///////////////////// END OF LOGICAL FILE 3 /////////////////////////
C/////////////////////////////////////////////////////////////////////
C////////////////////////////////////////////////////////////////////
C///////////////// ALGORITHM GENCOL ////////////////////////////////
C///////////////////////////////////////////////////////////////////
C>>>>>>>>> LOGICAL FILE 4 : RELATED SUBPROGRAMS TO S/R REGION<<<<<<<<
C////////////////////////////////////////////////////////////////////
      SUBROUTINE DOMAIN(XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,BCOORD,
     .                  GTYPE,XBOUND,YBOUND,PIECE,BPTYPE,BNEIGH,BGRID,
     .                  BPARAM,NBDIM)
C
C   ***** THIS SUBROUTINE PROCESSES THE RECTANGULAR GRID AND SPECIFIED
C  BOUNDARY. IT APPLIES TO ONE CLOSED LOOP OR ARC OF THE BOUNDARY AND
C  MAY BE CALLED SEVERAL TIMES FOR A COMPLEX DOMAIN
C
C  ------------------- INPUT INFORMATION FOR DOMAIN ------------------
C     THE INPUT INFORMATION IS IN THE COMMON BLOCKS AND ARGUMENTS
C               SEE THE MAIN DRIVER REGION FOR MORE DETAILS
C
C 1** OUTPUT LEVEL CONTROL *******
C         LEVEL          = CONTROL SETTING, DETAILS GIVEN BELOW
C
C 2** GRID DEFINITION ************
C
C      NGDIMX,NGDIMY = NGRIDX,NGRIDY = GRID LINES IN X AND Y COORDINATES
C      XGRID(IX),YGRID(JY)  FOR  IX = 1 TO NGRIDX, JY = 1 TO NGRIDY
C
C 3** BOUNDARY DEFINITION ********
C
C         N1BND  = NUMBER OF THE FIRST BOUNDARY PIECE
C                  WILL DIFFER FROM 1 AFTER FIRST CALL
C         NPDIM  = ARRAY DIMENSION FOR BOUNDARY PIECES
C         N1BDPT = NUMBER OF THE FIRST BOUNDARY POINT
C                  WILL DIFFER FROM 1 AFTER FIRST CALL
C         NBDIM  = ARRAY DIMENSION FOR BOUNDARY POINTS
C         BCOORD = PARAMETERIZED DEFINITION OF THE BOUNDARY.
C                  BCOORD(P,X,Y,IPIECE) GIVES THE X,Y VALUES OF THE
C                  POINT ON PIECE IPIECE WITH PARAMETER VALUE = P.
C         BRANGE(2,I) = FIRST AND LAST VALUES OF PARAMETERS DEFINING
C                       THE I-TH BOUNDARY PIECE
C         CLOCKW = SWITCH TO SPECIFY BOUNDARY ORIENTATION
C                  .TRUE. MEANS BOUNDARY IS CLOCKWISE
C                  .FALSE.MEANS BOUNDARY IS COUNTER-CLOCKWISE
C         ARC    = .TRUE. MEANS DOMAIN IS AN ARC WITH NO INTERIOR
C
C  -----------     THE OUTPUT OF DOMAIN IS IN TWO PARTS -------------
C
C  1** GRID SPECIFICATION **********
C
C    GTYPE(IX,JY) FOR  IX = 1 TO NGRIDX, JY = 1 TO NGRIDY
C      THE VALUES IN THIS ARRAY GIVE THE TYPE OF THE GRID POINTS AND
C      INFORMATION ABOUT THEIR RELATION TO THE BOUNDARY.
C      THERE IS A PACKING FACTOR IPACKB WHICH IS NORMALLY 1000. FOR
C      VERY LARGE PROBLEMS, IPACKB AND RELATED CONSTANT NSIDE =
C      IPACKB - 1 MUST BE INCREASED SO THAT NSIDE .GT. NBNDPT.
C      POSSIBLE VALUE ARE( ASSUMING IPACKB = 1000)
C
C    = INTEGER OVER 1000
C      GRID POINT IS NEXT TO THE BOUNDARY AND THE GTYPE VALUE IS
C      GTYPE = K + 1000*J   WHERE
C            K IS THE INDEX OF THE LOWEST NUMBERED BOUNDARY NEIGHBOR
C                MUST DOUBLE CHECK USE WHEN K = 1
C            J = FOUR BITS TO NOTE LOCATION OF BOUNDARY POINTS
C                 0001 - BOUNDARY NEIGHBOR TO NORTH  (NOON)
C                 0010 - BOUNDARY NEIGHBOR TO EAST   (3 O'CLOCK)
C                 0100 - BOUNDARY NEIGHBOR TO SOUTH  (6 O'CLOCK)
C                 1000 - BOUNDARY NEIGHBOR TO WEST   (9 O'CLOCK)
C            THUS J=9 (1001 IN BINARY) IMPLIES THAT THERE ARE BOUNDARY
C            NEIGHBORS TO THE NORTH AND WEST
C          EXAMPLES ( X = BNDRY PT., 0 = GRID PT )
C
C                        X          X       X  0   X
C
C                     X  0          0  X       X
C                       J=9        J=3        J=14
C
C          ***** NOTE THAT GTYPE IS INITIALLY SET NEG. AND THEN MADE
C                POSITIVE WHEN THE INTERIOR IS FILLED
C
C    = 999
C      MEANS GRID POINT IS INTERIOR TO THE REGION AND NOT CLOSE
C      TO THE BOUNDARY, NSIDE HAS BEEN SET TO 999 SO THAT THE PACKING
C
C    = INTEGER LESS THAN 1000  ( 1000 IS IPACKB VALUE)
C      GRID POINT IS ALSO BOUNDARY PT., GTYPE IS ITS INDEX
C
C    = 0
C      GRID POINT IS EXTERIOR FAR FROM THE BOUNDARY
C
C    = NEGATIVE INTEGER
C      GRID POINT IS EXTERIOR NEXT TO THE BOUNDARY, ITS LOCATION
C      RELATIVE TO THE BOUNDARY IS ENCODED AS FOR INTERIOR POINTS
C      NEAR THE BOUNDARY
C
C  2** BOUNDARY SPECIFICATION **********
C
C         NBNDPT = NUMBER OF THE LAST BOUNDARY POINT ACTUALLY FOUND
C                     ON THIS CALL TO DOMAIN
C         NBNDP1 = TOTAL NUMBER OF BOUNDARY POINTS SO FAR.
C                  THIS INCLUDES FIRST POINT OF ANY CLOSED BOUNDARY
C                  COPIED TO END OF THAT BOUNDARY - BUT NOT ARCS
C         XBOUND(I),YBOUND(I) = COORDINATES OF I-TH BOUNDARY POINT
C         BPARAM(I) = PARAMETER VALUE OF I-TH BOUNDARY POINT
C         PIECE(I)  = INDEX OF BOUNDARY PIECE TO WHICH PT. BELONGS
C                        SMALLEST NUMBER FOR CORNER POINTS
C         BPTYPE(I) = TYPE OF BOUNDARY POINT
C                   = HORZ,VERT,BOTH,INTE OR JUMP
C         BNEIGH(I) = POINTER TO THE INTERIOR POINTS FROM THE I-TH
C                     BOUNDARY POINT. SAME SCHEME IS USED TO ENCODE
C                     DIRECTIONS AS FOR THE J PART OF GTYPE ABOVE
C         BGRID(I) = IX + IPACKB*JY WHEN PT. I IS IN GRID SQUARE IX,JY
C
C         MAXIMUM NUMBER OF BOUNDARY POINTS = NBDIM
C         THIS IS ESTIMATED BY THE CALLING PROGRAM FOR DIMENSIONING
C         THE VARIOUS ARRAYS.
C
C ***********   DOMAIN PROCESSOR SUBPROGRAMS   ***********************
C
C   REGION    - MAIN DRIVER AND USER INTERFACE
C   HOLE      - ALTERNATE DRIVER, INSERTS HOLES IN FIRST DOMAIN
C      REMOVH - REMOVES HOLE FROM GRID TYPES, UPDATES ALL INFO
C   DOMAIN    - MAIN PROGRAM TO PROCESS A DOMAIN
C   BWALK     - WALK ALONG BOUNDARY TO FIND GRID INTERSECTION
C      CHKTAN - CHECK FOR BOUNDARY TANGENT TO A GRID LINE
C      CROSS2 - CHECK FOR DOUBLE CROSSING OF A GRID LINE
C      DBACK  - FOLLOW BOUNDARY FOR A DOUBLE CROSSING OF A GRID LINE
C      REGULA - MODIFIED REGULA FALSI METHOD
C      SECANT - SECANT METHOD
C   CHANGE    - MAKE CHANGE OF BOUNDARY PIECE
C   FILL      - LOCATE AND FILL INTERIOR OF THE DOMAIN
C      EXPAND - EXPAND INTERIOR OF DOMAIN BY 1 POINT
C   GVALUS    - SET TYPES OF ALL GRID POINTS
C      ISETGT - COMPUTE GTYPE VALUE FOR A GRID POINT, SET IT
C   INSIDE    - UTILITY: CHECKS POINT INSIDE SPECIFIED SUBGRID
C   LOCATE    - UTILITY: LOCATE POINT IN GRID, TYPE IT
C   NEIGH     - COMPUTE POINTERS FROM BOUNDARY TO GRID POINTS
C   TABLGT    - TABLE THE GTYPE VALUES FOR THE GRID
C
C ********************************************************************
C *************  OUTPUT CONTROL  *************************************
C
C    LEVEL = 0   ***** FATAL ERROR MESSAGES **************************
C     BWALK  - BOUNDARY GOES OUTSIDE DOMAIN
C             - UNABLE TO FIND GRID WHERE BOUNDARY GOES
C             - UNABLE TO FIND GRID INTERSECTION WITH BOUNDARY
C     CHANGE - BOUNDARY PIECES DO NOT JOIN UP
C     DOMAIN - BOUNDARY PARAMETER NOT INCREASING
C             - OVERFLOW OF STORAGE ALLOCATED FOR BOUNDARY POINTS
C                 MUST CHANGE DECLARATIONS IN PROGRAM CALLING REGION
C             - ABNORMAL EXIT FROM BOUNDARY PROCESSING
C                 PROBABLY CANNOT HAPPEN
C             - BOUNDARY DOES NOT CLOSE
C     FILL   - FAILURE TO FIND INTERIOR OF DOMAIN
C             - OVERFLOW IN QUEUE FOR FILLING INTERIOR OF DOMAIN
C                 MUST INCREASE QLIMIT AND RECOMPILE SUBROUTINE FILL
C     GVALUS - HAVE ILLEGAL POINT TYPE. THINGS ARE REALLY MESSED UP
C     LOCATE - ASKED TO FIND POINT OUTSIDE DOMAIN
C                 INPUT ERROR OR PROGRAM BUG
C     REMOVH - HOLE IS TOO CLOSE TO BOUNDARY OF CONTAINING DOMAIN
C                 NEED TO USE FINER GRID
C     REGION - GRID LINES ARE NOT INCREASING, INPUT ERROR
C             - LESS THAN 2 GRID LINES IN X OR Y DIRECTION, INPUT ERROR
C
C   LEVEL = 1   *****  MINIMAL MESSAGES *******************************
C
C     REGION - NEW PAGE, DOMAIN PROCESSING STARTS
C     HOLE   - NEW PAGE, HOLE PROCESSING STARTS
C     DOMAIN - NUMBER OF BOUNDARY/GRID INTERSECTIONS
C             - FATAL ERROR NOTE( IF PRESENT )
C
C   LEVEL = 2   ***** SUMMARY TRACE **********************************
C
C     CHANGE - CHANGE OF BOUNDARY PIECE
C     CHKTAN - POINT REPLACEMENT DUE TO TANGENCY TO GRID LINE
C     DOMAIN - INITIAL BOUNDARY POINT
C             - BOUNDARY POINT FOUND
C             - TABLE OF GRID POINT TYPES( FROM TABLGT )
C             - SUMMARY OUTPUT OF BOUNDARY POINT ARRAYS
C     FILL   - WARNING, FIRST TRY TO LOCATE INTERIOR FAILS
C             - FAILURE INFORMATION ( IF NEEDED )
C     REMOVH - TABLE OF GRID POINT TYPES ( FROM TABLGT )
C
C   LEVEL = 3   **** SOME DETAILS OF THE PROGRAM OPERATION ***********
C
C     BWALK  - START
C             - WARNING, POINT FOUND IS NOT IN EXPECTED GRID
C             - BOUNDARY IS VERTICAL/HORIZONTAL
C             - FINISH
C     CHANGE - BOUNDARY PIECE CHANGE DATA
C     CHKTAN - DATA FOR THE TANGENCY TEST( WHEN IT SUCCEEDS )
C     CROSS2 - SIGNS FOR DOUBLE CROSSING CHECK
C             - SUCCESSFUL RESULTS
C     DOMAIN - BOUNDARY PIECE CHANGE DATA
C             - HEADING FOR TYPING OF GRID POINTS
C     FILL   - DATA ON SEARCH FOR INTERIOR POINTS
C     NEIGH  - SUMMARY OF EXECUTION
C
C   LEVEL = 4   **** MORE DETAILS ABOUT OPERATION ********************
C
C     BWALK  - EXPANSION OF PARAMETER STEP
C     DBACK  - SUMMARY OF EXECUTION
C     EXPAND - DETAILS OF IDENTIFICATION OF THE INTERIOR POINTS
C     GVALUS - DATA USED FOR TYPING GRID POINTS
C     LOCATE - DATA FOR POINTS LOCATED
C     NEIGH  - DATA FOR NEIGHBORS LOCATED
C     REGULA - SUMMARY OF EXECUTION
C     SECANT - START INFORMATION
C             - FINISH INFORMATION
C
C   LEVEL = 5   **** DETAILS OF MANY LOOPS ,ETC. ********************
C
C     BWALK  - DATA ON INITIAL ROUGH GUESSES AT INTERSECTIONS
C     DBACK  - LOOP DETAILS
C     FILL   - DATA FOR LOOP TO FIND FIRST INTERIOR POINT
C             - DATA ON CLASSIFYING INTERIOR POINTS
C     GVALUS - DATA ON TYPING GRID POINTS
C     SECANT - BOUNDS USED ON VARIABLES
C             - DATA WHEN VARIABLE BOUNDS ARE USED
C************************************************* END OF OUTPUT SPECS
C
      EXTERNAL BCOORD

      CHARACTER *4 TYPEC,BPTYPE
      REAL XGRID(NGDIMX),YGRID(NGDIMY),BRANGE(2,NPDIM)
      DIMENSION GTYPE(NGDIMX,NGDIMY),XBOUND(NBDIM),YBOUND(NBDIM),
     .          BPTYPE(NBDIM),BNEIGH(NBDIM),BGRID(NBDIM),BPARAM(NBDIM),
     .          PIECE(NBDIM)
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
C --------------------------  PROCESS A DOMAIN  ------------------------
C
C        CHECK THAT THE BOUNDARY PARAMETERS ARE ALWAYS INCREASING
      FATAL = .FALSE.
      DO 10 IB = N1BND,NBOUND
         IF (BRANGE(2,IB).GT.BRANGE(1,IB)) GO TO 10
C             FATAL ERROR IN BOUNDARY PARAMETERIZATION
C
         IF (LEVEL.GE.0) WRITE (MOUTPT,9001) IB
         FATAL = .TRUE.
   10 CONTINUE
      IF (FATAL) GO TO 190
      IF (LEVEL.GE.4) WRITE (MOUTPT,9011) CLOCKW,ARC,INHOLE,NGRIDX,
     .    NGRIDY,N1BND,N1BDPT,NBOUND,IPACKB,NSIDE,MOUTPT,QLIMIT,EPSGRD,
     .    EPSTAN
C         LOCATE FIRST BOUNDARY POINT
      IPIECE = N1BND
      PARAM = BRANGE(1,IPIECE)
      CALL BCOORD(PARAM,XB,YB,IPIECE)
C
      CALL LOCATE(XB,YB,IX,JY,TYPE,XGRID,YGRID,NGDIMX,NGDIMY)
      IF (FATAL) GO TO 190
C               SET INFORMATION FOR FIRST BOUNDARY POINT
      NBNDPT = N1BDPT
      XBOUND(N1BDPT) = XB
      YBOUND(N1BDPT) = YB
      PIECE(N1BDPT) = N1BND
      BPTYPE(N1BDPT) = TYPE
      BGRID(N1BDPT) = IX + IPACKB*JY
      BPARAM(N1BDPT) = PARAM
      IF ((TYPE.EQ.VERT) .OR. (TYPE.EQ.BOTH)) XBOUND(N1BDPT) = XGRID(IX)
      IF ((TYPE.EQ.HORZ) .OR. (TYPE.EQ.BOTH)) YBOUND(N1BDPT) = YGRID(JY)
C         DEBUG
      IF (LEVEL.GE.2) WRITE (MOUTPT,9021) IPIECE,PARAM,XB,YB,TYPE
C
C ----------------------- START LOOP OVER BOUNDARY PIECES  ----------
C         FOLLOW EACH PIECE OF THE BOUNDARY AND DETERMINE ALL GRID
C         LINE INTERSECTIONS AND VARIOUS INFORMATION ASSOCIATED WITH
C         THE INTERSECTION POINTS( = BOUNDARY POINTS)
      DO 50 IB = N1BND,NBOUND
C
C             CHECK FOR HITTING END OF BOUNDARY IMMEDIATELY
         IF (PARAM.GT.BRANGE(2,IPIECE)-EPSGRD) GO TO 40
C
C -----------  LOOP WHICH WALKS ALONG THE IB-TH BOUNDARY PIECE  ----
         DO 30 IDUMB = 1,NBDIM
C         LOCATE THE NEXT GRID LINE INTERSECTION BY WALKING FROM XB,YB
C         ALONG THE BOUNDARY TO GRID LINES IXNU,JYNU AND PARAM = PNEXT
            CALL BWALK(XB,YB,IX,JY,PNEXT,XNEXT,YNEXT,IXNU,JYNU,XGRID,
     .                 YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,BCOORD)
            IF (FATAL) GO TO 190
C
C                  CHECK FOR HAVING REACHED THE END OF THE BOUNDARY
C     *** EXIT *** THIS IS THE NORMAL EXIT OF LOOP OVER BOUNDARY PIECE
            IF (IB.EQ.NBOUND .AND. PNEXT.GT.BRANGE(2,IB)-
     .          EPSGRD) GO TO 60
C
C             CHECK FOR GOING PAST THE END OF THIS BOUNDARY PIECES
            IF (PNEXT.GE.BRANGE(2,IPIECE)-.5E0*EPSGRD) GO TO 40
C
C                   PUT NEW POINT IN BOUNDARY POINT LIST
            NBNDPT = NBNDPT + 1
            IF (NBNDPT.LE.NBDIM) GO TO 20
C
C             HAVE EXCEEDED BOUNDARY POINT LIMIT
C
            IF (LEVEL.GE.0) WRITE (MOUTPT,9031) NBDIM
            FATAL = .TRUE.
            GO TO 190

   20       CONTINUE
            XBOUND(NBNDPT) = XNEXT
            YBOUND(NBNDPT) = YNEXT
            PIECE(NBNDPT) = IPIECE
            BPTYPE(NBNDPT) = TYPE
            BGRID(NBNDPT) = IXNU + IPACKB*JYNU
            BPARAM(NBNDPT) = PNEXT
            IX = IXNU
            JY = JYNU
            PARAM = PNEXT
            IF (LEVEL.GE.2) WRITE (MOUTPT,9041) NBNDPT,XBOUND(NBNDPT),
     .          YBOUND(NBNDPT)
            IF (LEVEL.GE.3) WRITE (MOUTPT,9051) IPIECE,TYPE,IX,JY
C
C         CHECK TO SEE IF THIS IS A POINT OF TANGENCY.
C         IF SO, REPLACE THE PREVIOUS POINT IN THE BOUNDARY LIST
C                BY THIS ONE IF IT IS VERY CLOSE BY
            CALL CHKTAN(XBOUND,YBOUND,BPARAM,BPTYPE,BGRID,NBDIM,NPDIM,
     .                  BRANGE,XGRID,YGRID,NGDIMX,NGDIMY,TYPE,IXNU,JYNU)
C               TANGENT CHECK MIGHT MOVE START OF NEXT SEARCH AWAY
C               FROM XNEXT,YNEXT.  THE VALUE OF PARAM IS ALWAYS
C               INCREASING AND IS USED TO FORCE SEARCH AWAY FROM THE
C               POINT OF TEANGENCY.
            CALL BCOORD(PARAM,XB,YB,IPIECE)
   30    CONTINUE
C
C     ---------------- END OF LOOP OVER POINTS ON 1 BOUNDARY PIECE ---
C
   40    CONTINUE
C
C             HAVE COME TO END OF PIECE, TAKE CARE OF LAST POINT
C             THIS ALSO INSERTS THE END POINT IN THE BOUNDARY LIST,
C             LOCATES THE NEXT POINT ON THE BOUNDARY AND PUTS IT IN
C             THE LIST TOO.
C
C                  DO NOT TRY TO CHANGE THE LAST PIECE
         IF (IB.EQ.NBOUND) GO TO 60
C
         CALL CHANGE(BRANGE,NPDIM,XGRID,YGRID,NGDIMX,NGDIMY,NBDIM,
     .               XBOUND,YBOUND,PIECE,BPTYPE,BGRID,BPARAM,BCOORD)
         IF (FATAL) GO TO 190
         CALL LOCATE(XBOUND(NBNDPT),YBOUND(NBNDPT),IX,JY,TYPE,XGRID,
     .               YGRID,NGDIMX,NGDIMY)
         CALL BCOORD(PARAM,XB,YB,IPIECE)
   50 CONTINUE
C
C                ABNORMAL EXIT FROM BOUNDARY PROCESSING
C           THIS CANNOT HAPPEN UNLESS THINGS ARE REALLY MESSED UP
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9061) PARAM,IPIECE,NBOUND,NBNDPT,
     .    IX,JY,NGRIDX,NGRIDY,NPDIM,NBDIM
      FATAL = .TRUE.
      GO TO 190
C
   60 CONTINUE
C             CHECK CLOSING OF THE BOUNDARY
      IF (ARC) GO TO 70
      CALL BCOORD(BRANGE(2,NBOUND),XLAST,YLAST,IPIECE)
      DIST = ABS(XLAST-XBOUND(N1BDPT)) + ABS(YLAST-YBOUND(N1BDPT))
      IF (DIST.LT.3.*EPSGRD) GO TO 70
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9071)
      FATAL = .TRUE.
      GO TO 190
C
C     -------------- END LOOP OVER BOUNDARY PIECES ---------------------
   70 CONTINUE
C
C               FIX END OF BOUNDARY, DUPLICATE FIRST POINT AT END
C               OF LIST, CHECK FOR TANGENCY
C         IF NOT ON AN ARC
      IF (ARC) GO TO 100
C
C         COPY FIRST BOUNDARY POINT IN BOUNDARY LIST AS LAST POINT
C                  ONLY IF THERE IS SPACE FOR IT
      IF (NBNDPT+1.LT.NBDIM) GO TO 80
C
      WRITE (MOUTPT,9031) NBDIM
      FATAL = .TRUE.
      RETURN
C
C           PRINT WARNING IF LAST PIECE HAS ONLY ONE BOUNDARY POINT
   80 IF (PIECE(NBNDPT).EQ.NBOUND) GO TO 90
      IF (LEVEL.GE.1) WRITE (MOUTPT,9081)
   90 XBOUND(NBNDPT+1) = XBOUND(N1BDPT)
      YBOUND(NBNDPT+1) = YBOUND(N1BDPT)
      PIECE(NBNDPT+1) = NBOUND
      BPTYPE(NBNDPT+1) = BPTYPE(N1BDPT)
      BGRID(NBNDPT+1) = BGRID(N1BDPT)
      BPARAM(NBNDPT+1) = BRANGE(2,NBOUND)
      NBNDP1 = NBNDPT + 1
      GO TO 110
C
C            ON AN ARC PUT IN THE LAST BOUNDARY POINT FOUND
  100 CONTINUE
      NBNDPT = NBNDPT + 1
      XBOUND(NBNDPT) = XNEXT
      YBOUND(NBNDPT) = YNEXT
      PIECE(NBNDPT) = IPIECE
      BPTYPE(NBNDPT) = TYPE
      BGRID(NBNDPT) = IXNU + IPACKB*JYNU
      BPARAM(NBNDPT) = PNEXT
      NBNDP1 = NBNDPT
C
C
C            CHECK FOR TANGENCY AS THE BOUNDARY CLOSES
C            TEMPORARILY INCREASE BOUNDARY POINT COUNT
  110 NBNDPT = NBNDPT + 1
      CALL CHKTAN(XBOUND,YBOUND,BPARAM,BPTYPE,BGRID,NBDIM,NPDIM,BRANGE,
     .            XGRID,YGRID,NGDIMX,NGDIMY,TYPE,IXNU,JYNU)
      NBNDPT = NBNDPT - 1
C            NO DUPLICATE AT THE END FOR AN ARC
      IF (ARC) NBNDPT = NBNDP1
C
      IF (LEVEL.GE.3) WRITE (MOUTPT,9091)
C
C                   FIND BOUNDARY NEIGHBORS GRID POINT INFORMATION
C        THIS SETS VALUES OF GTYPE FOR GRID POINTS NEXT TO THE BOUNDARY
C        EXCEPT THAT THE SIGN( = INTERIOR/EXTERIOR FLAG ) IS NOT
C        SET UNTIL FILL IS CALLED LATER.
C                   SET ALL GRID POINT TYPES = 0
      DO 130 IX = 1,NGRIDX
         DO 120 JY = 1,NGRIDY
            GTYPE(IX,JY) = 0
  120    CONTINUE
  130 CONTINUE
      CALL GVALUS(XGRID,YGRID,NGDIMX,NGDIMY,GTYPE,NBDIM,XBOUND,YBOUND,
     .            BPTYPE,BGRID)
C
C      NOW FILL IN THE TYPES OF GRID POINTS AWAY FROM THE BOUNDARY
C      IN THE INTERIOR OF THE REGION.  SET GTYPE POS. NEXT TO BOUNDARY
C
      IF ( .NOT. ARC) CALL FILL(GTYPE,NGDIMX,NGDIMY,NBDIM,BGRID,BPTYPE,
     .                          PIECE,XBOUND,YBOUND)
      IF (FATAL) GO TO 190
C
C          GO ALONG THE BOUNDARY AND SET BNEIGH VALUES FOR ALL
C          THE BOUNDARY POINTS.
      CALL NEIGH(NGDIMX,NGDIMY,GTYPE,NBDIM,BNEIGH,BPTYPE,BGRID)
C          SET BNEIGH FOR LAST+1 POINT IF NOT ON AN ARC
      IF ( .NOT. ARC) BNEIGH(NBNDP1) = BNEIGH(NBNDPT)
C
C                       SUMMARY OUTPUT
      NBPTF = NBNDPT - N1BDPT + 1
      NBDF = NBOUND - N1BND + 1
      IF (LEVEL.GE.1) WRITE (MOUTPT,9101) NBPTF,NBDF,NGRIDX,NGRIDY
      IF (LEVEL.LT.2) GO TO 150
      WRITE (MOUTPT,9111)
      DO 140 IB = N1BDPT,NBNDP1
         TYPEC = BPTYPE(IB)
         WRITE (MOUTPT,9121) IB,XBOUND(IB),YBOUND(IB),BPARAM(IB),
     .     PIECE(IB),TYPEC,BGRID(IB),BNEIGH(IB)
  140 CONTINUE
  150 CONTINUE
C                 PRINT TABLE OF GTYPE VALUES
      IF (LEVEL.GT.1) CALL TABLGT(MOUTPT,1,NGRIDX,NGRIDY,GTYPE)
C
C       CHECK FOR INTERIOR POINTS ON BOUNDARY OF GTYPE ARRAY
C       IF FOUND IT MEANS THERE IS SOMETHING WRONG
      DO 160 I = 1,NGDIMX
         IF (GTYPE(I,1).EQ.NSIDE .OR. GTYPE(I,NGDIMY).EQ.
     .       NSIDE) GO TO 180
  160 CONTINUE
      DO 170 J = 1,NGDIMY
         IF (GTYPE(1,J).EQ.NSIDE .OR. GTYPE(NGDIMX,J).EQ.
     .       NSIDE) GO TO 180
  170 CONTINUE
C
C         ******* NORMAL RETURN ******** NO PROBLEMS DETECTED
      RETURN
C
C                       SERIOUS PROBLEM DETECTED
  180 CONTINUE
C
      WRITE (MOUTPT,9131)
      FATAL = .TRUE.
C
C                    ALL FATAL ERRORS COME TO 300
  190 CONTINUE
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9141)
      IF (LEVEL.GE.2 .AND. N1BDPT.LT.NBNDPT) WRITE (MOUTPT,9101) NBPTF,
     .    NBDF,NGRIDX,NGRIDY
      IF (LEVEL.LT.2 .OR. N1BDPT.GE.NBNDPT) GO TO 210
      DO 200 IB = N1BDPT,NBNDPT
         TYPEC = BPTYPE(IB)
         WRITE (MOUTPT,9151) IB,XBOUND(IB),YBOUND(IB),BPARAM(IB),
     .     PIECE(IB),TYPEC,BGRID(IB)
  200 CONTINUE
  210 CONTINUE
      FATAL = .TRUE.
      RETURN

 9001 FORMAT (/5 (3H **),' FATAL ERROR IN BOUNDARY DEFINITION',/9X,
     .                 'PARAMETER OF PIECE ',I4,' IS DECREASING')
 9011 FORMAT (/20X,'DATA AND CONSTANTS FOR THIS DOMAIN PROCESSING'/3X,
     .  'SWITCHES CLOCKW, ARC AND INHOLE =',3L4/3X,'NGRIDX,NGRIDY =',
     .  2I4,' FIRST BOUNDARY POINT, PIECE =',2I4,' WITH ',I3,
     .  ' BOUNDARY PIECES'/3X,'CONSTANTS IPACKB,NSIDE,MOUTPT,QLIMIT = ',
     .  4I6/3X,'GEOMETRY TOLERANCES EPSGRD, EPSTAN =',2E15.5)
 9021 FORMAT (5X,'ON PIECE',I4,' PARAM ',F10.6,', GIVES COORDS ',2F10.6,
     .  ', TYPE ',A4)
 9031 FORMAT (/5 (3H **),' FATAL ERROR IN PROCESSING DOMAIN, ',
     .                 'ESTIMATED NUMBER',I4/9X,
     .                'OF BOUNDARY-GRID INTERSECTION POINTS IS TOO LOW,'
     .                 /9X,'THE REMEDY IS TO INCREASE NBDIM')
 9041 FORMAT (2X,3 ('--'),' FOUND BOUNDARY POINT ',I3,
     .                  ' WITH COORDINATES ',2F10.6)
 9051 FORMAT (15X,'ON PIECE ',I2,' OF TYPE ',A4,' IN GRID ',2I3)
 9061 FORMAT (/5 (3H **),' ABNORMAL EXIT FROM LOOP IN BOUNDARY ',
     .                 'PROCESSING, = FATAL ERROR'/9X,
     .               'HAPPENS ONLY IF BOUNDARY INFO IS REALLY MESSED UP'
     .                 /6X,' PARAM =',F12.7/6X,
     .     58H IPIECE,NBOUND,NBNDPT,IX,JY,NGRIDX,NGRIDY,NPDIM,NBDIM =
     .                 /6X,9I5)
 9071 FORMAT (/5 (3H **),' FATAL ERROR, BOUNDARY DOES NOT CLOSE')
 9081 FORMAT (/5 (3H **),' WARNING, LAST PIECE HAS ONLY 1 BOUNDARY',
     .                 ' POINT')
 9091 FORMAT (//20X,'SET GTYPE FOR NEIGHBORS OF BOUNDARY POINTS')
 9101 FORMAT (/5X,'BOUNDARY POINTS FOUND',I12/5X,
     .  'BOUNDARY PIECES FOUND',I12/5X,'GRID SIZE',I17,' BY',I4/5X,
     .  'EXECUTION SUCCESSFUL')
 9111 FORMAT (/3X,'BOUNDARY PT  XBOUND',8X,'YBOUND',8X,'BPARAM',
     .  '    PIECE BPTYPE  BGRID  BNEIGH')
 9121 FORMAT (I10,3F14.9,I5,4X,A4,I8,I6)
 9131 FORMAT (4X,'FOUND INTERIOR POINTS ON EDGE OF GRID'/4X,
     .  'PROBABLE CAUSES ARE:'/8X,'BOUNDARY ORIENTATION IS WRONG'/8X,
     .  'BOUNDARY STARTS WHERE INTERIOR IS TOO THIN'/8X,
     .  'BOUNDARY OSCILLATES TOO RAPIDLY SOMEPLACE')
 9141 FORMAT (/5 (3H **),' FATAL ERROR IN DOMAIN,'/9X,
     .                 'STOPS PROCESSING RUN')
 9151 FORMAT (I10,3F14.9,I5,4X,A4,I8)
      END
      SUBROUTINE BWALK(XSTART,YSTART,IX,JY,PNEW,XNEW,YNEW,INEW,JNEW,
     .                 XGRID,YGRID,NGDIMX,NGDIMY,BRANGE,NPDIM,BCOORD)
C
C     THIS SUBROUTINE STARTS AT (XSTART,YSTART) AND WALKS ALONG THE
C  BOUNDARY TO THE NEXT GRID LINE INTERSECTION AT THE POINT
C  (XNEW,YNEW) WITH PARAMETER PNEW IN GRID SQUARE INEW,JNEW.
C  THE PARAMETER START IS PARAM IN COMMON
C  THE INTERSECTION POINT TYPE IS FOUND AND PLACED IN COMMON AS TYPE
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
      EXTERNAL BCOORD

      LOGICAL XCASE,INSIDE,ALLIN
      REAL XGRID(NGDIMX),YGRID(NGDIMY),BRANGE(2,NPDIM),GESS1(4),
     .     GESS2(4),TSIGN(4),SIGNOW(4),X(4),Y(4),P(4)
      INTEGER LTRY(4)
      CHARACTER *4 TYPEC,OUTCM(4)
      CHARACTER *4 WBACK,NOTRY,SUCCES,FAIL,OUTSID
      DATA WBACK,NOTRY,SUCCES,FAIL,OUTSID/'BACK','NOTR','SUCC','FAIL',
     .     'OUTS'/
C
C         **********  VARIABLES IN BWALK  **********
C
C  *** INPUT  *  XSTART,YSTART,IX,JY
C             VIA COMMON = PARAM,EPSGRD,EPSTAN,IPIECE
C             THE POINT WHERE THE WALK STARTS IS (XSTART,YSTART)
C             IN GRID SQUARE (IX,JY).
C
C  *** OUTPUT *  PNEW = PARAM VALVE FOR NEXT GRID LINE INTERSECTION
C             INEW, JNEW = INDICES FOR NEXT X,Y GRID LINES
C             XNEW,YNEW  = COORDINATES FOR NEXT INTERSECTION
C             TYPE(VIA COMMON) = TYPE OF INTERSECTION
C                              = HORZ, VERT OR BOTH
C             FATAL(VIA COMMON)= SWITCH FOR FATAL ERRORS
C
C  *** MAIN VARIABLES ***
C
C         DELP,DELPNU        = CHANGES TRIED FOR PARAM
C         DXGRD,DYGRD        = X AND Y GRID SPACING
C         XTAR1,XTAR2,YTAR1,YTAR2 = TARGETS OF SECANT METHOD
C                            = COORDINATES OF GRID SQUARE
C         XTEST,YTEST,PTEST  = TEST POINT INSIDE TARGET GRID SQUARE
C         X(I), 1 TO 4       = X COORDINATES OF SECANT RESULTS
C         Y(I)  1 TO 4       = Y COORDINATES OF SECANT RESULTS
C         P(I)  1 TO 4       = PARAMETERS FOR SECANT RESULTS
C         OUTCM(I)           = FLAGS FOR SECANT RESULTS
C                              NOTRY - NOT TRIED
C                              WBACK - WALKED BACK TO START
C                              SUCCES- HIT THE TARGET
C                              FAIL  - FAILED TO HIT TARGET
C                              OUTSID- HIT TARGET, BUT OUTSIDE GRID
C         PEND,XEND,YEND     = VALUES OF PARAM,X,Y AT END OF BOUNDARY
C         PSTOP              = PARAMETER VALUE KNOWN TO BE BEYOND
C                              THE NEXT INTERSECTION WITH THE GRID
C         INEXT, JNEXT       = INDEXES OF OTHER TWO SIDES OF SQUARE
C                               BOUNDARY IS PASSING THROUGH
C
C  *** OTHER VARIABLES ***
C
C        NHORZ,NVERT         = LOOP COUNTERS
C        LTRY                = ORDER FOR TRIES TO HIT TARGETS
C        TESTX,TESTY         = TEMP VALUES OF X,Y
C        TSIGN,SIGNOW        = SIGNS OF ERROR FOR THE 4 TARGETS
C        KTAR,STEP           = COUNTER AND STEP FOR INITIAL SEARCH
C        DGMAX,DIST          = DISTANCES FROM START
C        ALLIN               = .TRUE. IF ALL PTS MET ARE INSIDE GRID
C
      PEND = BRANGE(2,IPIECE)
      CALL BCOORD(PEND,XEND,YEND,IPIECE)
C          SET INITIAL OUTCOME SWITCHES, VARIABLES TO INVALID VALUES
      DO 10 I = 1,4
         X(I) = -100.E0
         Y(I) = -100.E0
         P(I) = PEND + 1.E0
         OUTCM(I) = NOTRY
   10 CONTINUE
      ALLIN = .TRUE.
C
C               SET INDICES OF GRIDS THAT MUST CONTAIN NEXT POINT
      ISOUTH = JY
      IWEST = IX
      INORTH = MIN0(JY+1,NGRIDY)
      IEAST = MIN0(IX+1,NGRIDX)
      IF (TYPE.EQ.VERT .OR. TYPE.EQ.BOTH) IWEST = MAX0(1,IX-1)
      IF (TYPE.EQ.HORZ .OR. TYPE.EQ.BOTH) ISOUTH = MAX0(1,JY-1)
      DXGRD = XGRID(IEAST) - XGRID(IWEST)
      DYGRD = YGRID(INORTH) - YGRID(ISOUTH)
C
C            FIND PARAMETER CHANGE = PTEST TO GIVE POINT IN THIS GRID
      FACTOR = 10.E0
      PTEST = (PEND-PARAM)/FACTOR
      DO 30 I1 = 1,2
         DO 20 I2 = 1,8
            CALL BCOORD(PARAM+PTEST,XTEST,YTEST,IPIECE)
            IF (INSIDE(XTEST,YTEST,IWEST,ISOUTH,IEAST,INORTH,XGRID,
     .          NGDIMX,YGRID,NGDIMY)) GO TO 50
C               NOT INSIDE REQUIRED GRID AREA YET
C               CHECK FOR GOING OUTSIDE GRID DOMAIN
            IF ( .NOT. INSIDE(XTEST,YTEST,1,1,NGRIDX,NGRIDY,XGRID,
     .          NGDIMX,YGRID,NGDIMY)) GO TO 40
            ALLIN = .FALSE.
            PSTOP = PARAM + PTEST
            PTEST = PTEST/FACTOR
   20    CONTINUE
         FACTOR = 2.E0*FACTOR
   30 CONTINUE
C
C                     ****** FATAL ERROR *****
C       DID NOT FIND REQUIRED POINT = INEXPLICABLE SITUATION
C       PARAMETER CHANGE IS NOW LESS THAN (PEND-PARAM)/1E+12
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9001) XSTART,YSTART
      FATAL = .TRUE.
      GO TO 520
C
C            ***** FATAL ERROR, BOUNDARY GOES OUTSIDE GRID *****
   40 CONTINUE
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9011) XSTART,YSTART
      PTEST = PARAM + PTEST
      IF (LEVEL.GE.2) WRITE (MOUTPT,9021) XTEST,YTEST,PTEST
      FATAL = .TRUE.
      GO TO 520
C
C               NOW CHECK TO SEE IF TEST PT. IS TOO CLOSE TO START
   50 DELP2 = PTEST
      FACTOR = 1.6E0
      DO 60 I1 = 1,6
         DX = ABS(XSTART-XTEST)
         DY = ABS(YSTART-YTEST)
         IF (AMIN1(DX,DY).GT.4.E0*EPSGRD) GO TO 80
C                       DO NOT LET DELP2 GET OUT OF RANGE
         IF (DELP2*FACTOR.GT.PEND-PARAM) GO TO 80
         DELP2 = DELP2*FACTOR
         CALL BCOORD(PARAM+DELP2,XTEST,YTEST,IPIECE)
         IF ( .NOT. INSIDE(XTEST,YTEST,IWEST,ISOUTH,IEAST,INORTH,XGRID,
     .       NGDIMX,YGRID,NGDIMY)) GO TO 70
   60 CONTINUE
C          DID NOT MOVE AWAY FROM START IN BOTH COORDINATES
C          COULD HAVE VERT OR HORZ BOUNDARY WHICH IS DETECTED LATER
      GO TO 80
C          MOVED TEST POINT OUTSIDE GRID AREA, RESTORE PREVIOUS DELP
   70 DELP2 = DELP2/FACTOR
      ALLIN = .FALSE.
   80 PTEST = DELP2
C               LOCATE THE POINT INSIDE THE DESIRED GRID AREA
      CALL BCOORD(PARAM+PTEST,XTEST,YTEST,IPIECE)
      CALL LOCATE(XTEST,YTEST,IX,JY,TYPEC,XGRID,YGRID,NGDIMX,NGDIMY)
C          CHECK LOCATION FAILURE - SHOULD NEVER HAPPEN
      IF ( .NOT. FATAL) GO TO 90
C
      WRITE (MOUTPT,9031) XSTART,YSTART
      FATAL = .TRUE.
      GO TO 520

   90 INEXT = MIN0(IX+1,NGRIDX)
      JNEXT = MIN0(JY+1,NGRIDY)
C
C             HAVE THE GRID SQUARE LOCATED, SET TARGETS
      XTAR1 = XGRID(INEXT)
      XTAR2 = XGRID(IX)
      YTAR1 = YGRID(JNEXT)
      YTAR2 = YGRID(JY)
C             COMPUTE THE CORRECT GRID WIDTHS
      DXGRD = ABS(XGRID(IX)-XGRID(INEXT))
      DYGRD = ABS(YGRID(JY)-YGRID(JNEXT))
C
      IF (LEVEL.GE.3) WRITE (MOUTPT,9041) XSTART,YSTART,IX,JY,PTEST,
     .    XTEST,YTEST,PARAM,XTAR1,XTAR2,YTAR1,YTAR2
C
C         STEP ALONG BOUNDARY TO FIND WHERE IT GOES AND TO OBTAIN ROUGH
C         ESTIMATES OF INTERSECTIONS. ORDER THE TARGETS WITH MOST
C         LIKELY FIRST AND OBTAIN INITIAL GUESSES FOR EACH
      DGMAX = SQRT(DXGRD**2+DYGRD**2)
C                SET INITIAL SIGNS OF DIRECTION TO TARGETS
      TSIGN(1) = XTAR1 - XSTART
      TSIGN(2) = XTAR2 - XSTART
      TSIGN(3) = YTAR1 - YSTART
      TSIGN(4) = YTAR2 - YSTART
      PSTOP = PEND
C
C **************** RETURN HERE FOR POSSIBLE ADDITIONAL ATTEMPTS
C **************** TO BRACKET THE TARGET COORDINATES
C                  LOCAL  RETURN FROM INSIDE DO-LOOP 22
C                  GLOBAL RETURN FROM EXIT OF DO-LOOP 200
  100 KTAR = 0
      STEP = (PSTOP-PARAM)/5.E0
C                       SET TARGET ORDERS TO NULL
      DO 110 I = 1,4
         LTRY(I) = 0
  110 CONTINUE
C                       STEP ALONG THE BOUNDARY
      DO 160 I = 1,5
         STEPI = FLOAT(I)*STEP
         CALL BCOORD(PARAM+STEPI,TESTX,TESTY,IPIECE)
         DIST = SQRT((TESTX-XSTART)**2+ (TESTY-YSTART)**2)
         IF (ALLIN) ALLIN = INSIDE(TESTX,TESTY,IX,JY,INEXT,JNEXT,XGRID,
     .                      NGDIMX,YGRID,NGDIMY)
         IF (DIST.LE.DGMAX) GO TO 120
C               HAVE GONE TOO FAR, REDUCE PSTOP
         PSTOP = PARAM + STEPI
C               IF THIS IS 1ST OR 2ND POINT, REPEAT WHOLE PROCESS
         IF (I.LE.2) GO TO 100
C               FIND SIGNS TO TARGETS FROM THIS BOUNDARY POINT
  120    SIGNOW(1) = XTAR1 - TESTX
         SIGNOW(2) = XTAR2 - TESTX
         SIGNOW(3) = YTAR1 - TESTY
         SIGNOW(4) = YTAR2 - TESTY
C
C               CHECK TO SEE IF ANY SIGN CHANGES HAVE OCCURRED
C               IF SO, SET ORDER AND SAVE GUESSES
         DO 150 L = 1,4
C                   SKIP THOSE TARGETS ALREADY IN ORDER
            IF (KTAR.EQ.0) GO TO 140
            DO 130 K = 1,KTAR
               IF (LTRY(K).EQ.L) GO TO 150
  130       CONTINUE
  140       IF (SIGNOW(L)*TSIGN(L).GE.0.0) GO TO 150
            KTAR = KTAR + 1
            LTRY(KTAR) = L
            GESS1(KTAR) = STEPI - STEP
            GESS2(KTAR) = STEPI
  150    CONTINUE
  160 CONTINUE
C
C     CHECK IF THE CURRENT BOUNDARY PIECE ENDS INSIDE THE GRID SQUARE
C           AND NO POINT EXAMINED HAS BEEN OUTSIDE THE GRID SQUARE.
      IF (INSIDE(XEND,YEND,IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,
     .    NGDIMY) .AND. ALLIN) GO TO 500
C
C         SET GUESSES AND ORDER FOR TARGETS NOT CROSSED
      DO 190 L = 1,4
C                  SKIP THOSE TARGETS ALREADY IN ORDER
         IF (KTAR.EQ.0) GO TO 180
         DO 170 K = 1,KTAR
            IF (LTRY(K).EQ.L) GO TO 190
  170    CONTINUE
  180    KTAR = KTAR + 1
         LTRY(KTAR) = L
         GESS1(KTAR) = 0.
         GESS2(KTAR) = STEP
  190 CONTINUE
      IF (LEVEL.GE.5) WRITE (MOUTPT,9051) (I,LTRY(I),GESS1(I),GESS2(I),
     .    I=1,4)
C      -----------  CASE STATEMENT TO SEARCH FOR 4 TARGETS  ------------
      DO 460 I = 1,4
         KTAR = LTRY(I)
         GO TO (200,200,330,330),KTAR
C
C               X TARGETS - CHECK FOR VERTICAL BOUNDARY BEFORE SEARCH
  200    CONTINUE
         XCASE = .TRUE.
C               TEST FOR NEARLY VERTICAL BOUNDARY
C
C        LOOP TO INCREASE DELP AND SEE IF BOUNDARY STAYS
C        ALONG THE GRID LINE UNTIL DXGRD STEP PASSED
         DELP = (PSTOP-PARAM)/8.E0
         DO 210 NVERT = 1,8
            CALL BCOORD(PARAM+DELP*FLOAT(NVERT),TESTX,TESTY,IPIECE)
            DX = ABS(TESTX-XSTART)
            IF (DX.GT.10.E0*DXGRD*EPSGRD) GO TO 230
            IF (DX.GE.DXGRD) GO TO 220
  210    CONTINUE
C
C                   BOUNDARY IS VERTICAL, SKIP THE X SEARCHES
  220    IF (LEVEL.GE.3) WRITE (MOUTPT,9061) XSTART,YSTART
         OUTCM(1) = VERT
         OUTCM(2) = VERT
         GO TO 460
C                   SELECT PROPER X TARGET
  230    IF (KTAR.EQ.2) GO TO 270
         IF (INEXT.LE.NGRIDX) GO TO 240
C                   WALKING OUTSIDE GRID DOMAIN
         IF (LEVEL.GE.1) WRITE (MOUTPT,9071) XSTART,YSTART
         OUTCM(1) = OUTSID
         OUTCM(2) = OUTSID
         GO TO 460

  240    CONTINUE
C
CASE 1 -------------------- TRY THE FIRST X-TARGET
C               SKIP THIS TARGET IF = XSTART, LOOK FOR DOUBLE BACK
         IF (ABS(XSTART-XTAR1).LE.EPSGRD) GO TO 320
C              CHECK FOR TRYING TO WALK OUTSIDE IN X-DIRECTION
         IF (OUTCM(1).EQ.OUTSID) GO TO 460
C             USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS
         CALL REGULA(XTAR1,XCASE,GESS1(I),GESS2(I),PSTOP,EPSTAN,BCOORD)
C                NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY
         CALL SECANT(PARAM+GESS1(I),PARAM+GESS2(I),XTAR1,XCASE,P(1),
     .               X(1),Y(1),OUTCM(1),PSTOP,DXGRD,BCOORD)
         IF (OUTCM(1).EQ.FAIL) GO TO 460
C              CHECK FOR HAVING WALKED BACKWARD
         IF (ABS(XSTART-X(1)).LE.EPSGRD) OUTCM(1) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF (INSIDE(X(1),Y(1),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,
     .       NGDIMY)) GO TO 250
C             POINT FOUND IS OUTSIDE GRID SQUARE, CHECK
C             FOR HAVING CROSSED THE TARGET TWICE.
         OUTCM(1) = OUTSID
         CALL CROSS2(XSTART,XTAR1,XCASE,P(1),X(1),Y(1),OUTCM(1),BCOORD)
C              CHECK TO SEE IF THE POINT IS STILL OUTSIDE
         IF ( .NOT. INSIDE(X(1),Y(1),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,
     .       YGRID,NGDIMY)) OUTCM(1) = OUTSID
  250    IF (OUTCM(1).NE.SUCCES) GO TO 460
C
C             FOUND CROSSING AT THE NEXT X-GRID LINE
  260    PNEW = P(1)
         XNEW = XTAR1
         YNEW = Y(1)
         JNEW = JY
         INEW = INEXT
C            FINISH UP AS IN THE SECOND GRID TARGET
         GO TO 300
C
CASE 2 ---------------- TRY THE SECOND X-TARGET
  270    CONTINUE
C               SKIP THIS TARGET IF = XSTART, LOOK FOR DOUBLE BACK
         IF (ABS(XSTART-XTAR2).LE.EPSGRD) GO TO 320
C             USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS
         CALL REGULA(XTAR2,XCASE,GESS1(I),GESS2(I),PSTOP,EPSTAN,BCOORD)
C                NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY
         CALL SECANT(PARAM+GESS1(I),PARAM+GESS2(I),XTAR2,XCASE,P(2),
     .               X(2),Y(2),OUTCM(2),PSTOP,DXGRD,BCOORD)
         IF (OUTCM(2).EQ.FAIL) GO TO 460
C              CHECK FOR HAVING WALKED BACKWARD
         IF (ABS(XSTART-X(2)).LE.EPSGRD) OUTCM(2) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF (INSIDE(X(2),Y(2),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,
     .       NGDIMY)) GO TO 280
C             POINT FOUND IS OUTSIDE GRID SQUARE, CHECK
C             FOR HAVING CROSSED THE TARGET TWICE.
         OUTCM(2) = OUTSID
         CALL CROSS2(XSTART,XTAR2,XCASE,P(2),X(2),Y(2),OUTCM(2),BCOORD)
C              CHECK TO SEE IF THE POINT IS STILL OUTSIDE
         IF ( .NOT. INSIDE(X(2),Y(2),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,
     .       YGRID,NGDIMY)) OUTCM(2) = OUTSID
  280    IF (OUTCM(2).NE.SUCCES) GO TO 460
C
C             FOUND CROSSING POINT AT THIS X-GRID LINE
  290    CONTINUE
         JNEW = JY
         INEW = IX
         PNEW = P(2)
         XNEW = XTAR2
         YNEW = Y(2)
  300    TYPE = VERT
C             CHECK POSSIBILITY OF TYPE BOTH
         IF (ABS(YNEW-YTAR1).GT.EPSGRD) GO TO 310
         TYPE = BOTH
         JNEW = JNEXT
         YNEW = YTAR1
         GO TO 510

  310    IF (ABS(YNEW-YTAR2).GT.EPSGRD) GO TO 510
         TYPE = BOTH
         JNEW = JY
         YNEW = YTAR2
         GO TO 510
C
C              SPECIAL CALCULATION FOR DOUBLING BACK ON TARGET
  320    CONTINUE
         CALL DBACK(XSTART,XTEST,PTEST,XCASE,DXGRD,OUTCM(KTAR),P(KTAR),
     .              X(KTAR),Y(KTAR),BCOORD,PSTOP)
         IF (OUTCM(KTAR).EQ.FAIL) GO TO 460
C              CHECK FOR HAVING WALKED BACKWARD TO STARTING POINT
         IF (ABS(YSTART-Y(KTAR)).LE.EPSGRD) OUTCM(KTAR) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF ( .NOT. INSIDE(X(KTAR),Y(KTAR),IX,JY,INEXT,JNEXT,XGRID,
     .       NGDIMX,YGRID,NGDIMY)) OUTCM(KTAR) = OUTSID
         IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.1) GO TO 260
         IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.2) GO TO 290
         GO TO 460
C
C               Y TARGETS - CHECK FOR HORIZONTAL BOUNDARY FIRST
  330    CONTINUE
         XCASE = .FALSE.
C
C        ---------- CHECK FOR NEARLY HORIZONTAL BOUNDARY
C        LOOP TO INCREASE DELP AND SEE IF BOUNDARY STAYS
C        ALONG THE GRID LINE UNTIL DYGRD IS PASSED
         DELP = (PSTOP-PARAM)/8.
         DO 340 NHORZ = 1,8
            CALL BCOORD(PARAM+FLOAT(NHORZ)*DELP,TESTX,TESTY,IPIECE)
            DY = ABS(TESTY-YSTART)
            IF (DY.GT.10.E0*DYGRD*EPSGRD) GO TO 360
            IF (DY.GT.DYGRD) GO TO 350
  340    CONTINUE
C
C                   BOUNDARY HORIZONTAL
  350    IF (LEVEL.GE.3) WRITE (MOUTPT,9081) XSTART,YSTART
         OUTCM(3) = HORZ
         OUTCM(4) = HORZ
         GO TO 460
C               SELECT THE PROPER Y TARGET
  360    IF (KTAR.EQ.4) GO TO 400
         IF (JNEXT.LE.NGRIDY) GO TO 370
C                   WALKING OUTSIDE GRID DOMAIN
         OUTCM(3) = OUTSID
         OUTCM(4) = OUTSID
         IF (LEVEL.GE.1) WRITE (MOUTPT,9071) XSTART,YSTART
         GO TO 460

  370    CONTINUE
C
CASE 3 ------------------ TRY THE FIRST Y-TARGET
C                  SKIP THIS TARGET IF = YSTART, LOOK FOR DOUBLE BACK
         IF (ABS(YSTART-YTAR1).LE.EPSGRD) GO TO 450
C              CHECK FOR TRYING TO WALK OUTSIDE IN Y-DIRECTION
         IF (OUTCM(3).EQ.OUTSID) GO TO 460
C             USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS
         CALL REGULA(YTAR1,XCASE,GESS1(I),GESS2(I),PSTOP,EPSTAN,BCOORD)
C                NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY
         CALL SECANT(PARAM+GESS1(I),PARAM+GESS2(I),YTAR1,XCASE,P(3),
     .               X(3),Y(3),OUTCM(3),PSTOP,DYGRD,BCOORD)
         IF (OUTCM(3).EQ.FAIL) GO TO 460
C              CHECK FOR HAVING WALKED BACKWARD
         IF (ABS(YSTART-Y(3)).LE.EPSGRD) OUTCM(3) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF (INSIDE(X(3),Y(3),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,
     .       NGDIMY)) GO TO 380
C             POINT FOUND IS OUTSIDE GRID SQUARE, CHECK
C             FOR HAVING CROSSED THE TARGET TWICE.
         OUTCM(3) = OUTSID
         CALL CROSS2(YSTART,YTAR1,XCASE,P(3),X(3),Y(3),OUTCM(3),BCOORD)
C              CHECK TO SEE IF THE POINT IS STILL OUTSIDE
         IF ( .NOT. INSIDE(X(3),Y(3),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,
     .       YGRID,NGDIMY)) OUTCM(3) = OUTSID
  380    IF (OUTCM(3).NE.SUCCES) GO TO 460
C
C             FOUND CROSSING POINT AT THE NEXT Y-GRID LINE
  390    CONTINUE
         JNEW = JNEXT
         INEW = IX
         PNEW = P(3)
         XNEW = X(3)
         YNEW = YTAR1
C               FINISH UP AS AT THE FOURTH TARGET
         GO TO 430
C
CASE 4 ------------ TRY THE SECOND Y-TARGET
  400    CONTINUE
C                  SKIP THIS TARGET IF = YSTART, LOOK FOR DOUBLE BACK
         IF (ABS(YSTART-YTAR2).LE.EPSGRD) GO TO 450
C         USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS
         CALL REGULA(YTAR2,XCASE,GESS1(I),GESS2(I),PSTOP,EPSTAN,BCOORD)
C                NOW USE THE SECANT METHOD TO HIT TARGET ACCURATELY
         CALL SECANT(PARAM+GESS1(I),PARAM+GESS2(I),YTAR2,XCASE,P(4),
     .               X(4),Y(4),OUTCM(4),PSTOP,DYGRD,BCOORD)
         IF (OUTCM(4).EQ.FAIL) GO TO 460
C                  CHECK FOR HAVING WALKED BACKWARD
         IF (ABS(YSTART-Y(4)).LE.EPSGRD) OUTCM(4) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF (INSIDE(X(4),Y(4),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,YGRID,
     .       NGDIMY)) GO TO 410
C             POINT FOUND IS OUTSIDE GRID SQUARE, CHECK
C             FOR HAVING CROSSED THE TARGET TWICE.
         OUTCM(4) = OUTSID
         CALL CROSS2(YSTART,YTAR2,XCASE,P(4),X(4),Y(4),OUTCM(4),BCOORD)
C              CHECK TO SEE IF THE POINT IS STILL OUTSIDE
         IF ( .NOT. INSIDE(X(4),Y(4),IX,JY,INEXT,JNEXT,XGRID,NGDIMX,
     .       YGRID,NGDIMY)) OUTCM(4) = OUTSID
  410    IF (OUTCM(4).NE.SUCCES) GO TO 460
C
C             FOUND CROSSING POINT AT THIS Y-GRID LINE
  420    CONTINUE
         PNEW = P(4)
         XNEW = X(4)
         YNEW = YTAR2
         INEW = IX
         JNEW = JY
  430    TYPE = HORZ
C             CHECK POSSIBILITY OF TYPE BOTH
C             WE GET HERE ONLY IN UNUSUAL CASES
         IF (ABS(XNEW-XTAR1).GT.EPSGRD) GO TO 440
         TYPE = BOTH
         INEW = INEXT
         XNEW = XTAR1
         GO TO 510

  440    IF (ABS(XNEW-XTAR2).GT.EPSGRD) GO TO 510
         TYPE = BOTH
         INEW = IX
         XNEW = XTAR2
         GO TO 510
C
C              SPECIAL CALCULATION FOR DOUBLING BACK ON TARGET
  450    CONTINUE
         CALL DBACK(YSTART,YTEST,PTEST,XCASE,DYGRD,OUTCM(KTAR),P(KTAR),
     .              X(KTAR),Y(KTAR),BCOORD,PSTOP)
         IF (OUTCM(KTAR).EQ.FAIL) GO TO 460
C              CHECK FOR HAVING WALKED BACKWARD TO STARTING POINT
         IF (ABS(XSTART-X(KTAR)).LE.EPSGRD) OUTCM(KTAR) = WBACK
C                  TEST FOR GETTING A RESULT IN THE RIGHT GRID SQUARE
         IF ( .NOT. INSIDE(X(KTAR),Y(KTAR),IX,JY,INEXT,JNEXT,XGRID,
     .       NGDIMX,YGRID,NGDIMY)) OUTCM(KTAR) = OUTSID
         IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.3) GO TO 390
         IF (OUTCM(KTAR).EQ.SUCCES .AND. KTAR.EQ.4) GO TO 420
C
  460 CONTINUE
C
C         HAVE CHECKED ALL POSSIBLE GRID INTERSECTIONS - NOTHING FOUND
C              CHECK THE OUTSIDE POINTS FOR ACCEPTABILITY
C              OR FOR REDUCING PSTOP AND MAKING THE WHOLE SEARCH OVER
C              AGAIN - GOING BACK UP TO 12
      POUT = PSTOP
      DO 470 I1 = 1,4
         IF (OUTCM(I1).NE.OUTSID) GO TO 470
         POUT = AMIN1(P(I1),POUT)
         IF ( .NOT. INSIDE(X(I1),Y(I1),IWEST,ISOUTH,IEAST,INORTH,XGRID,
     .       NGDIMX,YGRID,NGDIMY)) GO TO 470
C               HAVE ACCEPTABLE POINT, TAKE IT
         IF (LEVEL.GE.3) WRITE (MOUTPT,9091) X(I1),Y(I1)
C          LOCATE NEW POINT,SET COORDINATES AND EXIT BWALK
         CALL LOCATE(X(I1),Y(I1),INEW,JNEW,TYPE,XGRID,YGRID,NGDIMX,
     .               NGDIMY)
C          CHECK LOCATION FAILURE - SHOULD NEVER HAPPEN
         IF (FATAL) GO TO 470
         XNEW = X(I1)
         YNEW = Y(I1)
         PNEW = P(I1)
         GO TO 510

  470 CONTINUE
C
C               CHECK TO SEE IF PSTOP HAS BEEN REDUCED.  IF SO REDO THE
C               WHOLE SEARCH - GOING BACK UP TO 12 TO RESTART
      IF (POUT.GE.PSTOP) GO TO 490
      PSTOP = POUT
      IF (LEVEL.GE.3) WRITE (MOUTPT,9101) PSTOP,OUTCM(1),XTAR1,X(1),
     .    Y(1),P(1),OUTCM(2),XTAR2,X(2),Y(2),P(2),OUTCM(3),YTAR1,X(3),
     .    Y(3),P(3),OUTCM(4),YTAR2,X(4),Y(4),P(4)
      DO 480 I = 1,4
         X(I) = -100.E0
         Y(I) = -100.E0
         P(I) = PEND + 1.E0
         OUTCM(I) = NOTRY
  480 CONTINUE
      ALLIN = .TRUE.
      GO TO 100
C                   STILL NO POINT FOUND
C         *************************************
C         ***********  FATAL ERROR  ***********
C         *************************************
  490 CONTINUE
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9111) IPIECE,XSTART,YSTART,IX,JY,
     .    OUTCM(1),XTAR1,X(1),Y(1),P(1),OUTCM(2),XTAR2,X(2),Y(2),P(2),
     .    OUTCM(3),YTAR1,X(3),Y(3),P(3),OUTCM(4),YTAR2,X(4),Y(4),P(4)
C
      FATAL = .TRUE.
      GO TO 520
C
C                  CURRENT PIECE STOPS IN THE GRID SQUARE
  500 CONTINUE
      TYPE = INTER
      IF (LEVEL.GE.2) WRITE (MOUTPT,9121) IPIECE,XEND,YEND,IX,JY
      XNEW = XEND
      YNEW = YEND
      PNEW = PEND + .1E0*EPSGRD
C
C                     NORMAL SUBPROGRAM EXIT
  510 CONTINUE
C                   PRINT POSSIBLE BEFORE RETURN
      IF (LEVEL.GE.3) WRITE (MOUTPT,9131) PNEW,XNEW,YNEW,TYPE,INEW,
     .    JNEW, (OUTCM(I),X(I),Y(I),I=1,4)
      RETURN
C               EXIT FOR FATAL ERROR DETECTED AT 200 ABOVE
  520 FATAL = .TRUE.
      RETURN

 9001 FORMAT (/5 (3H **),' FATAL ERROR, DID NOT FIND BOUNDARY'/9X,
     .                 'POINT CLOSE TO',2F15.9/9X,
     .          'INEXPLICABLE, IT OCCURS IN FINDING INITIAL GRID SQUARE'
     .                 )
 9011 FORMAT (/5 (3H **),' FATAL ERROR, BOUNDARY GOES'/9X,
     .                 'OUTSIDE GRID DOMAIN AT',2F12.6)
 9021 FORMAT (9X,'TO POINT ',2F12.6,' WITH PARAMETER ',F12.8)
 9031 FORMAT (/5 (3H **),' UNABLE TO START SEARCH FOR GRID',
     .                 ' INTERSECTION AT ',2F12.6)
 9041 FORMAT (5X,'BOUNDARY WALK START WITH XSTART,YSTART,IX,JY,PTEST',
     .  ',XTEST,YTEST = ',2F11.6,2I4,3F11.6/9X,'PARAM = ',F10.6,
     .  ', X AND Y TARGETS ARE ',4F10.6)
 9051 FORMAT (9X,'TARGET ORDERS,GUESSES =',4 (I3,I2,2F9.5))
 9061 FORMAT (9X,'BOUNDARY IS VERTICAL AT ',2F10.6)
 9071 FORMAT (5 (' **'),'POSSIBLE FATAL ERROR, BOUNDARY GOES',
     .                ' OUTSIDE GRID DOMAIN AT ',2F11.6)
 9081 FORMAT (9X,'BOUNDARY IS HORIZONTAL AT ',2F10.6)
 9091 FORMAT (1X,8 ('+'),'THE POINT ',2F10.6,
     .                 ' TAKEN EVEN THOUGH OUTSIDE',' EXPECTED GRID')
 9101 FORMAT (9X,'ALL POINTS FOUND ARE TOO FAR AWAY. RESTART',
     .  ' SEARCH WITH SMALLER PARAMETER LIMIT',F11.6/3X,
     .  'THE RESULTS OF 4 TRIES TO FOLLOW BOUNDARY'/
     .  2 (2X,A4,'  X=',F10.6,2F14.8,F12.6)/
     .  2 (2X,A4,'  Y=',F10.6,2F14.8,F12.6))
 9111 FORMAT (/5 (3H **),' FATAL ERROR IN BOUNDARY '/3X,
     .                 'OR GRID DEFINITION FOR PIECE ',I3,' AT POINT ',
     .                 2F10.6/2X,' IN GRID ',2I3,3X,
     .                 'THE RESULTS OF 4 TRIES TO FOLLOW',
     .                 ' BOUNDARY'/2X,
     .     'NOTR - TARGET NOT TRIED         OUTS - TARGET OUTSIDE GRID '
     .                 /2X,
     .     'FAIL - TARGET NOT ACHIEVED      BACK - FOUND PREVIOUS PT.  '
     .                 /2X,'OUTCOME',4X,'TARGET',12X,'X,Y-COORDS',10X,
     .                 'PARAMETER'/2 (2X,A4,'  X=',F10.6,2F14.8,F12.6/),
     .                 2 (2X,A4,'  Y=',F10.6,2F14.8,F12.6/))
 9121 FORMAT (2X,3 ('--'),' PIECE',I3,' ENDS AT POINT',2F10.6,
     .                  '  NEAR GRID ',2I4)
 9131 FORMAT (5X,'--- END OF BOUNDARY WALK, PNEW =',F12.7/5X,
     .  'GIVES POINT ',2F10.6,' OF TYPE ',A4,' IN GRID ',2I3,
     .  ' OUTCOMES,PTS FOR 4 CASES ='/3X,4 (2X,A4,2F12.6))
      END
      SUBROUTINE CHANGE(BRANGE,NPDIM,XGRID,YGRID,NGDIMX,NGDIMY,NBDIM,
     .                  XBOUND,YBOUND,PIECE,BPTYPE,BGRID,BPARAM,BCOORD)
C
C   ***  THIS PROGRAM CHANGES BOUNDARY PIECES AND CHECKS FOR CONTINUTIY
C        IT ALSO PUTS THE END OF THE BOUNDARY PLUS THE NEXT POINT IN THE
C        BOUNDARY LIST.
C
C     INPUT = BRANGE,  (VIA COMMON) IPIECE,PARAM
C             PLUS DIMENSION VARIABLES = NPDIM,NGDIMX,NGDIMY,NBDIM
C     OUTPUT= (VIA COMMON) TYPE,PARAM,IPIECE,NBNDPT
C              XBOUND(NBNDPT)     PIECE(NBNDPT)      BPTYPE(NBNDPT)
C              YBOUND(NBNDPT)     BGRID(NBNDPT)      BPARAM(NBNDPT)
C
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
      EXTERNAL BCOORD

      CHARACTER *4 TYPEC
      REAL XGRID(NGDIMX),YGRID(NGDIMY),BRANGE(2,NPDIM)
      DIMENSION XBOUND(NBDIM),YBOUND(NBDIM),PIECE(NBDIM),BPTYPE(NBDIM),
     .          BGRID(NBDIM),BPARAM(NBDIM)
C                   COORDINATES OF END OF CURRENT PIECE
      CALL BCOORD(BRANGE(2,IPIECE),XOLD,YOLD,IPIECE)
C
C            INCREMENT COUNTER IPIECE OF BOUNDARY PIECES
C            FIND COORDINATES (XCOR,YCOR) THAT START THE NEW PIECE
      IPIECE = IPIECE + 1
      PARAM = BRANGE(1,IPIECE)
      CALL BCOORD(PARAM,XCOR,YCOR,IPIECE)
C
C             CHECK THAT THE NEXT PIECE JOINS UP WITH CURRENT PIECE
      DIST = SQRT((XOLD-XCOR)**2+ (YOLD-YCOR)**2)
      IF (DIST.LE.4.E0*EPSGRD) GO TO 10
C
      IF (LEVEL.GE.0) WRITE (MOUTPT,9001) IPIECE,XCOR,YCOR,DIST
      FATAL = .TRUE.
      RETURN
C
   10 CONTINUE
C
C            LOCATE AND TYPE CORNER POINT
      CALL LOCATE(XCOR,YCOR,IXCOR,JYCOR,TYPEC,XGRID,YGRID,NGDIMX,NGDIMY)
C
C                        DEBUG
      IF (LEVEL.GE.3) WRITE (MOUTPT,9011) IPIECE,NBNDPT,PARAM,XCOR,
     .    YCOR,TYPEC,IXCOR,JYCOR
C
C         CHECK TO SEE IF THE END POINT IS ALREADY IN BOUNDARY LIST
C         IF SO, GO ON TO THE NEXT POINT
      IF (ABS(BRANGE(2,IPIECE-1)-BPARAM(NBNDPT)).LE.EPSGRD) GO TO 20
C
C             HAVE NEW POINT (XCOR,YCOR) TO INSERT IN BOUNDARY LIST
      NBNDPT = NBNDPT + 1
      XBOUND(NBNDPT) = XCOR
      YBOUND(NBNDPT) = YCOR
      BPTYPE(NBNDPT) = TYPEC
      BGRID(NBNDPT) = IXCOR + IPACKB*JYCOR
      BPARAM(NBNDPT) = BRANGE(2,IPIECE-1)
      PIECE(NBNDPT) = IPIECE - 1
C           CHECK FOR TANGENCY AT CORNER POINT AS CHANGE IS MADE
      CALL CHKTAN(XBOUND,YBOUND,BPARAM,BPTYPE,BGRID,NBDIM,NPDIM,BRANGE,
     .            XGRID,YGRID,NGDIMX,NGDIMY,TYPEC,IXCOR,JYCOR)
C
C         *******   SUBPROGRAM EXIT AND DEBUG OUTPUT   *******
   20 CONTINUE
      IF (LEVEL.LE.1) RETURN
      TYPEC = BPTYPE(NBNDPT)
      WRITE (MOUTPT,9021) IPIECE,XCOR,YCOR,TYPEC
      WRITE (MOUTPT,9031) NBNDPT,XBOUND(NBNDPT),YBOUND(NBNDPT)
      IF (LEVEL.EQ.2) RETURN
      WRITE (MOUTPT,9041) IXCOR,JYCOR,TYPEC,BPARAM(NBNDPT),DIST
      RETURN

 9001 FORMAT (/5 (3H **),' FATAL ERROR, PIECE ',I3,' STARTS AT',
     .                 2F15.7/9X,'PIECES DO NOT JOIN UP, BREAK IS ',
     .                 1PE15.6)
 9011 FORMAT (1X,'---CHANGE TO PIECE',I3,' AT BOUNDARY POINT',I3,
     .  ' WITH PARAM = ',F10.6,', CORNER POINT IS ',2F10.6,', OF TYPE ',
     .  A4/4X,'IT IS IN GRID ',2I3)
 9021 FORMAT (/9X,'CHANGE TO BOUNDARY PIECE ',I3,' AT ',2F10.6,
     .  ' WITH TYPE = ',A4)
 9031 FORMAT (2X,3 ('--'),' IT IS BOUNDARY POINT ',I3,
     .                  ' WITH COORDINATES ',2F10.6)
 9041 FORMAT (5X,'CORNER POINT DATA IX,JY,TYPE,PARAM,DISTANCE',2I4,1X,
     .  A4,F10.6,1PE10.2)
      END
      SUBROUTINE CHKTAN(XBOUND,YBOUND,BPARAM,BPTYPE,BGRID,NBDIM,NPDIM,
     .                  BRANGE,XGRID,YGRID,NGDIMX,NGDIMY,TYPEPT,IXNU,
     .                  JYNU)
C
C  ***** THIS PROGRAM CHECKS FOR THE BOUNDARY BEING TANGENT TO A GRID
C        LINE.  IF TWO BOUNDARY POINTS ARE VERY CLOSE, THAT IS ABOUT
C        SQRT(EPSGRD) APART. THEN THEY ARE ASSUMED TO BE THE SAME POINT
C        AND THEY ARE MERGED INTO ONE POINT.  NORMALLY THE SECOND POINT
C        FOUND( HAS INDEX = NBNDPT ) IS TAKEN, BUT WHEN THE TANGENT
C        IS THE FIRST POINT ON A BOUNDARY PIECE WE TAKE IT INSTEAD.
C
      REAL XBOUND(NBDIM),YBOUND(NBDIM),BPARAM(NBDIM),XGRID(NGDIMX),
     .     YGRID(NGDIMY),BRANGE(2,NPDIM)
      CHARACTER *4 BPTYPE,TYPEC,TYPEPT
      DIMENSION BPTYPE(NBDIM),BGRID(NBDIM)
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
      DELX = ABS(XBOUND(NBNDPT)-XBOUND(NBNDPT-1))
      DELY = ABS(YBOUND(NBNDPT)-YBOUND(NBNDPT-1))
C        CHECK FOR TANGENCY TO HORIZONTAL GRID LINE
C              IF POINT IS ON HORIZONTAL GRID LINE
      IF ( .NOT. (TYPEPT.EQ.HORZ.OR.TYPEPT.EQ.BOTH)) GO TO 10
      IF ( .NOT. (DELX.LE.EPSTAN.AND.DELY.LE.EPSGRD)) GO TO 10
      DELGY2 = AMIN1(ABS(YBOUND(NBNDPT)-YGRID(JYNU)),
     .         ABS(YBOUND(NBNDPT)-YGRID(JYNU+1)))
      DELGY1 = AMIN1(ABS(YBOUND(NBNDPT-1)-YGRID(JYNU)),
     .         ABS(YBOUND(NBNDPT-1)-YGRID(JYNU+1)))
      GO TO 20
C
C        CHECK FOR TANGENCY TO VERTICAL GRID LINE
C              IF POINT IS ON VERTICAL GRID LINE
   10 IF ( .NOT. (TYPEPT.EQ.VERT.OR.TYPEPT.EQ.BOTH)) GO TO 50
      IF ( .NOT. (DELY.LE.EPSTAN.AND.DELX.LE.EPSGRD)) GO TO 50
      DELGX2 = AMIN1(ABS(XBOUND(NBNDPT)-XGRID(IXNU)),
     .         ABS(XBOUND(NBNDPT)-XGRID(IXNU+1)))
      DELGX1 = AMIN1(ABS(XBOUND(NBNDPT-1)-XGRID(IXNU)),
     .         ABS(XBOUND(NBNDPT-1)-XGRID(IXNU+1)))
C
C              HAVE FOUND TANGENCY, MERGE THE 2 POINTS
C     CHOOSE THE SECOND POINT UNLESS THE FIRST IS CLOSER TO THE
C     GRID LINES OR THE START OF A PIECE.
C     FURTHER CHECK TO SEE IF TWO INTERSECTIONS ARE SO CLOSE TO A
C     CORNER THAT THEY LOOK LIKE TANGENCY
   20 CONTINUE
      IF (DELGX2.LE.EPSGRD .AND. DELGY1.LE.EPSGRD) GO TO 50
      IF (DELGX1.LE.EPSGRD .AND. DELGY2.LE.EPSGRD) GO TO 50
C               DONT HAVE FALSE TANGENCY NEAR A CORNER
      DELGRD = AMIN1(DELGX2,DELGY2)
      DELOLD = AMIN1(DELGX1,DELGY1)
      IF (ABS(BRANGE(1,IPIECE)-BPARAM(I1NBPT-1)).LE.EPSGRD) GO TO 30
      IF (DELGRD.GT.DELOLD) GO TO 30
C              USE THE SECOND POINT
C        ITEMS NOT CHANGED HERE ARE THE SAME FOR BOTH POINTS
      XBOUND(NBNDPT-1) = XBOUND(NBNDPT)
      YBOUND(NBNDPT-1) = YBOUND(NBNDPT)
      BPTYPE(NBNDPT-1) = BPTYPE(NBNDPT)
      BPARAM(NBNDPT-1) = BPARAM(NBNDPT)
      BGRID(NBNDPT-1) = BGRID(NBNDPT)
      NBNDPT = NBNDPT - 1
      TYPEC = BPTYPE(NBNDPT)
      IF (LEVEL.GE.2) WRITE (MOUTPT,9001) NBNDPT,XBOUND(NBNDPT),
     .    YBOUND(NBNDPT),TYPEC
      GO TO 40
C
C              USE THE FIRST POINT - BUT INCREASE PARAM
   30 NBNDPT = NBNDPT - 1
      PARAM = PARAM + EPSGRD
      TYPEC = BPTYPE(NBNDPT)
      IF (LEVEL.GE.2) WRITE (MOUTPT,9011) NBNDPT,XBOUND(NBNDPT),
     .    YBOUND(NBNDPT),TYPEC
   40 IF (LEVEL.GE.3) WRITE (MOUTPT,9021) DELX,DELY
   50 RETURN

 9001 FORMAT (2X,4 ('+++'),' FOUND BOUNDARY TANGENT TO GRID AT POINT',
     .                   I4/9X,'= ',2F10.6,' OF TYPE ',A4,
     .                   ' REPLACES PREVIOUS POINT')
 9011 FORMAT (2X,4 ('+++'),' FOUND BOUNDARY TANGENT TO GRID AT POINT',
     .                   I4/9X,'= ',2F10.6,' OF TYPE ',A4,
     .                   ', MOVED PAST IT')
 9021 FORMAT (8X,' TEST VARIABLES DELX,DELY =',4E12.3)
      END
      SUBROUTINE CROSS2(START,TARGET,XCASE,PEND,XEND,YEND,OUTCOM,BCOORD)
C
C ****  THIS ROUTINE IS TO CHECK FOR THE POSSIBILITY THAT THE BOUNDARY
C       CROSSED THE TARGET GRID LINE MORE THAT ONCE AND THAT THE
C       FIRST CROSSING WAS NOT LOCATED.  THE IDEA IS TO CHECK TO SEE
C       IF THE BOUNDARY IS MOVING TOWARD OR AWAY FROM THE START.  IF
C       A DOUBLE( OR MORE ) CROSSING IS DETECTED THEN REGULA ( THE
C       MODIFIED REGULA FALSI METHOD) IS USED TO LOCATE A CROSSING
C       CLOSER TO THE START.  NOTE THIS LOGIC WILL NOT NOT DETECT A
C       TRIPLE CROSSING.
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
      CHARACTER *4 OUTCOM,SUCCES
      EXTERNAL BCOORD

      LOGICAL XCASE
      DATA SUCCES/'SUCC'/
C
C            FIRST CHECK FOR A DOUBLE CROSSING
      PCHECK = PEND - EPSTAN
      CALL BCOORD(PCHECK,XCHECK,YCHECK,IPIECE)
      CHECK = YCHECK
      IF (XCASE) CHECK = XCHECK
      SIGN = TARGET - START
      CSIGN = TARGET - CHECK
      IF (LEVEL.GE.3) WRITE (MOUTPT,9001) PEND,SIGN,CSIGN
      IF (SIGN*CSIGN.GT.0.) RETURN
C             FOUND A SIGN CHANGE = THERE IS CLOSER CROSSING
C             USE REGULA TO FIND IT ACCURATELY
      PLEFT = 0.
      PRIGHT = PCHECK - PARAM
      CALL REGULA(TARGET,XCASE,PLEFT,PRIGHT,PCHECK-PARAM,EPSGRD,BCOORD)
C          REGULA SHOULD ALWAYS CONVERGE SINCE THE TARGET IS BRACKETED
C          CHECK CONVERGENCE
      CALL BCOORD(PARAM+PRIGHT,X,Y,IPIECE)
      V = Y
      IF (XCASE) V = X
      IF (ABS(TARGET-V).GT.EPSGRD) RETURN
C
C          HAVE CONVERGENCE, SET VALUES AND OUTCOM
      OUTCOM = SUCCES
      IF (LEVEL.GE.3) WRITE (MOUTPT,9011) XEND,YEND,X,Y
      PEND = PRIGHT
      XEND = X
      YEND = Y
      RETURN

 9001 FORMAT (9X,'***** ','DOUBLE CROSSING CHECK AT',F12.6,' HAS SIGNS',
     .  2E13.3)
 9011 FORMAT (9X,'***** ','DOUBLE CROSSING CHECK FOUND CLOSER POINT'/
     .  15X,'ORIG OUTSIDE PT =',2F12.6,' NEW PT =',2F12.6)
      END
      SUBROUTINE DBACK(TARGET,VTEST,DELPST,XCASE,DVGRD,OUTCOM,P,X,Y,
     .                 BCOORD,PSTOP)
C
C  ***** THIS PROGRAM ATTEMPTS TO FIND A CROSSING DUE TO THE
C        BOUNDARY DOUBLING BACK OVER A GRID LINE WITHIN THE GRID
C        SQUARE. THE BOUNDARY IS FOLLOWED FROM (XTEST,YTEST) UNTIL
C        IT CROSSES THE TARGET AGAIN.  THIS GIVES PARAMETERS P1,P2
C        WHICH ARE USED TO START SECANT SEARCH FOR THE RESULTS
C        P,X,Y = SOLUTION.
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
      EXTERNAL BCOORD

      CHARACTER *4 OUTCOM,FAIL
      LOGICAL XCASE
      DATA FAIL/'FAIL'/
C               ALWAYS LIMIT PARAMETER TO CURRENT PIECE
      DELPMX = PSTOP - PARAM
      FACTOR = 1.5E0
C         DEVST = SIZE OF POINT-TARGET AT START
C         DEVNU = SIZE OF POINT-TARGET AT NEXT STEP
      DEVST = VTEST - TARGET
      DELP = DELPST
      KEXPND = 0
      IF (LEVEL.GE.5) WRITE (MOUTPT,9001) KEXPND,DELP,DEVST,VTEST
C
C               EXPAND PARAMETER STEP UNTIL BOUNDARY
C               CROSSES THE TARGET GRID LINE
      DO 10 KEXPND = 1,20
C               STEP AWAY FROM PREVIOUS POINT TO CROSS THE TARGET
         DELPNU = AMIN1(FACTOR*DELP,DELPMX)
         CALL BCOORD(PARAM+DELPNU,XNU,YNU,IPIECE)
         VNU = YNU
         IF (XCASE) VNU = XNU
         DEVNU = VNU - TARGET
C             TEST FOR FINDING A SIGN CHANGE IN POINT-TARGET DIFFERENCE
C             THIS MEANS THAT WE HAVE CROSSED THE TARGET GRID LINE
C             SO WE START TO USE SECANT METHOD
         IF (DEVST*DEVNU.LT.0.E0) GO TO 30
         DELP = DELPNU
C                       INCREASE FACTOR OCCAISIONALLY
         IF (MOD(KEXPND,4).EQ.0) FACTOR = FACTOR*1.2E0
C                   CHECK FOR GETTING TOO FAR AWAY = FAILURE
         IF (ABS(DEVNU).GE.DVGRD .OR. DELPNU.GE.DELPMX) GO TO 20
         IF (LEVEL.GE.5) WRITE (MOUTPT,9001) KEXPND,DELPNU,DEVNU,XNU,
     .       YNU
   10 CONTINUE
C
C             FAILED TO FIND A DOUBLE CROSSING = FAILURE EXIT
   20 OUTCOM = FAIL
      IF (LEVEL.GE.4) WRITE (MOUTPT,9011) TARGET,VTEST,XNU,YNU
      RETURN
C
C             FOUND CROSSING, SET PARAMETER GUESSES P1 AND P2
C             CALL SECANT TO FINISH FINDING CROSSING
   30 CONTINUE
      P1 = PARAM + DELPNU
      P2 = PARAM + DELP
      IF (LEVEL.GE.5) WRITE (MOUTPT,9021) P1,P2
C             USE SPECIAL REGULA FALSI TO IMPROVE THE INITIAL GUESS
      CALL REGULA(TARGET,XCASE,DELPNU,DELP,PSTOP,EPSTAN,BCOORD)
C             USE SECANT TO GET THE FINAL RESULTS
      CALL SECANT(PARAM+DELPNU,PARAM+DELP,TARGET,XCASE,P,X,Y,OUTCOM,
     .            PSTOP,DVGRD,BCOORD)
      IF (OUTCOM.EQ.FAIL .AND. LEVEL.GE.4) WRITE (MOUTPT,9011) TARGET,
     .    VTEST,XNU,YNU
      RETURN

 9001 FORMAT (9X,'DBACK LOOP',I3,' GIVES DELPNU,DEVNU,X,Y =',4F11.6)
 9011 FORMAT (9X,'DBACK ROUTINE FAILS FOR TARGET ',F10.6,', START AT ',
     .  F11.7,', END AT ',2F11.7)
 9021 FORMAT (9X,'DBACK SUCCEEDED TO BRACKET TARGET',
     .  ' WITH PARAMETER = ',2F10.6)
      END
      SUBROUTINE EXPAND(ICENT,JCENT,NUPTS,INEW,JNEW,GTYPE,NGDIMX,NGDIMY)
C
C  *****  THIS PROGRAM IS GIVEN THE POINT ICENT,JCENT AS AN INTERIOR
C     POINT.  IT THEN DETERMINES IF ITS 4 NEIGHBORS ARE INTERIOR
C     POINTS AND, IF SO, MARKS THEM BY
C             GTYPE =  NSIDE FOR PTS AWAY FROM THE BOUNDARY
C             GTYPE =  -GTYPE FOR PTS NEXT TO THE BOUNDARY
C
C
      INTEGER N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,LEVEL,
     .        QLIMIT,QX(250),QY(250),GTYPE,PIECE,BGRID,BNEIGH
      REAL PARAM,EPSGRD,EPSTAN
      LOGICAL CLOCKW,ARC,FATAL,INHOLE
      CHARACTER *4 TYPE,BPTYPE,HORZ,VERT,BOTH,INTER,JUMP
      COMMON /DMCINT/N1BND,N1BDPT,NSIDE,IPACKB,MOUTPT,NGRIDX,NGRIDY,
     .       LEVEL,QLIMIT,QX,QY
      COMMON /DMCREL/PARAM,EPSGRD,EPSTAN
      COMMON /BNDRYI/IPIECE,NBOUND,NBNDPT
      COMMON /BNDRYL/CLOCKW,ARC,FATAL,INHOLE
      COMMON /DMCHAR/TYPE,HORZ,VERT,BOTH,INTER,JUMP
C
      DIMENSION GTYPE(NGDIMX,NGDIMY)
C
      INTEGER INEW(4),JNEW(4),PTYPE,IXN(4),JYN(4),CENTER
      LOGICAL PTNAY,BCROSS
      DATA IXN(1),IXN(2),IXN(3),IXN(4)/0,-1,0,1/
      DATA JYN(1),JYN(2),JYN(3),JYN(4)/-1,0,1,0/
C
C     ARITHMETIC STATEMENT FUNCTIONS
C             PTNAY IS TRUE FOR BOUNDARY POINT LOCATED
C                  ABOVE   J=1
C                  RIGHT   J=2
C                  BELOW   J=3
C                  LEFT    J=4
      PTNAY(K,J) = MOD(K/ (IPACKB*2** (MOD(J-1,4))),2) .EQ. 1
C
C         BCROSS IS TRUE IF THE LINE FROM THE J-TH POINT TO THE
C         CENTER POINT CROSSES THE BOUNDARY
      BCROSS(J) = PTNAY(PTYPE,J) .AND. PTNAY(CENTER,J+2)
C
      CENTER = GTYPE(ICENT,JCENT)
      NUPTS = 0
      DO 10 JNAY = 1,4
         IX = ICENT + IXN(JNAY)
         JY = JCENT + JYN(JNAY)
         PTYPE = -GTYPE(IX,JY)
C             PTYPE IS NEGATIVE FOR BOUNDARY POINTS
         IF (PTYPE.LT.0) GO TO 10
         IF (BCROSS(JNAY)) GO TO 10
C                  HAVE AN INTERIOR POINT
C                  CHECK FOR LEGALITY OF INDEXES
         IF (IX.LT.1 .OR. IX.GT.NGRIDX) GO TO 10
         IF (JY.LT.1 .OR. JY.GT.NGRIDY) GO TO 10
         NUPTS = NUPTS + 1
         INEW(NUPTS) = IX
         JNEW(NUPTS) = JY
         IF (PTYPE.EQ.0) PTYPE = NSIDE
         GTYPE(IX,JY) = PTYPE
   10 CONTINUE
      IF (LEVEL.GE.4) WRITE (MOUTPT,9001) ICENT,JCENT,PTYPE
      IF (LEVEL.GE.5 .AND. NUPTS.GT.0) WRITE (MOUTPT,9011) (I,INEW(I),
     .    JNEW(I),I=1,NUPTS)
      RETURN

 9001 FORMAT (5X,'EXPAND INTERIOR AT ',2I3,' FINAL TYPE =',I7)
 9011 FORMAT (9X,'NEW POINTS =  ',4 (I7,' = ',2I3))
      END
      SUBROUTINE FILL(GTYPE,NGDIMX,NGDIMY,NBDIM,BGRID,BPTYPE,PIECE,
     .                XBOUND,YBOUND)
C
C  *****  THIS PROGRAM LOCATES AN INTERIOR POINT NEAR THE START OF THE
C     BOUNDARY AND THEN EXPANDS THE INTERIOR, MARKING ALL INTERIOR
C     POINTS BY CHANGING GTYPE TO PLUS OR SETTING GTYPE = INTER.
C     THE METHOD USES A QUEUE OF POINTS ON THE FRINGE OF THE CURRENTLY
C     KNOWN INTERIOR AND THIS QUEUE IS PROCESSED VIA EXPAND WHICH
C     LOCATES NEW INTERIOR POINTS ADJACENT TO ANY GIVEN ONE.
C
C     INPUT.  XBOUND,YBOUND = BOUNDARY POINT COORDINATES
C             BGRID         =    -       -   GRID LOCATIONS
C