      PROGRAM  DRIVER
C     DUMMY DRIVER FOR MESH90,  JULY 1990
      PARAMETER ( KFLT = 20000, KFIX = 15000, MINNPT = 1,
     1            MINEL = 1 )
      DIMENSION AD(KFLT), ID(KFIX)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
C     DEFINE I/O UNITS: CARDS, PRINTER, PUNCH, PLOT
      NCRD = 5
      NPRT = 6
      NPUN = 7
      OPEN (UNIT=NPUN, FILE='PUNCH', STATUS='UNKNOWN')
      NPLT = 8
      OPEN (UNIT=NPLT, FILE='CALCOMP', STATUS='UNKNOWN')
C     PROBLEM UNITS: SUMMARY, PTS, ELEMS, BCS, BCVALUE,
C             MATERIAL CODES, IN PRINT-PLOTS, OUT PRINT-PLOTS
      NPIO = 10
      OPEN (UNIT=NPIO, FILE='SUMMARY', STATUS='UNKNOWN')
      NPPT = 11
      OPEN (UNIT=NPPT, FILE='POINTS', STATUS='UNKNOWN')
      NPEL = 12
      OPEN (UNIT=NPEL, FILE='ELEM', STATUS='UNKNOWN')
      NPBC = 13
      OPEN (UNIT=NPBC, FILE='BCNULL', STATUS='UNKNOWN')
      NPBV = 14
      OPEN (UNIT=NPBV, FILE='BCVAL', STATUS='UNKNOWN')
      NPLM = 15
      OPEN (UNIT=NPLM, FILE='MAT', STATUS='UNKNOWN')
      NPIP = 16
      OPEN (UNIT=NPIP, FILE='INPLOT', STATUS='UNKNOWN')
      NPOP = 17
      OPEN (UNIT=NPOP, FILE='OUTPLOT', STATUS='UNKNOWN')
C     INITIALIZE
      DO 10 K = 1,KFLT
 10   AD(K) = 0.
      DO 20 K = 1,KFIX
 20   ID(K) = 0
      CALL  SETUP  (KFLT,KFIX,AD,ID,MINNPT,MINEL)
C     CLOSE THE UNITS HERE
      STOP
      END
      SUBROUTINE  SETUP (KFLT,KFIX,AD,ID,MINNPT,MINEL)
C      * * * * * * * * * * * * * * * * * * * * * * * * * * *
C      FEM MESH GENERATION FOR TRIANGLES AND QUADRILATERALS
C      * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     TITLE IS THE NAME OF THE PROBLEM,  ETC.
C     NPTCD = NUMBER OF NODAL POINT INPUT CARDS
C     MAXR = MAX. VALUE OF J IN DATA (1 + MULTIPLE OF NEDGE)
C     MAXC = MAX. VALUE OF I IN DATA (1 + MULTIPLE OF NEDGE)
C     FORMATS FOR PUNCHED OUTPUT (IF ANY)
C         IOPT = 0, AKIN MODEL CODE
C              = 1  AKIN FINITE CODE
C              = 2  PAFEC
C              = 3  NAVIER-STOKES CODE
C              > 0, USER SUPPLIED FORMATS
C     OUTPT AND OUTEL GOVERN THE METHOD OF OUTPUT FOR
C        THE POINTS AND ELEMENTS. THE OPTIONS ARE:
C        = 0, PRINT
C        = 1, PRINT AND PUNCH
C        = 2, PUNCH ONLY
C        = 3, NO OUTPUT
C     LTYPE = ELEMENT TYPE SELECTION
C        = 3, 3 NODE TRIANGLE     NEDGE = 2
C        = 4, 4 NODE QUAD                 2
C        = 6, 6 NODE TRIANGLE             3
C        = 8, 8 NODE QUAD                 3
C        = 9, 9 NODE QUAD                 3
C        = 10, 10 NODE TRIANGLE           4
C        = 12, 12 NODE QUAD               4
C        = 15, 15 NODE TRIANGLE           5
C        = 16, 16 NODE LAGRANGIAN QUAD    5
C        =-16, 16 NODE QUAD               6
C        = 17, 17 NODE QUAD               6
C        = 25, 25 NODE QUAD               6
C     THE UNSPECIFIED POINTS ARE DETERMINED BY:
C        MAP = 0 , MAPPING USING LAPLACE'S EQ (ITERATIVE)
C            = 1 , TWO POINT LINEAR INTERPOLATION BETWEEN
C                  BOUNDING POINTS ON THE J LINE
C            = 2 , FOUR POINT INTERPOLATION BETWEEN GIVEN
C                  POINTS ON THE I AND J LINES
C            = 3 , EIGHT POINT INTERPOLATION BETWEEN GIVEN
C                  POINTS ON THE I AND J LINES PLUS FOUR
C                  'CORNER' POINTS (COON'S PATCH)
C     NMTCD = NUMBER OF ELEMENT CODE INPUT CARDS
C     IF ECHO = 0 THEN INPUT DATA ARE LISTED
C     LPLOT = PLOTTING OPTIONS
C        = 0, NO PLOTS
C        = 1, PLOT R,Z MESH
C        = 2, PLOT I,J MESH
C        = 3, PLOT BOTH R,Z AND I,J MESHES
C     NEDCD = NO. CARDS FOR EDGE FLUX CALCULATIONS = 0
C----------------------------------------------------------
C     LOGICAL UNITS USED BY MESH90:
C        5 = CARD READER,  6 = LINE PRINTER
C        7 = CARD PUNCH, WHEN OUTPT OR OUTEL =  1 OR 2
C----------------------------------------------------------
      INTEGER  OUTPT, OUTEL, ECHO
      DIMENSION  TITLE(20), AD(KFLT), ID(KFIX)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / BUGS / MAXBUG
      DATA IJFLAG, IJ2JI  / -999, 1 /
C            *** PRINT AUTHOR CREDIT ***
C     CALL  TOOT
      IERR = 0
C-->     *** READ,PRINT, AND CHECK INPUT DATA ***
      READ (NCRD,5000) TITLE
 5000 FORMAT (20A4)
      READ (NCRD,5010) MAXR, MAXC, NPTCD, NMTCD, LTYPE, LPLOT,
     1                 MAP, OUTPT, OUTEL, IOPT, ECHO, NEDCD
 5010 FORMAT (16I5)
      IF ( LTYPE .LT. 3 .AND. LTYPE . NE. -16 )  LTYPE = 3
C      TEST FOR DEBUG OUTPUT REQUEST
      MAXBUG = 0
      IF ( MAXR .LT. 0 )  THEN
        MAXBUG = 1
        MAXR   = IABS ( MAXR )
      ENDIF
      WRITE (NPRT,5020)  TITLE
 5020 FORMAT ( /, 1X, 20A4,/ )
C      CHECK FOR POOR BAND OR SKYLINE
      IJ2J1 = 1
      IF ( MAXR .LT. MAXC )  THEN   
        IJ2JI = 0
        CALL SWAP (MAXC,MAXR)
        WRITE (NPRT,*) ' NOTE: I,J REVERSED TO REDUCE BANDWIDTH'
      ENDIF   
      WRITE (NPRT,5030) NPTCD, NMTCD, NEDCD, MAXC, MAXR, OUTPT,
     1                  MAP, OUTEL, LPLOT, IOPT, ECHO, LTYPE
 5030 FORMAT (' ****  M E S H 9 0  ****',//,
     1' NUMBER OF INPUT CARDS FOR:',/,
     2'            NODAL POINTS..............',I5,/,
     3'            ELEMENT CODES.............',I5,/,
     4'            EDGE CONDITIONS...........',I5,//,
     5' MAXIMUM INPUT VALUES:  I =',I3,', J =',I3,//,
     6' CONTROL PARAMETERS:',/,
     7'     OUTPT....',I3,10X,'MAP.......',I3,/,
     8'     OUTEL....',I3,10X,'LPLOT.....',I3,/,
     9'     IOPT.....',I3,10X,'ECHO......',I3,/,
     1'     ELTYPE...',I3,/)
C      DETERMINE I,J INCREMENTS FOR LTYPE
      LINC = 1
      IF ( LTYPE.EQ.6  .OR. LTYPE.EQ.8  .OR. LTYPE.EQ.9  ) LINC = 2
      IF ( LTYPE.EQ.10 .OR. LTYPE.EQ.12 .OR. LTYPE.EQ.16 ) LINC = 3
      IF ( LTYPE.EQ.-16 )                                  LINC = 4
      IF ( LTYPE.EQ.15 .OR. LTYPE.EQ.17 .OR. LTYPE.EQ.25 ) LINC = 4
C      SET NUMBER OF NODES PER EDGE
      NPERE = LINC + 1
C      CHECK DATA
      MAXRL1 = MAXR - 1
      IF ( (MAXRL1/LINC)*LINC .NE. MAXRL1 )  THEN
        WRITE (NPRT,*) 'ERROR: MAXR NOT VALID FOR TYPE', LTYPE
        IERR = 1
      ENDIF
      MAXCL1 = MAXC - 1
      IF ( (MAXCL1/LINC)*LINC .NE. MAXCL1 )  THEN
        WRITE (NPRT,*) 'ERROR: MAXC NOT VALID FOR TYPE', LTYPE
        IERR = 1
      ENDIF
C
C     FIXED PT ARRAY STORAGE LI: 1-INNER 2-MAT 3-IMAX
C            4-IMIN 5-LIMIN 6-LIMAX 7-IGIVEN 8=1-NPTNUM
C     FLOATING PT ARRAY STORAGE MI: 1-R 2-Z 3-CODE 4-ALINE
C
C-->         *** CALCULATE ARRAY POINTERS ***
      CALL PTER (L1,L2,L3,L4,L5,L6,L7,L8,LTOTAL,
     1           M1,M2,M3,M4,MTOTAL,MAXR,MAXC)
C               *** POINTERS COMPLETED ***
      WRITE (NPRT, 5040)  KFLT, MTOTAL, KFIX, LTOTAL
 5040 FORMAT (' ARRAY STORAGE: AVAILABLE   REQUIRED',/,
     1        ' FLOATING POINT',I10,1X,I10,/,
     2        ' FIXED POINT   ',I10,1X,I10)
      IF ( KFLT .LT. MTOTAL .OR. KFIX .LT. LTOTAL)  THEN
        WRITE (NPRT,*) 'DIMENSION STATEMENT EXCEEDED.  STOP'
        STOP     
      ENDIF   
      CALL  MESH90 (MAXR, MAXC, MINNPT, MINEL, NPTCD,
     1              NMTCD, NEDCD, ECHO, IOPT, LINC,
     2              LTYPE, MAP, OUTEL, OUTPT, LPLOT,
     3              ID(L1), ID(L2), ID(L3), ID(L4),
     4              ID(L5), ID(L6), ID(L7), ID(L8), 
     5              AD(M1), AD(M2), AD(M3), AD(M4),
     6              AD, ID, KFLO, KFIX, IJ2JI, IJFLAG)
      RETURN
      END
      SUBROUTINE  MESH90 (MAXR, MAXC, MINNPT, MINEL, NPTCD,
     1                    NMTCD, NBCCD, ECHO, IOPT, LINC,
     2                    LTYPE, MAP, OUTEL, OUTPT, LPLOT,
     3                    INNER, MAT, IMAX, IMIN,
     4                    LIMIN, LIMAX, IGIVEN, NPTNUM,
     5                    R, Z, CODE, ALINE, 
     6                    AD, ID, KFLO, KFIX, IJ2JI, IJFLAG)
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     FEM MESH GENERATION FOR TRIANGLES AND QUADRILATERALS
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
      INTEGER  OUTPT, OUTEL, ECHO
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / BUGS / MAXBUG
      DIMENSION  ALINE(MAXR), CODE(MAXC,MAXR), R(MAXC,MAXR),
     1           Z(MAXC,MAXR), IGIVEN(MAXC,MAXR), IMIN(MAXR),
     2           IMAX(MAXR), LIMAX(MAXR), LIMIN(MAXR),
     3           MAT(MAXC,MAXR), NPTNUM(MAXC,MAXR),
     4           INNER(MAXC,MAXR)
      DIMENSION  AD(KFLO), ID(KFIX)
C     MAXR  - Max. number of input rows (J max). See notes.
C     MAXC  - Max. number of input columns (I max).
C     MINNPT- Minimum node point number for this mesh
C     MINEL - Minimum element number for this mesh
C     NPTCD - Number of nodal point card generation sets.
C     NMTCD - Number of material code change sets.
C     LTYPE - Element type. That is, nodes per element: 
C              3 - linear triangle
C              4 - bilinear quadrilateral
C              6 - quadratic triangle
C              8 - serendipity biquadratic quadrilateral
C              9 - lagrangian biquadratic quadrilaterial
C             10 - full cubic triangle ( Pafec 9t )
C             12 - serendipity bicubic quadrilateral
C             15 - full quartic triangle
C             16 - lagrangian bicubic quadrilateral
C            -16 - serendipity biquartic quadrilateral
C             17 - serendipity biquartic quadrilateral w center
C             25 - lagrangian biquartic quadrilateral
C     LPLOT - Element plotting option to file poltter.
C              0 = no plots
C              1 = plot physical mesh
C              2 = plot parent mesh
C              3 = plot both meshes
C     MAP  -  Mesh mapping option.
C              0 = Iterative Laplace smoothing (default).
C              1 = Linear interpolation on J = constant.
C              2 = Four point bilinear interpolation.
C              3 = Eight point blending (Coons) interpolation.
C     OUTP  -  Nodal output options.
C              0 = list only
C              1 = list and punch for IOPT program
C              2 = punch for IOPT program only
C              3 = no output
C             This also applies to default essential bc values.
C     OUTL  -  Element output option, same as OUTP.
C     IOPT  -  Format option
C              0 = MODEL Code
C              1 = FINITE Code
C              2 = PAFEC Code
C              3 = FLOWNS Navier-Stokes
C     ECHO  -  Echo input data if 0.
C     NBCCD -  Number of essential bc value generation sets.
C     ALINE -  A line for print-plot
C     CODE  -  Boundary condition code for each I, J node
C     R     -  X coordinate array
C     Z     -  Y coordinate array
C     IGIVEN-  Flags all input nodes on the grid
C     INNER -  Array defineing interior nodes and holes
C     MAT   -  Material codes at I, J
C     IMAX  -  Max. I for each J node
C     IMIN  -  Min. I for each J node
C     LIMAX -  Max. I for each J element
C     LIMIN -  Min. I for each J element
C     NPTNUM-  Generated node point numbers array
C             *** INITIALIZE ARRAYS ***
      CALL  INITAL  (MAXR,MAXC,MAT,IMIN,LIMIN)
C
C-->          *** READ INPUT POINTS ***
C
      CALL  INPTS (MAXC,MAXR,R,Z,CODE,INNER,IMIN,IMAX,
     1             LIMIN,LIMAX,ECHO,NPTCD,MAP,IJ2JI)
      WRITE (NPRT,*) ' END OF GEOMETRY INPUT'
C             *** PRINT MAP OF INPUT POINTS ***
      CALL IJMESH (MAXC,MAXR,ALINE,INNER,IMAX,IMIN,1)
C             *** PRINT MAP OF INPUT BC FLAGS ***
      CALL IJMESH (MAXC,MAXR,ALINE,CODE,IMAX,IMIN,2)
C
C-->        *** READ CHANGES IN ELEMENT CODES ***
C
      IF ( NMTCD .GT. 0 )  THEN
        CALL INEL (MAXC,MAXR,MAT,LIMIN,LIMAX, ECHO,NMTCD,
     1             IJ2JI,IJFLAG, LINC,IMIN,IMAX)
        WRITE (NPRT,*) ' END OF MATERIAL CHANGES'
      ENDIF
C
C-->                *** MAP POINTS ***
C
      IF ( LPLOT .EQ. 2 .AND. OUTPT .EQ. 3 )  GO TO 15
      IF ( MAP .EQ. 0 )  CALL LAPLAC (MAXC,MAXR,R,Z,INNER,
     1                                IMIN,IMAX,IGIVEN,N1)
      IF ( MAP .EQ. 1 )  CALL TWOPT  (MAXC,MAXR,R,Z,INNER,
     1                                IMIN,IMAX,N1)
      IF ( MAP .EQ. 2 )  CALL FOURPT (MAXC,MAXR,R,Z,INNER,
     1                                IMIN,IMAX,IGIVEN,N1)
      IF ( MAP .EQ. 3 )  CALL EIGHTP (MAXC,MAXR,R,Z,INNER,
     1                                IMIN,IMAX,IGIVEN,N1)
      IF ( OUTPT .LT. 2 )  WRITE (NPRT,*)  
     1    'COORDINATE ITERATIONS USED =', N1
C
C-->                  *** NUMBER POINTS ***
C
   15 CONTINUE
      IF ( LTYPE .EQ. 8 .OR. LTYPE .EQ. 12 .OR.
     1     LTYPE .EQ. 17 )  CALL OMITP (MAXC,MAXR,MAT,LTYPE,
     2                                  LINC,IJFLAG,R,Z)
      CALL  PCOUNT (MAXC,MAXR,R,Z,CODE,MAT,NPTNUM,IMIN,IMAX,
     1              IJFLAG,IOPT,OUTPT,MINNPT,LINC)
C                   DEBUG LISTS
      IF ( MAXBUG .NE. 0 )  THEN
        WRITE (NPRT,*) ' Material Numbers'
        CALL IPRINT (MAT,MAXC,MAXR)
        WRITE (NPRT,*) ' Nodal Numbers'
        CALL IPRINT (NPTNUM,MAXC,MAXR)
      ENDIF
C              *** PRINT MAP OF ACTUAL NODES ***
      CALL IJNODE (MAXC,MAXR,ALINE,NPTNUM,IMAX,IMIN)
      CALL IJMESH (MAXC,MAXR,ALINE,NPTNUM,IMAX,IMIN,3)
C
C-->               *** GENERATE ELEMENTS ***
C             (AND OVERWRITE ARRAY IGIVEN WITH LEMNUM)
      IF ( OUTEL .LT. 4 )  CALL  LCOUNT (MAXC,MAXR,R,Z,LIMAX,
     1                     LIMIN,MAT,NPTNUM,LTYPE,LINC,MINEL,
     2                     OUTEL,IOPT,IGIVEN)
C          *** PRINT MAP OF ACTUAL ELEMENT NUMBERS***
      CALL IJELEM (MAXC,MAXR,ALINE,IGIVEN,LIMAX,LIMIN)
C            *** PRINT MAP OF MATERIAL CODES ***
      CALL IJMAT (MAXC,MAXR,ALINE,IGIVEN,LIMAX,LIMIN,MAT,LINC)
C
C-->            *** GENERATE BOUNDARY CONDITIONS ***
C
C             *** PUNCH DEFAULT ESSENTIAL BC FLAGS ***
      IF (  OUTPT .EQ. 1 .OR. OUTPT .EQ. 2 )
     1      CALL  PBCOUT (MAXC,MAXR,NPTNUM,CODE,IOPT)
C-->        *** PUNCH GIVEN ESSENTIAL BC VALUES ***
C                ( AND OVERWRITE INNER )
      IF (  NBCCD.GT.0 .AND. OUTPT .EQ. 1 .OR. OUTPT .EQ. 2 )
     1      CALL  VALUE  (MAXC,MAXR,ALINE,INNER,IMIN,IMAX,
     2                    NPTNUM,ECHO,NBCCD,IJ2JI)
C
C-->                 *** PLOT RESULTS ***
C
      IF ( LPLOT .EQ. 0 )  GO TO 20
C        INITIALIZE PLOT FILE AND SET FIRST ORIGIN
      CALL  STARTP
C        PLOT ACTUAL R-Z MESH
      IF ( LPLOT .EQ. 1 .OR. LPLOT .EQ. 3 )  
     1     CALL  PLOTIJ (MAXC,MAXR,R,Z,NPTNUM,MAT,IMAX,IMIN,
     2                   LIMIN,LIMAX,LINC,LTYPE,LPLOT)
      IF ( LPLOT .LT. 2 )  GO TO 20
C        ADVANCE PLOTTER FOR NEXT FRAME
      IF (  LPLOT .EQ. 3 )  CALL ADVANC
C        REPLACE R,Z BY I,J VALUES AND PLOT
      CALL  IJTORZ (MAXC,MAXR,R,Z)
      CALL  PLOTIJ (MAXC,MAXR,R,Z,NPTNUM,MAT,IMAX,
     1              IMIN,LIMIN,LIMAX,LINC,LTYPE,LPLOT)
   20 CONTINUE
C        CLOSE PLOT FILE
      IF (  LPLOT .GT. 0 )  CALL  FINSHP
      WRITE (NPR,*) ' NORMAL ENDING OF MESH/2D'
      RETURN
      END
      SUBROUTINE  ADVANC
C     - - - - - - - - - - - - - - - - - - - - - - - - - - -
C      INSTALLATION DEPENDENT PLOT FRAME ADVANCE
C     - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      EXTERNAL PLOT
C     CALCOMP CONTROL
      CALL PLOT(12.,0.,-3)
      RETURN
      END
      SUBROUTINE  EIGHTP (MC,MR,R,Z,INNER,IMIN,IMAX,IGIVEN,N1)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     EIGHT POINT COON'S NONLINEAR INTERPOLATION
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(MC,MR), Z(MC,MR), INNER(MC,MR), IMIN(MR),
     1           IMAX(MR), IGIVEN(MC)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C        DEFINE COONS BLENDING FUNCTIONS
      F0(X) = 10.*X**3 - 15.*X**4 + 6.*X**5
      F1(X) = 1.-10.*X**3+15.*X**4-6.*X**5
      JLESS = MR - 1
      N1 = 1
      IF ( JLESS .LE. 1 ) GO TO 90
      DO 80 J = 2,JLESS
      NG = 0
      I1 = IMIN(J)
      I2 = IMAX(J)
      DO 10 I = I1,I2
      IF ( INNER(I,J) .NE. 1 ) GO TO 10
      NG = NG + 1
      IGIVEN(NG) = I
   10 CONTINUE
C        SKIP TO NEXT ROW IF ALL POINTS ARE GIVEN
      IF ( NG .EQ. (I2-I1+1) ) GO TO 80
      NGLESS = NG - 1
      DO 70 IG = 1,NGLESS
      IGPLUS = IG + 1
C        FIND LIMITING I POINTS
      IF ( (IGIVEN(IGPLUS)-IGIVEN(IG)) .EQ. 1 ) GO TO 70
      IBACK = IGIVEN(IG)
      IFRONT = IGIVEN(IGPLUS)
      I1 = IBACK + 1
      I2 = IFRONT - 1
      DO 60 I = I1,I2
C        FIND LIMITING J POINTS
      J1 = J + 1
      DO 20  JJ = J1,MR
      IF ( INNER(I,JJ) .EQ. 1 )  GO TO 30
   20 CONTINUE
      JJ = MR
   30 JTOP = JJ
      J2 = J - 1
      DO 40 JJ = 1,J2
      JJJ = J - JJ
      IF ( INNER(I,JJJ) .EQ. 1 ) GO TO 50
   40 CONTINUE
   50 JBOT = JJJ
C-->     COON'S BLENDING FUNCTIONS
      U = FLOAT(I-IBACK)/FLOAT(IFRONT-IBACK)
      V = FLOAT(J-JBOT)/FLOAT(JTOP-JBOT)
C        F0U GOES FROM 0 TO 1 AS U INCREASES
C        F1U GOES FROM 1 TO 0 AS U INCREASES. ETC.
      F0U = F0(U)
      F1U = F1(U)
      F0V = F0(V)
      F1V = F1(V)
      R(I,J) = R(IBACK,J)*F1U + R(IFRONT,J)*F0U + R(I,JBOT)*F1V
     1 + R(I,JTOP)*F0V - R(IBACK,JBOT)*F1U*F1V - R(IFRONT,JBOT)*F0U*F1V
     2 - R(IBACK,JTOP)*F1U*F0V - R(IFRONT,JTOP)*F0U*F0V
      Z(I,J) = Z(IBACK,J)*F1U + Z(IFRONT,J)*F0U + Z(I,JBOT)*F1V
     1 + Z(I,JTOP)*F0V - Z(IBACK,JBOT)*F1U*F1V - Z(IFRONT,JBOT)*F0U*F1V
     2 - Z(IBACK,JTOP)*F1U*F0V - Z(IFRONT,JTOP)*F0U*F0V
   60 CONTINUE
   70 CONTINUE
   80 CONTINUE
   90 CONTINUE
      RETURN
      END
      SUBROUTINE  ELNODE (I,J,NPTNUM,MAXCOL,MAXROW,LTYPE,
     1                    IPLUS,JPLUS,NODES)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     BUILD ELEMENT INCIDENCES LIST FOR ELEMENT AT POINT (I,J)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION NPTNUM(MAXCOL,MAXROW), IPLUS(LTYPE), JPLUS(LTYPE),
     1          NODES(LTYPE)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      DO 5  K = 1,LTYPE
        IK = I + IPLUS(K)
        JK = J + JPLUS(K)
        NODES(K) = NPTNUM(IK,JK)
        IF ( NODES(K) .LE. 0 )  THEN
          WRITE ( NPRT,*) ' STOP, ZERO NODE IN TOPOLOGY'
          CALL  IPRINT (NPTNUM,MAXCOL,MAXROW)
          STOP
        ENDIF
    5 CONTINUE
      RETURN
      END
      SUBROUTINE  FINSHP
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INSTALLATION DEPENDENT PLOT CLOSE
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      EXTERNAL PLOT
C     CALCOMP CONTROL
      WRITE (NPRT,1000 )
 1000 FORMAT (' END PLOTS')
      CALL PLOT (0.0,0.0,999)
      RETURN
      END
      SUBROUTINE FOURPT (MC,MR,R,Z,INNER,IMIN,IMAX,IGIVEN,N1)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C        FOUR POINT BI-LINEAR INTERPOLATION
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(MC,MR), Z(MC,MR), INNER(MC,MR), IMIN(MR),
     1           IMAX(MR), IGIVEN(MC)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      JLESS = MR - 1
      N1 = 1
      IF ( JLESS .LE. 1 )  RETURN
      DO 80  J = 2,JLESS
C        COUNT GIVEN POINTS IN ROW
      NG = 0
      I1 = IMIN(J)
      I2 = IMAX(J)
      DO 10 I = I1,I2
      IF ( INNER(I,J) .NE. 1 )  GO TO 10
      NG = NG + 1
      IGIVEN(NG) = I
   10 CONTINUE
C        SKIP TO NEXT ROW IF ALL POINTS ARE GIVEN
      IF ( NG .EQ. (I2-I1+1) )  GO TO 80
      NGLESS = NG - 1
      DO 70 IG = 1,NGLESS
      IGPLUS = IG + 1
C        FIND LIMITING I POINTS
      IF ( (IGIVEN(IGPLUS)-IGIVEN(IG)) .EQ. 1 )  GO TO 70
      IBACK = IGIVEN(IG)
      IFRONT = IGIVEN(IGPLUS)
      I1 = IBACK + 1
      I2 = IFRONT - 1
      DO 60  I = I1,I2
C        FIND LIMITING J POINTS
      J1 = J + 1
      DO 20  JJ = J1,MR
      IF ( INNER(I,JJ) .EQ. 1 )  GO TO 30
   20 CONTINUE
   30 JTOP = JJ
      J2 = J - 1
      DO 40  JJ = 1, J2
      JJJ = J - JJ
C        CALCULATE COORDINATES
      IF ( INNER(I,JJJ) .EQ. 1 )  GO TO 50
   40 CONTINUE
   50 JBOT = JJJ
      A = IFRONT - I
      B = JTOP - J
      C = I - IBACK
      D = J - JBOT
C        CALCULATE R COORDINATE OF POINT I,J
      R(I,J) = ( R(IFRONT,J)/(A*(A+C)) + R(I,JTOP)/(B*(B+D))
     1 + R(IBACK,J)/(C*(C+A)) + R(I,JBOT)/(D*(D+B)))*A*B*C*D/(A*C+B*D)
C        CALCULATE Z COORDINATE OF POINT I,J
      Z(I,J) = ( Z(IFRONT,J)/(A*(A+C)) + Z(I,JTOP)/(B*(B+D))
     1 + Z(IBACK,J)/(C*(C+A)) + Z(I,JBOT)/(D*(D+B)))*A*B*C*D/(A*C+B*D)
   60 CONTINUE
   70 CONTINUE
   80 CONTINUE
      RETURN
      END
      SUBROUTINE  IJELEM (MC,MR,ALINE,LEMNUM,IMAX,IMIN)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PRINT SKETCH OF GIVEN I-J ELEMENT DIGITS FOR MESH
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  ALINE(MC), LEMNUM(MC,MR), IMAX(MR), IMIN(MR),
     1           ANUM(10)
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA  BLANK, STAR/' ','*'/
      DATA  ANUM /'0','1','2','3','4','5','6','7','8','9'/
      JPLUS = MR + 1
      WRITE (NPOP,5010)
 5010 FORMAT (/,'   J, 1ST ELEM, OTHER ELEMENT NUMBER DIGITS')
      DO 30  J = 1,MR
        JJ = JPLUS - J
        I1 = IMIN(JJ)
        L  = LEMNUM(I1,JJ)
        I2 = IMAX(JJ)
        DO 10  I = 1,MC
   10     ALINE(I) = BLANK
        DO 20 I = I1,I2
          ITEST = LEMNUM(I,JJ)
          ND    = ITEST - ((ITEST/10)*10)
          IF ( ITEST .GT. 0 )  THEN
            ALINE(I) = ANUM(ND + 1)
          ELSE
            ALINE(I) = BLANK
          ENDIF
   20     CONTINUE
        IF ( L .GT. 0 ) THEN
          WRITE (NPOP,5115) JJ, L, ALINE
 5115     FORMAT ( I5, I6, 5X, (99A1) )
        ELSE
          WRITE (NPOP,5120) JJ, ALINE
 5120     FORMAT ( I5, 6X, 5X, (99A1) )
        ENDIF
   30 CONTINUE
      RETURN
      END
      SUBROUTINE  IJMAT (MC,MR,ALINE,LEMNUM,LIMAX,LIMIN,MAT,LINC)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PRINT SKETCH OF ELEMENT MATERIAL CODES
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  ALINE(MC), LEMNUM(MC,MR), LIMAX(MR), LIMIN(MR),
     1           MAT(MC,MR), ANUM(10)
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA  BLANK, STAR/' ','*'/
      DATA  ANUM /'0','1','2','3','4','5','6','7','8','9'/
      JPLUS = MR + 1
      WRITE (NPLM,5010)
 5010 FORMAT (/,'   J, 1ST ELEM, MATERIAL NUMBER DIGITS')
      DO 30  J = 1,MR
        JJ = JPLUS - J
        DO 10  I = 1,MC
   10     ALINE(I) = BLANK
        I1 = LIMIN(JJ)
        L  = LEMNUM(I1,JJ)
        I2 = LIMAX(JJ)
        DO 20 I = I1,I2,LINC
          LI = LEMNUM(I,JJ)
          IF ( LI .GT. O )  THEN
            ITEST = MAT(I,JJ)
            ND    = ITEST - ((ITEST/10)*10)
            IF ( ITEST .GT. 0 )  THEN
              ALINE(I) = ANUM(ND + 1)
            ELSE
              ALINE(I) = BLANK
            ENDIF
          ENDIF
   20     CONTINUE
        IF ( L .GT. 0 ) THEN
          WRITE (NPLM,5115) JJ, L, ALINE
 5115     FORMAT ( I5, I6, 5X, (99A1) )
        ELSE
          WRITE (NPLM,5120) JJ, ALINE
 5120     FORMAT ( I5, 6X, 5X, (99A1) )
        ENDIF
   30 CONTINUE
      RETURN
      END
      SUBROUTINE  IJMESH (MC,MR,ALINE,INNER,IMAX,IMIN,ISTEP)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PRINT SKETCH OF GIVEN I-J POINTS FOR MESH
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  ALINE(MC), INNER(MC,MR), IMAX(MR), IMIN(MR)
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA  BLANK,STAR/' ','*'/
      JPLUS = MR + 1
      IF ( ISTEP .EQ. 1 )  THEN
        WRITE (NPIP,5000)
 5000   FORMAT (/,'    J  IMIN  IMAX     GIVEN I-J POINTS')
      ELSEIF ( ISTEP .EQ. 2 )  THEN
        WRITE (NPIP,5005)
 5005   FORMAT (/,'    J  IMIN  IMAX     GIVEN BC FLAGS')
      ELSEIF ( ISTEP .EQ. 3 )  THEN
        WRITE (NPOP,5010)
 5010   FORMAT (/,'    J  IMIN  IMAX     GENERATED NODES')
      ELSEIF ( ISTEP .EQ. 4 )  THEN
        WRITE (NPIP,5015)
 5015   FORMAT (/,'    J  IMIN  IMAX     GIVEN BC VALUES')
      ELSEIF ( ISTEP .EQ. 5 )  THEN
        WRITE (NPIP,5020)
 5020   FORMAT (/,'    J  IMIN  IMAX     GIVEN FLUX VALUES')
      ENDIF
      DO 30  J = 1,MR
        JJ = JPLUS - J
        I1 = IMIN(JJ)
        I2 = IMAX(JJ)
        DO 10  I = 1,MC
   10     ALINE(I) = BLANK
        DO 20 I = I1,I2
          ITEST = INNER(I,JJ)
          IF ( ITEST .GT. 0 )  ALINE(I) = STAR
   20     CONTINUE
        IF ( ISTEP .EQ. 3 )  THEN
          WRITE (NPOP,5115) JJ,I1,I2,ALINE
 5115     FORMAT ( I5,2I6,5X,(99A1) )
        ELSE
          WRITE (NPIP,5115) JJ,I1,I2,ALINE
        ENDIF
   30 CONTINUE
      RETURN
      END
      SUBROUTINE  IJNODE (MC,MR,ALINE,NPTNUM,IMAX,IMIN)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PRINT SKETCH OF GIVEN I-J NODE DIGITS FOR MESH
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  ALINE(MC), NPTNUM(MC,MR), IMAX(MR), IMIN(MR),
     1           ANUM(10)
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA  BLANK,STAR/' ','*'/
      DATA  ANUM /'0','1','2','3','4','5','6','7','8','9'/
      JPLUS = MR + 1
      WRITE (NPOP,5010)
 5010 FORMAT (/,'   J, 1ST NODE, OTHER NODE NUMBER DIGITS')
      DO 30  J = 1,MR
        JJ = JPLUS - J
        I1 = IMIN(JJ)
        N1 = NPTNUM(I1,JJ)
        I2 = IMAX(JJ)
        DO 10  I = 1,MC
   10     ALINE(I) = BLANK
        DO 20 I = I1,I2
          ITEST = NPTNUM(I,JJ)
          ND    = ITEST - ((ITEST/10)*10)
          IF ( ITEST .GT. 0 )  ALINE(I) = ANUM(ND + 1)
   20     CONTINUE
        WRITE (NPOP,5115) JJ,N1,ALINE
 5115   FORMAT ( I5, I6, 5X, (99A1) )
   30 CONTINUE
      RETURN
      END
      SUBROUTINE  IJQUAD (LTYPE,LN,IQ,JQ)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     FORM (I,J) COORDS OF NODES OF QUADRILATERAL
C     NODE K OF ELEMENT AT I,J IS AT (I+IQ(K),J+JQ(K))
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  IQ(LN), JQ(LN)
      DIMENSION I4(4), J4(4), I8(8), J8(8), I12(12), J12(12),
     1          I9(9), J9(9), I16(16), J16(16), I17(17), J17(17),
     2          I25(25), J25(25)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C      LAGRANGIAN Q4
      DATA I4, J4  / 0,1,1,0           ,   0,0,1,1 /
C      SENDIPITY Q8
      DATA I8, J8  / 0,2,2,0,1,2,1,0   ,   0,0,2,2,0,1,2,1 /
C      LAGRANGIAN Q9
      DATA I9, J9  / 0,2,2,0,1,2,1,0,1 ,   0,0,2,2,0,1,2,1,1 /
C      SENDIPITY Q12
      DATA I12,J12 / 0,3,3,0,1,3,2,0,2,3,1,0 ,
     1               0,0,3,3,0,1,3,2,0,2,3,1 /
C      LAGRANGIAN Q16
      DATA I16,J16 / 0,3,3,0,1,3,2,0,2,3,1,0,1,2,2,1 ,
     1               0,0,3,3,0,1,3,2,0,2,3,1,1,1,2,2 /
C      SENDIPITY Q17
      DATA I17,J17 / 0,4,4,0,1,4,3,0,2,4,2,0,3,4,1,0,2 ,
     1               0,0,4,4,0,1,4,3,0,2,4,2,0,3,4,1,2 /
C      LAGRANGIAN Q25
      DATA I25,J25 / 0,4,4,0,1,4,3,0,2,4,2,0,3,4,1,0,
     1               1,2,3,3,3,2,1,1,2 ,
     2               0,0,4,4,0,1,4,3,0,2,4,2,0,3,4,1,
     3               1,1,1,2,3,3,3,2,2 /
C      LAGRANGIAN Q4
      IF ( LTYPE .EQ. 4 )  THEN    
        DO 4 I = 1,LTYPE
        IQ(I) = I4(I)
   4    JQ(I) = J4(I)
        RETURN
      ENDIF
C      LAGRANGIAN Q9
      IF ( LTYPE .EQ. 9 )  THEN    
        DO 9 I = 1,LTYPE
        IQ(I) = I9(I)
   9    JQ(I) = J9(I)
        RETURN
      ENDIF
C      LAGRANGIAN Q16
      IF ( LTYPE .EQ. 16 )  THEN    
        DO 16 I = 1,LTYPE
        IQ(I) = I16(I)
   16   JQ(I) = J16(I)
        RETURN
      ENDIF
C      LAGRANGIAN Q25
      IF ( LTYPE .EQ. 25 )  THEN    
        DO 25 I = 1,LTYPE
        IQ(I) = I25(I)
   25   JQ(I) = J25(I)
        RETURN
      ENDIF
C      SERINDIPITY Q8
      IF ( LTYPE .EQ. 8 )  THEN    
        DO 8  I = 1,LTYPE
        IQ(I) = I8(I)
   8    JQ(I) = J8(I)
        RETURN
      ENDIF
C      SERINDIPITY Q12
      IF ( LTYPE .EQ. 12 )  THEN    
        DO 12  I = 1,LTYPE
        IQ(I) = I12(I)
   12   JQ(I) = J12(I)
        RETURN
      ENDIF
C      SERINDIPITY Q16
      IF ( LTYPE .EQ.-16 )  THEN    
        DO 216  I = 1,16   
        IQ(I) = I17(I)
  216   JQ(I) = J17(I)
        RETURN
      ENDIF
C      SERINDIPITY Q17
      IF ( LTYPE .EQ. 17 )  THEN    
        DO 17  I = 1,LTYPE
        IQ(I) = I17(I)
   17   JQ(I) = J17(I)
        RETURN
      ENDIF
      WRITE (NPRT,*) ' INVALID QUAD TYPE ', LTYPE
      STOP  
      END
      SUBROUTINE  IJTORZ (MC,MR,R,Z)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     REPLACE R,Z COORDINATES BY CORRESPONDING I,J VALUES
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(MC,MR), Z(MC,MR)
      DO 10  J = 1,MR
        DO 20  I = 1,MC
        R(I,J) = I
   20   Z(I,J) = J
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  IJTRI (LTYPE,IT13,JT13,IT24,JT24)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     SELECT I,J VALUES OF NODES FOR TRIANGLES BY SPLITTING QUAD
C               INTO COMPLETE LAGRANGIAN TRIANGLES
C     NODE K OF ELEMENT AT I,J IS AT (I+IT(K),J+JT(K))
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION I313(6), I324(6), I613(12), I624(12),
     1          I1013(20), I1024(20), I1513(30), I1524(30),
     2          J313(6), J324(6), J613(12), J624(12),
     3          J1013(20), J1024(20), J1513(30), J1524(30)
C     DIMENSION IT(2*LTYPE), JT(2*LTYPE)
      DIMENSION IT13(1), JT13(1), IT24(1), JT24(1)
C     IT#13, JT#13 give increments to nodes when split on 1-3 diagonal
C     IT#24, JT#24 give increments to nodes when split on 2-4 diagonal
C     nodes are numbered ccw. corners 1st, 1st node on each side, then
C       2nd on each side, etc.  Then to interior loop working to center.
      DATA I313, J313 / 0,1,1  ,  0,1,0   ,
     1                  0,0,1  ,  0,1,1 /
      DATA I324, J324 / 0,1,0  ,  1,1,0   ,
     1                  0,0,1  ,  0,1,1 /
      DATA I613, J613 / 0,2,2,1,2,1  ,  0,2,0,1,1,0   ,
     1                  0,0,2,0,1,1  ,  0,2,2,1,2,1 /
      DATA I624, J624 / 0,2,0,1,1,0  ,  2,2,0,2,1,1   ,
     1                  0,0,2,0,1,1  ,  0,2,2,1,2,1 /
      DATA I1013, J1013 / 0,3,3,1,3,2,2,3,1,2 , 0,3,0,1,2,0,2,1,0,1 ,
     1                    0,0,3,0,1,2,0,2,1,1 , 0,3,3,1,3,2,2,3,1,2 /
      DATA I1024, J1024 / 0,3,0,1,2,0,2,1,0,1 , 3,3,0,3,2,1,3,1,2,2 ,
     1                    0,0,3,0,1,2,0,2,1,1 , 0,3,3,1,3,2,2,3,1,2 /
      DATA I1513, J1513 / 0,4,4,1,4,3,2,4,2,3,4,1,2,3,3,
     1                    0,4,0,1,3,0,2,2,0,3,1,0,1,2,1,
     2                    0,0,4,0,1,3,2,2,2,0,3,1,1,1,2,
     3                    0,4,4,1,4,3,2,4,2,3,4,1,2,3,3  /
      DATA I1524, J1524 / 0,4,0,1,3,0,2,2,0,3,1,0,1,2,1,
     1                    4,4,0,4,3,1,4,2,2,4,1,3,3,3,2,
     2                    0,0,4,0,1,3,0,2,2,0,3,1,1,1,2,
     3                    0,4,4,1,4,3,2,4,2,3,4,1,2,3,3  /
      LIMIT = LTYPE + LTYPE
      IF (  LTYPE .EQ. 3 )   THEN    
        DO 3  I = 1,LIMIT
        IT13(I) = I313(I)
        JT13(I) = J313(I)
        IT24(I) = I324(I)
   3    JT24(I) = J324(I)
        RETURN
      ENDIF
      IF (  LTYPE .EQ. 6 )   THEN    
        DO 6  I = 1,LIMIT
        IT13(I) = I613(I)
        JT13(I) = J613(I)
        IT24(I) = I624(I)
   6    JT24(I) = J624(I)
        RETURN
      ENDIF
      IF (  LTYPE .EQ. 10 )   THEN    
        DO 10  I = 1,LIMIT
        IT13(I) = I1013(I)
        JT13(I) = J1013(I)
        IT24(I) = I1024(I)
   10   JT24(I) = J1024(I)
        RETURN
      ENDIF
      IF (  LTYPE .EQ. 15 )   THEN    
        DO 15  I = 1,LIMIT
        IT13(I) = I1513(I)
        JT13(I) = J1513(I)
        IT24(I) = I1524(I)
   15   JT24(I) = J1524(I)
        RETURN
      ENDIF
      WRITE (NPRT,*) ' INVALID TRIANGLE TYPE ', LTYPE
      STOP
      END
      SUBROUTINE  INEL (MC,MR,MAT,LIMIN,LIMAX,ECHO,NMTCD,
     1                  IJ2JI,IJFLAG,LINC,IMIN,IMAX)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INPUT CHANGES IN ELEMENT CODES
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER ECHO
      DIMENSION  MAT(MC,MR), LIMIN(MR), LIMAX(MR),
     1           IMIN(MR), IMAX(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     MAT(I,J) = CODE OF QUAD WITH BOTTOM LEFT CORNER AT I,J
C     IJFLAG   = FLAG TO MARK NODE NOT TO BE NUMBERED
C     IF MAT   < 0, ELEMENT IS PART OF A HOLE
      IERR = 0
      IF ( ECHO .EQ. 0 )  WRITE (NPRT,5000)
 5000 FORMAT ( //,' CHANGES IN BLOCKS OF ELEMENT CODES',/,
     1 ' BLOCK  MAT.NO.  LIMITS:   I1   J1   I2   J2',/)
C             ..........< I2,J2   MM IS ASSIGNED TO ALL ELEMENT
C             .        .          CODES FOR ELEMENTS WITH THEIR
C             .        .          BOTTOM LHC WITHIN THE BLOCK
C      I1,J1 >..........  MAT(I,J)=MM FOR (I1,J1)=<(I,J)<(I2,J2)
      DO 80  N = 1,NMTCD
        READ (NCRD,5010)  I1, J1, I2, J2, MM
 5010   FORMAT ( 5I3 )
        IF ( I2 .EQ. 0 ) I2 = I1
        IF ( J2 .EQ. 0 ) J2 = J1
        IF ( IJ2JI .EQ. 0 )  CALL SWAP(I1,J1)
        IF ( IJ2JI .EQ. 0 )  CALL SWAP(I2,J2)
        IF ( ECHO .EQ. 0 )  WRITE (NPRT,5020) N, MM, I1, J1, I2, J2
 5020   FORMAT (I5, I8, 11X, 4I5)
        IF ( I1.LT.1 .OR. I2.GT.MC .OR. 
     1       J1.LT.1 .OR. J2.GT.MR )   WRITE (NPRT,5040)
 5040   FORMAT (' WARNING, NEXT DATA HAVE NONEXISTANT ELEMENTS')
        DO 50  J = J1,J2-1
          IF ( I1 .LT. LIMIN(J) .OR.
     1         I2 .GE. LIMAX(J) )  WRITE (NPRT,*) 
     2         ' NONEXISTANT ELEMENTS AT J =', J
          DO 40  I = I1,I2-1
            IF ( MM .LT. 1 )  THEN
C            * FLAG POINTS AND ELEMENTS TO SKIP FOR A HOLE *
              IF ( J.GT.J1 .AND. J.LT.J2 )  THEN                   
C               LEFT GRID EDGE
                IF ( I.EQ.IMIN(J) )            MAT(I,J) = IJFLAG
C               RIGHT GRID EDGE
                IF ( I.EQ.IMAX(J) )            MAT(I,J) = IJFLAG
C               INTERIOR TO HOLE
                IF ( I.GT.I1 .AND. I.LT.I2 )   MAT(I,J) = IJFLAG
              ENDIF
C             BOTTOM OR TOP ROW OF GRID
              IF ( I.GT.I1 .AND. I.LT.I2 .AND.                
     1           ( J.EQ.1  .OR.  J.EQ.MR ))    MAT(I,J) = IJFLAG
C             LEFT OR BOTTOM SIDE OF HOLE
              IF ( I.EQ.I1 .AND. J.LT.J2 )     MAT(I,J) = MM
              IF ( I.LT.I2 .AND. J.EQ.J1 )     MAT(I,J) = MM
            ELSE
C             * REAL MATERIAL *
              MAT(I,J) = MM
            ENDIF
   40     CONTINUE
   50   CONTINUE
        WRITE (NPRT,*) '---- END OF MATERIAL SET ', N
   80 CONTINUE
C      CHECK FOR NODES DELETED ON THE "EDGES"
      DO 90  J = 1,MR
        LJ  = LIMIN(J)
        IF ( (MAT(LJ,J) .EQ. -1 ) .AND. 
     1       (MAT(LJ+1,J) .EQ. IJFLAG ) )  THEN
          LIMIN(J) = LIMIN(J) + LINC
          IMIN(J)  = IMIN(J) + LINC
        ENDIF
        LJ  = LIMAX(J)
        IF ( (MAT(LJ,J) .EQ. IJFLAG ) .AND. 
     1       (MAT(LJ+1,J) .EQ. IJFLAG ) )  THEN
          LIMAX(J) = LIMAX(J) - LINC
          IMAX(J)  = IMAX(J) - LINC
        ENDIF
   90 CONTINUE
      RETURN
c     IF ( MM .LT. 1 )  MAT(I2,J1) = MM
      END
      SUBROUTINE  INIJ (I1,I2,I3,J1,J2,J3,R1,R2,R3,Z1,Z2,Z3,CC,
     1                  MC,MR,R,Z,CODE,INNER,IMIN,IMAX,INTEST,
     2                  ECHO,POLAR)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INPUT CURVES OR LINES IN I-J AND R-Z SPACES
C     I1---I3---I2   ===TO==>>   R1---R3---R2
C     J1   J3   J2               Z1   Z3   Z2
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER ECHO, POLAR
      DIMENSION R(MC,MR), Z(MC,MR), CODE(MC,MR), INNER(MC,MR),
     1          IMIN(MR), IMAX(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      IF ( I2 .EQ. 0 )  I2 = I1
      IF ( J2 .EQ. 0 )  J2 = J1
C        FIND BIGGEST LOOP REQUIRED
      LOOP = MAX0 ( IABS(I2-I1), IABS(J2-J1) ) + 1
      DL = LOOP - 1
      IF ( LOOP .EQ. 1 )  DL = 1.0
C        LOOP OVER POSSIBLE I,J PAIRS
      A1 = I1
      A2 = I2
      A3 = I3
      IF ( I3 .EQ. 0 )  A3 = FLOAT(I1+I2)/2.
      B1 = J1
      B2 = J2
      B3 = J3
      IF ( J3 .EQ. 0 )  B3 = FLOAT(J1+J2)/2.
      DO 10  L = 1, LOOP
C         COMPUTE UNIT COORDINATE
        X = FLOAT(L-1)/DL
C         LOCAL QUADRATIC INTERPOLATIONS
        H1 = (X-1.)*(X+X-1.)
        H2 = X*(X+X-1.)
        H3 = 4.*X*(1.-X)
C        INTERPOLATE IN I,J R,Z
        I = A1*H1 + A2*H2 + A3*H3 + 0.1
        IF ( I .LT. 1 )   I = 1
        IF ( I .GT. MC )  I = MC
        J = B1*H1 + B2*H2 + B3*H3 + 0.1
        IF ( J .LT. 1 )   J = 1
        IF ( J .GT. MR )  J = MR
C        RECOVER OLD BC FLAG ?
        IF ( INTEST .GT. 2 ) CC = CODE(I,J)
C        CIRCULAR ARC
        IF ( POLAR .GT. 0 )  THEN    
          RAD = R1*(1.-X) + R2*X
          ANG = Z1*(1.-X) + Z2*X
          R(I,J) = RAD*COS(ANG/57.29578) + R3
          Z(I,J) = RAD*SIN(ANG/57.29578) + Z3
          IF ( ECHO .EQ. 0 ) WRITE (NPRT,5080)
     1                       I, J, CC, RAD, ANG, R3, Z3
 5080     FORMAT ( 2I4, F7.0, 1PE14.5, 3 (1PE16.5) )
        ELSE
C          GENERAL PARABOLIC CHANGE OR LINEAR CHANGE
          R(I,J) = R1*H1 + R2*H2 + R3*H3
          Z(I,J) = Z1*H1 + Z2*H2 + Z3*H3
          IF ( ECHO .EQ. 0 ) WRITE (NPRT,5080)
     1                       I, J, CC, R(I,J), Z(I,J)
        ENDIF
C                FLAG OTHER DATA
        INNER(I,J) = 1
        CODE(I,J) = CC
        IF ( I .LT. IMIN(J) )  IMIN(J) = I
        IF ( I .GT. IMAX(J) )  IMAX(J) = I
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  INITAL (MR, MC, MAT, IMIN, LIMIN)
C     - - - - - - - - - - - - - - - - - - - - - - -
C           INITIALIZE VALUES OF ARRAYS
C     - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  MAT(MC,MR), IMIN(MR), LIMIN(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     MAT   = ELEMENT MATERIAL CODE
C     IMIN  = MIN I FOR ANY NODE ON ROW J
C     LIMIN = MIN I FOR ANY ELEMENT IN ROW J
C      SET INITIAL MINIMUMS = MAXIMUMS
      DO 30  I = 1, MR
        IMIN(I) = MC
   30 LIMIN(I)  = MC
C      SET MAT=1, EXCEPT TOP ROW
      DO 60 J = 1, MR
        DO 40 I = 1, MC-1  
   40   MAT(I,J) = 1
        DO 50 I = 1,MC    
   50   MAT(I,MR) = -1
   60 CONTINUE
      RETURN
      END
      SUBROUTINE  INPTS (MC,MR,R,Z,CODE,INNER,IMIN,IMAX,
     1                   LIMIN,LIMAX,ECHO,NPTCD,MAP,IJ2JI)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INPUT COORDINATES AND NODAL POINT CODES
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER ECHO, POLAR
      DIMENSION  R(MC,MR), Z(MC,MR), CODE(MC,MR), INNER(MC,MR),
     1            IMIN(MR), LIMIN(MR), LIMAX(MR), IMAX(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
C     R(I,J),Z(I,J)     = COORDS OF PT CORRESPONDING TO (I,J)
C     CODE(I,J)         = PT BOUNDARY CODE ASSIGNED TO (I,J)
C     INNER(I,J)        > 0 IMPLIES PT (I,J) WAS INPUT
C     IMIN(J),IMAX(J)   = MIN  MAX I FOR ANY PT IN ROW J
C     LIMIN(J),LIMAX(J) = MIN, MAX I FOR A PT IN AN ELEM IN ROW J
      IF ( ECHO .EQ. 0 )  WRITE (NPRT,5000)
 5000   FORMAT ( //,' NODAL POINT INPUT DATA',/,
     1  '   I   J    IBC   X OR R          Y OR ANGLE',
     2  '      X-CENTER       Y-CENTER')
      DO 160  N = 1,NPTCD
      ISKIP = 0
      R1 = 0.
      R2 = 0.
      R3 = 0.
      Z1 = 0.
      Z2 = 0.
      Z3 = 0.
      CC = 0.
C==>>  MEANING IF INTEST:
C      =0, SINGLE I,J AND R,Z POINT
C      =1, LINE IN I,J TO LINE IN R,Z
C      =2, LINE OR ARC IN I,J LINE OR ARC IN R,Z
C      =3, LINE IN I,J TO CONNECT PREVIOUS R,Z ENDS
C      =4, LINE OR ARC IN I,J TO CONNECT PREVIOUS R,Z POINTS
C     INTEST   REQUIRED ?  (Y-YES, N-NO, O-OPTIONAL, V-VOID)
C     VALUE  I1,J1  I2,J2  I3,J3  R1,Z1  R2,Z2  R3,Z3 POLAR
C     0        Y      N      N      Y      N      N     O
C     1        Y      Y      N      Y      Y      N     O
C     2        Y      Y      O      Y      Y      O     V
C     3        Y      Y      N      N      N      N     V
C     4        Y      Y      O      N      N      N     V
C     IF POLAR = ZERO THIS INDICATES INPUT IS IN CARTESIAN COORDS.
C       R1,Z1 IS A CARTESIAN COORD PAIR ( ETC. FOR R2,Z2 AND R3,Z3)
C     IF POLAR .NE. ZERO  THIS INDICATES INPUT IS IN POLAR COORDS.
C       RELATIVE TO A CIRCLE WHOSE CENTER IS (RC,ZC)
C       R1,Z1 = INITIAL RADIUS & ANGLE,  R2,Z2 = FINAL RADIUS & ANGLE
C       AND R3,Z3 = CENTER OF CIRCLE
      READ (NCRD,5010) INTEST, I1, J1, I2, J2, I3, J3, POLAR
 5010 FORMAT ( 8I4 )
      IF ( INTEST .GT. 1 .AND. POLAR .NE. 0 ) 
     1  WRITE (NPRT,*) 'POLAR OPTION INVALID, IGNORED.'
      IF ( INTEST .EQ. 0 )  THEN
        I2 = I1
        I3 = I1
        J2 = J1
        J3 = J1
      ENDIF
      IF ( IJ2JI .EQ. 0 )  THEN    
        CALL SWAP(I1,J1)
        CALL SWAP(I2,J2)
        CALL SWAP(I3,J3)
      ENDIF
      IF ( INTEST .LT. 3 )  
     1     READ (NCRD,5011) R1, Z1, R2, Z2, R3, Z3, CC
 5011      FORMAT ( 7F10.0 )
 2    IF ( INTEST .EQ. 1 .AND. POLAR .EQ. 0 )  THEN	
C        STRAIGHT LINE IN R-Z
        R3 = (R1 + R2)*0.5
        Z3 = (Z2 + Z1)*0.5
      ENDIF
 6    IF ( INTEST .EQ. 3 )  THEN   
C        TWO POINTS MUST EXIST
        IF ( (INNER(I1,J1)*INNER(I2,J2)) .EQ. 0 )  THEN
          WRITE (NPRT,*) 'ERROR, POINT NOT DEFINED, SKIP'
          ISKIP = 1
        ELSE
          R1 = R(I1,J1)
          Z1 = Z(I1,J1)
          R2 = R(I2,J2)
          Z2 = Z(I2,J2)
C          STRAIGHT LINE IN R-Z
          R3 = (R1+R2)*0.5
          Z3 = (Z1+Z2)*0.5
        ENDIF
      ENDIF
 8    IF ( INTEST .EQ. 4 )  THEN
C        THREE POINTS MUST EXIST
        IF ((INNER(I1,J2)*INNER(I2,J2)*INNER(I3,J3)).EQ.0) THEN
          WRITE (NPRT,*) 'ERROR, POINT NOT DEFINED, SKIP'
          ISKIP = 1
        ELSE
C        POINTS PREVIOUSLY GIVEN
          R1 = R(I1,J1)
          Z1 = Z(I1,J1)
          R2 = R(I2,J2)
          Z2 = Z(I2,J2)
          R3 = R(I3,J3)
          Z3 = Z(I3,J3)
        ENDIF
      ENDIF
      IF ( ISKIP .EQ. 0 ) CALL INIJ (I1,I2,I3,J1,J2,J3,
     1                    R1,R2,R3,Z1,Z2,Z3,CC,
     2                    MC,MR,R,Z,CODE,INNER,IMIN,IMAX,
     3                    INTEST,ECHO,POLAR)
      IF ( ECHO .EQ. 0 )  WRITE (NPRT,*) 
     1     '------ END OF POINT GENERATION CARD ', N
 160  CONTINUE
C        FIND BEGINNING AND END OF POINTS AND ELEMENTS FOR EACH J
      IF ( MR .GT. 2 )  THEN
        J2 = MR - 1
        DO 170  J = 1,J2
          LIMIN(J) = MAX0 ( IMIN(J), IMIN(J+1) )
  170   LIMAX(J) = MIN0 ( IMAX(J), IMAX(J+1) )
      ENDIF
      LIMIN(MR) = IMIN(MR)
      LIMAX(MR) = IMAX(MR)
C        DETERMINE MAXIMUM NODAL CONNECTIVITY
      MAX = 0
      DO 180  J = 1,J2
        INTEST = IMAX(J) - IMIN(J)
        IF ( INTEST .GT. MAX )  MAX = INTEST
  180 CONTINUE
      MAX = MAX + 2
      IF ( ECHO.EQ.0 )  WRITE (NPRT,5120)  MAX
      IF ( ECHO.EQ.0 )  WRITE (NPIO,5120)  MAX
 5120 FORMAT ( /,' NODAL POINT CONNECTIVITY: ',/,
     1 ' MAXIMUM DIFFERENCE IN NODE NUMBERS, ND =',I3,/,
     2 ' HALF BANDWIDTH = (ND+1)*(DEGREES OF FREEDOM PER NODE)',/)
      RETURN
      END
      SUBROUTINE  IPRINT (M,NR,NC)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
C      PRINTING OF AN INTEGER ARRAY M(NR,NC)
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
      PARAMETER ( LIMIT = 12 )
      DIMENSION  M(1), NCOL(LIMIT)
      DO 40  J = 1,NC,LIMIT
      JL1 = J - 1
      MAXCOL = 1
      K = NC - JL1
      MAXCOL = MIN0 (K,LIMIT)
      MXCLL1 = MAXCOL - 1
      DO 10  L = 1,MAXCOL
   10 NCOL(L) = L + JL1
      WRITE (NPRT,5000)  ( NCOL(N),N=1,MAXCOL )
 5000 FORMAT (' J,I:', 12I5 )  
      DO 30  N = 1,NR
        NL = N + (J-1)*NR
        NH = NL + MXCLL1*NR
        WRITE (NPRT,5010)  N,( M(I),I=NL,NH,NR )
 5010   FORMAT ( I5, 12I5 )
   30   CONTINUE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE LAPLAC (K1,K2,R,Z,INNER,IMIN,IMAX,IGIVEN,N1)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     MESH BY LAPLACE'S EQUATION
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(K1,K2), Z(K1,K2), INNER(K1,K2),
     1           IMIN(K2), IMAX(K2), IGIVEN(K1)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     ITER IS THE MAXIMUM NUMBER OF ITERATIONS ALLOWED
C     OMEGA IS A FACTOR USED TO OPTIMIZE THE NUMBER OF ITERATIONS
      DATA OMEGA, SMALL, TEST, ITER, NMIN / 1.85, 1.E-7, 0.05, 30, 5/
C        USE A FOURPT MAP FOR THE FIRST ITERATION
      CALL FOURPT (K1,K2,R,Z,INNER,IMIN,IMAX,IGIVEN,N1)
C        PREFORM NMIN ITERATIONS BEFORE CHECKING ACCURACY
      J2 = K2 - 1
      DO 40  N = 1,ITER
      N1 = N
      RESID = 0.0
      DO 30  J = 2,J2
      I1 = IMIN(J) + 1
      I2 = IMAX(J) - 1
C        SKIP TO NEXT ROW IF THERE ARE NO INTERIOR POINTS
      IF ( I2 .LT. I1 )  GO TO 30
      DO 20 I = I1,I2
      IF ( INNER(I,J) .NE. 0 )  GO TO 10
      DR = 0.25*( R(I,J-1)+R(I+1,J)+R(I,J+1)+R(I-1,J) ) - R(I,J)
      DZ = 0.25*( Z(I,J-1)+Z(I+1,J)+Z(I,J+1)+Z(I-1,J) ) - Z(I,J)
      R(I,J) = R(I,J) + OMEGA*DR
      Z(I,J) = Z(I,J) + OMEGA*DZ
      RESID = RESID + ABS(DR) + ABS(DZ)
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
      IF ( N .LT. NMIN )  GO TO 40
      IF ( N .EQ. NMIN )  RESID1 = RESID
      RATIO =  RESID/(RESID1 + SMALL)
      IF ( RATIO .LE. TEST )  RETURN
   40 CONTINUE
      RETURN
      END
      SUBROUTINE  LCOUNT (MC,MR,R,Z,LIMAX,LIMIN,MAT,NPTNUM,LTYPE,
     1                    LINC,MINEL,OUTEL,IOPT,LEMENT)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     COUNT AND OUTPUT GENERATED ELEMENTS
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER OUTEL
      DIMENSION  R(MC,MR), Z(MC,MR), LIMAX(MR), LEMENT(MC,MR),
     1           LIMIN(MR), MAT(MC,MR), NPTNUM(MC,MR)
      DIMENSION IQ(25), JQ(25), I13(30), J13(30), I24(30),
     1          J24(30), NODES(30)
      COMMON / IOUNIT / NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      EQUIVALENCE (I13(1),IQ(1)), (J13(1),JQ(1))
C      OVERWRITE THE OLD IGIVEN ARRAY TO RECEICE ELEMENT NUMBERS
      DO 10 J = 1,MR
        DO 20  I = 1,MC
 20     LEMENT(I,J) = -2
 10   CONTINUE
      L  = MINEL - 1
      LN = IABS( LTYPE )
      M  = LN + 1
C        EXTRACT I,J VALUES OF NODES FOR TRIANGLES
      IF ( LTYPE .EQ. 3 .OR. LTYPE .EQ. 6 .OR. LTYPE .EQ. 10 
     1                  .OR. LTYPE .EQ. 15 ) THEN
        CALL  IJTRI (LTYPE,I13,J13,I24,J24)
      ELSE
C       FOR QUADRILATERALS
        CALL  IJQUAD  (LTYPE,LN,IQ,JQ)
      ENDIF
      JLESS = MR - LINC
      DO 50  J = 1,JLESS,LINC
      I1 = LIMIN(J)
      I2 = LIMAX(J) - LINC
      DO 40  I = I1,I2,LINC
C        DO NOT NUMBER ELEMENTS IN A HOLE
      IF (  MAT(I,J) .LT. 0 )  GO TO 40
      MIJ = MAT(I,J)
      L = L + 1
      LEMENT(I,J) = L
C        CHECK AREAS
      IL = I + LINC
      JL = J + LINC
      A123 = R(I,J)*(Z(IL,J)-Z(IL,JL)) + R(IL,J)*(Z(IL,JL)-Z(I,J))
     1     + R(IL,JL)*(Z(I,J)-Z(IL,J))
      A134 = R(I,J)*(Z(IL,JL)-Z(I,JL)) + R(IL,JL)*(Z(I,JL)-Z(I,J))
     1     + R(I,JL)*(Z(I,J)-Z(IL,JL))
      IF ( A123 .LE. 0. .OR. A134 .LE. 0. )  WRITE (NPRT,*)  
     1         ' POSSIBLE NEGATIVE AREA NEAR ELEMENT ', L
      IF ( LTYPE .EQ. 3 .OR. LTYPE .EQ. 6 .OR. LTYPE .EQ. 10 .OR.
     1     LTYPE .EQ. 15 )  THEN                                   
C           TRIANGLES
      D13SQ = (R(I,J)-R(IL,JL))**2+(Z(I,J)-Z(IL,JL))**2
      D24SQ = (R(I,JL)-R(IL,J))**2+(Z(I,JL)-Z(IL,J))**2
      D24SQ = D24SQ*1.2
      IF ( D24SQ .LT. D13SQ ) GO TO 70
C            USE 1-3 DIAGONAL OF QUAD
      CALL ELNODE (I,J,NPTNUM,MC,MR,LN,I13,J13,NODES)
      CALL  LOUT (L,I,J,LN,OUTEL,IOPT,NODES,MIJ,LTYPE)
      L = L + 1
      CALL ELNODE (I,J,NPTNUM,MC,MR,LN,I13(M),J13(M),NODES)
      CALL  LOUT (L,I,J,LN,OUTEL,IOPT,NODES,MIJ,LTYPE)
      GO TO 40
C            USE SHORTER 2-4 DIAGONAL OF QUAD
   70 CALL ELNODE (I,J,NPTNUM,MC,MR,LN,I24,J24,NODES)
      CALL  LOUT (L,I,J,LN,OUTEL,IOPT,NODES,MIJ,LTYPE)
      L = L + 1
      CALL ELNODE (I,J,NPTNUM,MC,MR,LN,I24(M),J24(M),NODES)
      CALL  LOUT (L,I,J,LN,OUTEL,IOPT,NODES,MIJ,LTYPE)
      ELSE
C           QUADRILATERALS
      CALL ELNODE (I,J,NPTNUM,MC,MR,LN,IQ,JQ,NODES)
      CALL  LOUT (L,I,J,LN,OUTEL,IOPT,NODES,MIJ,LTYPE)
      ENDIF
   40 CONTINUE
   50 CONTINUE
      WRITE (NPRT,*) ' MAX ELEMENT NUMBER =', L
      WRITE (NPIO,*) ' MAX ELEMENT NUMBER =', L
      RETURN
      END
      SUBROUTINE LINETO (R,Z)
C     - - - - - - - - - - - - - - - - - - - -
C     MOVE TO A POINT, PEN DOWN
C     - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     CALCOMP CONTROL
      COMMON /CALCOM/ FIRSTR, FIRSTZ, DELTA
      EXTERNAL PLOT
      IF ( DELTA .NE. 0.0 )  THEN
        X = (R-FIRSTR)/DELTA
        Y = (Z-FIRSTZ)/DELTA
        CALL PLOT (X,Y,2)
      ELSE
        STOP 'ZERO DELTA IN LINETO'
      ENDIF
      RETURN
      END
      SUBROUTINE  LOUT (L,I,J,M,OUTEL,IOPT,NODES,MAT,LTYPE)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     OUTPUT ELEMENT INCIDENCES LIST AND INTEGER CODE
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER OUTEL
      DIMENSION NODES(M)
      DIMENSION L4(4), L8(8), L12(12), L17(17), LPAF(25)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA KALL /1/
C        PAFEC CONVERSION OF QUADS DATA
      DATA  L4  / 1,2,4,3 /
      DATA  L8  / 1,2,4,3, 5,8,6,7 /
      DATA  L12 / 1,2,4,3, 5,9 ,12,8, 6,10, 11,7 /
      DATA  L17 / 1,2,4,3, 5,9,13, 16,12,8, 6,10,14, 15,11,7, 17 /
C        PAFEC QUAD CODES
      DATA LPAF / 1,2,36100,36200,5,36110,7,36210,9,
     1            36120,11,36220,13,14,15,16,36230,18,
     2            19,20,21,22,23,24,25 /
      IF ( KALL .EQ. 0 )  GO TO 5
C        PRINT HEADING ON FIRST CALL
      KALL = 0
      IF ( OUTEL .LT. 2 )  WRITE (NPRT,1000) (K, K = 1,M)
 1000 FORMAT (//,'  E L E M E N T    O U T P U T',/,
     1 ' ELEMENT   (I,J)    MAT', 9('   N', I2),
     2 (/, 23X, 9('   N', I2)) )
      IF ( IOPT .EQ. 2 )  THEN
        WRITE (NPUN,500)  LPAF(M)
        WRITE (NPEL,500)  LPAF(M)
  500   FORMAT ('elements.from.mesh90', /, 'elem=', I5, /,
     1          'number property topology' )
      ENDIF
    5 CONTINUE
C        PRINT ELEMENT DATA
      IF ( OUTEL .LT. 3 )  WRITE(NPRT,1010)
     1                     L,I,J,MAT,(NODES(K),K=1,M)
 1010 FORMAT (I6,I7,',',I3,I6, 9(I6), (/,23X, 9I6) )
      IF ( OUTEL .EQ. 0 .OR. OUTEL .EQ. 3 )  RETURN
C        PUNCH SPECIAL FORMATTED FILES
      IF ( IOPT .EQ. 0 )  THEN    
C-->    MODEL CODE FORMATS
        WRITE (NPUN,1020)  L,NODES
        WRITE (NPEL,1020)  L,NODES
 1020   FORMAT ( 16(I5) )
        RETURN
      ENDIF
      IF ( IOPT .EQ. 1 )  THEN    
C-->    FINITE CODE FORMATS
        WRITE (NPUN,1040)  L,MAT,NODES
        WRITE (NPEL,1040)  L,MAT,NODES
 1040   FORMAT ( 16(I5) )
        RETURN
      ENDIF   
      IF ( IOPT .EQ. 2 )  THEN    
C-->    PAFEC CODE FORMATS
 1050   FORMAT ( 16(I5) )
C              TRIANGLES
        IF ( LTYPE .EQ. 3 .OR. LTYPE .EQ. 6 .OR. LTYPE .EQ. 15 ) THEN
          WRITE (NPUN,1050)  L,MAT,NODES
          WRITE (NPEL,1050)  L,MAT,NODES
        ENDIF
        IF ( LTYPE .EQ. 10 ) THEN
          WRITE (NPUN,1050) L,MAT,( NODES(K), K = 1,9)
          WRITE (NPEL,1050) L,MAT,( NODES(K), K = 1,9)
        ENDIF
C            QUADRILATERALS
        IF ( LTYPE .EQ. 9 .OR. LTYPE .EQ. 16 .OR. LTYPE .EQ. 25 )
     1    WRITE (NPRT,*) ' NO PAFEC TYPE', LTYPE
        IF ( LTYPE .EQ. 4  )  THEN
          WRITE (NPUN,1050) L,MAT,( NODES( L4(K) ), K = 1,4)
          WRITE (NPEL,1050) L,MAT,( NODES( L4(K) ), K = 1,4)
        ENDIF
        IF ( LTYPE .EQ. 8  )  THEN
          WRITE (NPUN,1050) L,MAT,( NODES( L8(K) ), K = 1,8)
          WRITE (NPEL,1050) L,MAT,( NODES( L8(K) ), K = 1,8)
        ENDIF
        IF ( LTYPE .EQ. 12  )  THEN                               
          WRITE (NPUN,1050) L,MAT,( NODES( L12(K) ), K = 1,12)
          WRITE (NPEL,1050) L,MAT,( NODES( L12(K) ), K = 1,12)
        ENDIF
        IF ( LTYPE .EQ. 17  )  THEN                               
          WRITE (NPUN,1050) L,MAT,( NODES( L17(K) ), K = 1,17)
          WRITE (NPEL,1050) L,MAT,( NODES( L17(K) ), K = 1,17)
        ENDIF
        RETURN
      ENDIF
      IF ( IOPT .EQ. 3 )  THEN    
C-->    BALA CODE FORMATS
        WRITE (NPUN,1050)  L, NODES
        WRITE (NPEL,1050)  L, NODES
        WRITE (NPLM,1050)  L, MAT  
        RETURN
      ENDIF
C        ANOTHER CODE
      WRITE (NPRT,*) 'INVALID IOPT, LOUT. STOP.'
      STOP
      END
      SUBROUTINE MOVETO (R,Z)
C     - - - - - - - - - - - - - - - - -
C     MOVE TO A POINT, PEN UP
C     - - - - - - - - - - - - - - - - -
C     CALCOMP CONTROL
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON /CALCOM/ FIRSTR, FIRSTZ, DELTA
      EXTERNAL PLOT
      IF ( DELTA .NE. 0.0 )  THEN
        X = (R-FIRSTR)/DELTA
        Y = (Z-FIRSTZ)/DELTA
        CALL PLOT (X,Y,3)
      ELSE
        STOP 'ZERO DELTA IN MOVETO'
      ENDIF
      RETURN
      END
      SUBROUTINE OMITP (MAXCOL,MAXROW,MAT,LTYPE,LINC,IJFLAG,R,Z)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     FLAG NODES TO BE OMITTED IN HIGH ORDER ELEMENTS
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION R(MAXCOL,MAXROW), Z(MAXCOL,MAXROW),
     1          MAT(MAXCOL,MAXROW)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      IF ( LTYPE .LT. 8 )  RETURN
      I2 = MAXCOL - LINC
      J2 = MAXROW - LINC
      DO 10  I = 1,I2,LINC
        IL = I + LINC
        DO 20  J = 1,J2,LINC
        IF ( MAT(I,J) .LT. 0 )  GO TO 20
        JL = J + LINC
        IF ( LTYPE .EQ. 8 )   MAT(I+1,J+1) = IJFLAG
        IF ( LTYPE .EQ. 12 )  THEN    
          MAT(I+2,J+2) = IJFLAG
          MAT(I+1,J+2) = IJFLAG
          MAT(I+2,J+1) = IJFLAG
          MAT(I+1,J+1) = IJFLAG
        ENDIF
        IF ( LTYPE .EQ.-16 )  THEN    
          MAT(I+3,J+1) = IJFLAG
          MAT(I+3,J+2) = IJFLAG
          MAT(I+3,J+3) = IJFLAG
          MAT(I+2,J+1) = IJFLAG
          MAT(I+2,J+2) = IJFLAG  
          MAT(I+2,J+3) = IJFLAG
          MAT(I+1,J+1) = IJFLAG
          MAT(I+1,J+2) = IJFLAG
          MAT(I+1,J+3) = IJFLAG
        ENDIF
          IF ( LTYPE .EQ. 17 )  THEN    
          MAT(I+3,J+1) = IJFLAG
          MAT(I+3,J+2) = IJFLAG
          MAT(I+3,J+3) = IJFLAG
          MAT(I+2,J+1) = IJFLAG
          MAT(I+2,J+2) = MAT(I,J)
          MAT(I+2,J+3) = IJFLAG
          MAT(I+1,J+1) = IJFLAG
          MAT(I+1,J+2) = IJFLAG
          MAT(I+1,J+3) = IJFLAG
        ENDIF
   20   CONTINUE
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  PBCOUT (MC,MR,NPTS,CODE,IOPT)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PUNCH NODAL POINT B.C. FLAGS AND DEFAULTS
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  NPTS(MC,MR), CODE(MC,MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DIMENSION KODES(6)
      ZERO  = 0.0
      MAX   = 0
      KOUNT = 0
      NG    = 1
      IF ( IOPT .EQ. 2 )  THEN
C       PAFEC TITLE
        WRITE (NPUN,500)
        WRITE (NPBC,500)
  500   FORMAT ('restraints',/,'node direction value')
      ENDIF
      DO 1 J = 1,MC
        DO 2 I = 1,MR
          IF ( CODE(I,J) .LE. 0.0 )  GO TO 2
          NGTEST = CODE(I,J) + 0.1
          IF ( NGTEST .GT. MAX )  MAX = NGTEST
 2      CONTINUE
 1    CONTINUE
      IF ( MAX .GT. 9 )     NG = 2
      IF ( MAX .GT. 99 )    NG = 3
      IF ( MAX .GT. 999 )   NG = 4
      IF ( MAX .GT. 9999 )  NG = 5
      IF ( MAX .GT. 99999 ) NG = 6
      DO 5  J = 1,MC
        DO 10  I = 1,MR
C       DOES NODE EXIST
        IF ( NPTS(I,J) .EQ. 0 )  GO TO 10
        NCODE = CODE(I,J) + 0.1
        IF ( NCODE .EQ. 0 )  GO TO 10
        NPTNUM = NPTS(I,J)
        CALL  PTCODE (NPTNUM,NG,NCODE,KODES)
        DO 20 K = 1, NG
          IF ( KODES(K) .EQ. 0 )  GO TO 20
          KOUNT = KOUNT + 1
          IF ( IOPT  .EQ.  0 ) THEN
C           MODEL  OPTIONS
            WRITE (NPUN,5040)  NPTNUM, K, ZERO
            WRITE (NPBC,5040)  NPTNUM, K, ZERO
 5040       FORMAT ( I5, I5, F5.0 )
          ENDIF
          IF ( IOPT .EQ. 1 .OR. IOPT .EQ. 2 ) THEN
C           FINITE, PAFEC OPTIONS
            WRITE (NPUN,5050)  NPTNUM, K, ZERO
            WRITE (NPBC,5050)  NPTNUM, K, ZERO
 5050       FORMAT ( 2I10,  F5.0 )
          ENDIF
          IF ( IOPT  .EQ.  3 ) THEN
C           FLOW90 OPTIONS
            WRITE (NPUN,5060)  NPTNUM, K, ZERO
            WRITE (NPBC,5060)  NPTNUM, K, ZERO
 5060       FORMAT ( I5, I5, F5.0 )
          ENDIF
   20   CONTINUE
   10 CONTINUE
    5 CONTINUE
      WRITE (NPRT,*) ' NUMBER OF BC FLAGS CREATED =', KOUNT
      WRITE (NPIO,*) ' NUMBER OF BC FLAGS CREATED =', KOUNT
      RETURN
      END
      SUBROUTINE  PCOUNT (MC,MR,R,Z,CODE,MAT,NPTNUM,IMIN,IMAX,
     1                    IJFLAG,IOPT,OUTPT,MINNPT,LINC)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     NUMBER AND OUTPUT NODAL COORDINATES AND POINT CODES
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER  OUTPT
      DIMENSION  R(MC,MR), Z(MC,MR), CODE(MC,MR), MAT(MC,MR),
     1           NPTNUM(MC,MR), IMIN(MR), IMAX(MR)
      COMMON / IOUNIT / NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      NODE = MINNPT - 1
      DO 70  J = 1,MR
        I1 = IMIN(J)
        I2 = IMAX(J)
        DO 60 I = I1,I2
C          CHECK FOR OMITTED POINTS IN HIGH ORDER ELEMENTS
          IF ( MAT(I,J) .NE. IJFLAG ) THEN                 
            NODE = NODE + 1
            NPTNUM(I,J) = NODE
C            OUTPUT POINTS
           CALL POUT (NODE,R(I,J),Z(I,J),CODE(I,J),
     1                I,J,OUTPT,IOPT)
          ENDIF
   60   CONTINUE
   70 CONTINUE
      WRITE (NPRT,*) ' MAX NODE NUMBER =', NODE
      WRITE (NPIO,*) ' MAX NODE NUMBER =', NODE
      RETURN
      END
      SUBROUTINE  PLOTIJ (MC,MR,R,Z,NPTNUM,MAT,IMAX,IMIN,
     1                    LIMIN,LIMAX,LINC,LTYPE,IPLOT)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     PLOT A GRID OF POINTS DEFINED BY I,J COORDINATES
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  NPTNUM(MC,MR), MAT(MC,MR), IMIN(MR), IMAX(MR),
     1           LIMIN(MR), LIMAX(MR)
      DIMENSION  R(MC,MR), Z(MC,MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      DATA KOUNT /0/
C        FIND LIMITS OF COORDINATES ALONE EDGES OF MESH
      CALL  RZSIZE (MC,MR,R,Z,IMIN,IMAX,NPTNUM,RMIN,RMAX,ZMIN,ZMAX)
C        ADD A BORDER TO REGION
      RDIFF = RMAX - RMIN
      ZDIFF = ZMAX - ZMIN
      BORDER = AMAX1(RDIFF,ZDIFF)/100.
      BORDER = ABS(BORDER)
      RMIN = RMIN - BORDER
      ZMIN = ZMIN - BORDER
      RMAX = RMAX + BORDER
      ZMAX = ZMAX + BORDER
C        ADVANCE FRAME AND SET WINDOW LIMITS
      CALL WINDO (RMIN,RMAX,ZMIN,ZMAX)
C        PLOT THE ROW LINES
      DO 10  J = 1,MR,LINC
      I1 = IMIN(J)
      I2 = IMAX(J) - 1
      JLL = J - LINC
      DO 20  I = I1,I2
        IF ( NPTNUM(I,J).EQ.0 .OR. NPTNUM(I+1,J).EQ.0 )  GO TO 20
        ILL = I - LINC
C          POSSIBLE LINE, MOVE TO POINT, PEN UP
        CALL MOVETO (R(I,J),Z(I,J))
C          IS AN ELEMENT ON EITHER SIDE
        IF ( MAT(I,J).GE.0 )  GO TO 5
        IF ( J.GT.1 .AND. MAT(I,JLL).LT.0 )  GO TO 20
C          LINE EXISTS, MOVE TO NEXT POINT, PEN DOWN
    5   CALL LINETO (R(I+1,J),Z(I+1,J))
   20   CONTINUE
   10 CONTINUE
C        PLOT THE COLUMN LINES
      JLESS = MR - 1
      DO 40  I = 1,MC,LINC
      ILL = I - LINC
      DO 50  J = 1,JLESS
        IF ( NPTNUM(I,J).EQ.0 .OR. NPTNUM(I,J+1).EQ.0 )  GO TO 50
C          POSSIBLE LINE, MOVE TO POINT, PEN UP
        CALL MOVETO (R(I,J),Z(I,J))
C          IS AN ELEMENT ON EITHER SIDE
        IF ( MAT(I,J).GE.0 )  GO TO 35
        IF ( I.GT.1 .AND. MAT(ILL,J).LT.0 )  GO TO 50
C          LINE EXISTS, MOVE TO NEXT POINT, PEN DOWN
   35   CALL LINETO (R(I,J+1),Z(I,J+1))
   50   CONTINUE
   40 CONTINUE
      IF ( IPLOT.EQ.2 )  RETURN
      KOUNT = KOUNT + 1
      IF ( IPLOT.EQ.3 .AND. KOUNT.EQ.2 )  RETURN
C        CHECK FOR TRIANGLES
      IF ( LTYPE.EQ.4  .OR. LTYPE.EQ.8  .OR. LTYPE.EQ.9 .OR.
     1     LTYPE.EQ.12 .OR. LTYPE.EQ.16 .OR. 
     2     LTYPE.EQ.17 .OR. LTYPE.EQ.25 )  RETURN
C        LOOP OVER QUADRILATERALS AND ADD DIAGONAL LINES
      JLESS = MR - LINC
      DO 60  J = 1,JLESS,LINC
      I1 = LIMIN(J)
      I2 = LIMAX(J) - LINC
      JL = J + LINC
      DO 80  I = I1,I2,LINC
      IF ( MAT(I,J).LT.0 )  GO TO 80
      IL = I + LINC
C        ELEMENT EXISTS, FIND SHORTEST DIAGONAL
      D13SQ = (R(I,J)-R(IL,JL))**2+(Z(I,J)-Z(IL,JL))**2
      D24SQ = (R(I,JL)-R(IL,J))**2+(Z(I,JL)-Z(IL,J))**2
      D24SQ = D24SQ*1.2
      IF ( D13SQ .GT. D24SQ )  GO TO 65
C          USE 1-3 DIAGONAL
C        MOVE TO CORNER 1, PEN UP
      CALL MOVETO (R(I,J),Z(I,J))
C        MOVE UP DIAGONAL, PEN DOWN
      CALL LINETO (R(I+1,J+1),Z(I+1,J+1))
      IF ( LINC.GT.1 ) CALL LINETO (R(I+2,J+2),Z(I+2,J+2))
      IF ( LINC.GT.2 ) CALL LINETO (R(I+3,J+3),Z(I+3,J+3))
      IF ( LINC.GE.3 ) CALL LINETO (R(I+4,J+4),Z(I+4,J+4))
      GO TO 80
C          USE SHORTER 2-4 DIAGONAL
C        MOVE TO CORNER 2, PEN UP
   65 CALL MOVETO (R(IL,J),Z(IL,J))
C        MOVE UP DIAGONAL, PEN DOWN
      CALL LINETO (R(IL-1,J+1),Z(IL-1,J+1))
      IF ( LINC.GT.1 ) CALL LINETO (R(IL-2,J+2),Z(IL-2,J+2))
      IF ( LINC.GT.2 ) CALL LINETO (R(IL-3,J+3),Z(IL-3,J+3))
      IF ( LINC.GE.3 ) CALL LINETO (R(IL-4,J+4),Z(IL-4,J+4))
   80 CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE  POUT (NODE,R,Z,CODE,I,J,OUTPT,IOPT)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     OUTPUT NODES, BC CODE, COORDINATES
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER OUTPT
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
      DATA KALL / 0 /
      IF ( KALL .EQ. 0 )  THEN
C        PRINT HEADING ON FIRST CALL
        KALL = 1
        IF ( IOPT .EQ. 2 )  THEN
C          PAFEC TITLE
          WRITE(NPUN,502)
          WRITE(NPPT,502)
  502     FORMAT ('nodes', /,'node   x   y')
        ENDIF
        IF ( IOPT .EQ. 4 )  THEN
C          CAEDS NODE CONTROL UNIVERSAL FILE
          WRITE(NPRT,*) 'WARNING, CAEDS NOT COMPLETE'
          WRITE(NPUN,504)
          WRITE(NPPT,504)
  504     FORMAT ('nodes', /,'node   x   y')
        ENDIF
        IF ( OUTPT .LT. 2 )  WRITE (NPRT,1000)
 1000     FORMAT (/, '  N O D A L   P O I N T   O U T P U T',/,
     1    ' POINT   BC CODE         X-COORD',
     2    '            Y-COORD    (I,J)')
      ENDIF    
C                FOR ALL CALLS
      IBC = CODE
C        PRINT NODAL DATA
      IF ( OUTPT .LT. 2 )  WRITE (NPRT,1010)  NODE, IBC, R, Z, I, J
 1010      FORMAT ( I6, I9, 2(1PE18.7), I6, ',', I3 )
      IF ( OUTPT .EQ. 0 .OR. OUTPT .EQ. 3 )  RETURN
C        PUNCH SPECIAL FORMATTED FILES
      IF ( IOPT .EQ. 0 )  THEN     
C-->    MODEL CODE FORMATS
        WRITE (NPUN,1020)  NODE, IBC, R, Z
        WRITE (NPPT,1020)  NODE, IBC, R, Z
 1020   FORMAT ( 2(I5), 2(F10.4) )
        RETURN
      ENDIF
      IF ( IOPT .EQ. 1 )  THEN     
C-->    FINITE CODE FORMATS
        WRITE (NPUN,1030)  NODE, IBC, R, Z
        WRITE (NPPT,1030)  NODE, IBC, R, Z
 1030   FORMAT ( 2(I10), 2(F10.4) )
        RETURN
      ENDIF
      IF ( IOPT .EQ. 2 )  THEN     
C-->    PAFEC FORMATS
        WRITE (NPUN,1040)  NODE, R, Z
        WRITE (NPPT,1040)  NODE, R, Z
 1040   FORMAT ( I10, 2(F10.4) )
        RETURN
      ENDIF
      IF ( IOPT .EQ. 3 )  THEN     
C-->    FLOW90 FORMATS
        WRITE (NPUN,1050)  NODE, R, Z
        WRITE (NPPT,1050)  NODE, R, Z
 1050   FORMAT ( I10, 2(F10.4) )
        RETURN
      ENDIF
      IF ( IOPT .EQ. 4 )  THEN     
C-->    CAEDS FORMATS
        WRITE (NPUN,1060)  NODE, R, Z
        WRITE (NPPT,1060)  NODE, R, Z
 1060   FORMAT ( I10, 2(F10.4) )
        RETURN
      ENDIF
C        ANOTHER CODE
      RETURN
      END
      SUBROUTINE PTCODE (JPT,NG,KODE,KODES)
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
C     EXTRACT B.C. INDICATORS AT NODE NUMBER JPT
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
      DIMENSION  KODES(NG)
C     JPT=NODE NO.,  NG=NO. PARAMETERS PER NODE
C     KODE=(NG) DIGIT INTEGER CONTAINING BC INDICATORS
C     KODES=VECTOR CONTAINING NG INTEGER CODES (0 OR I)
C        0  IMPLIES NO B. C.
C        I  IMPLIES A B. C.  OF TYPE I
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      NGPLUS = NG + 1
      IOLD = KODE
      ISUM = 0
      DO 10  I = 1,NG
      II = NGPLUS - I
      INEW = IOLD/10
      IK = IOLD - INEW*10
      ISUM = ISUM + IK*10**(I-1)
      IOLD = INEW
 10   KODES(II) = IK
C        WAS DATA RIGHT JUSTIFIED?
      IF ( KODE.GT.ISUM )  WRITE (NPRT,*) 
     1   ' BC NOT RIGHT JUSTIFIED AT NODE', JPT
      RETURN
      END
      SUBROUTINE  PTER (L1,L2,L3,L4,L5,L6,L7,L8,LTOTAL,
     1                  M1,M2,M3,M4,MTOTAL,MAXROW,MAXCOL)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     SET UP ARRAY STORAGE POINTERS
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     FIXED PT ARRAY STORAGE LI: 1-INNER 2-MAT 3-IMAX
C            4-IMIN 5-LIMIN 6-LIMAX 7-IGIVEN 8-NPTNUM
C     FLOATING PT ARRAY STORAGE MI: 1-R 2-Z 3-CODE 4-ALINE
      IRC = MAXROW*MAXCOL
C        INTEGERS: INNER, MAT, IMAX, IMIN, LIMIN, LIMAX,
C                  IGIVEN, NPTNUM
      L1 = 1
      L2 = L1 + IRC
      L3 = L2 + IRC
      L4 = L3 + MAXROW
      L5 = L4 + MAXROW
      L6 = L5 + MAXROW
      L7 = L6 + MAXROW
      L8 = L7 + IRC
      LTOTAL = L8 + IRC
C        REALS: R, Z, CODE, ALINE
      M1 = 1
      M2 = M1 + IRC
      M3 = M2 + IRC
      M4 = M3 + IRC
      MTOTAL = M4 + MAXCOL
      RETURN
      END
      SUBROUTINE  RPRINT (A,NR,NC,IOPT)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
C              PRINTING OF REAL MATRIX A(NR,NC)
C     * * * * * * * * * * * * * * * * * * * * * * * * * *
      DIMENSION  A(1), NCOL(8)
C     IF IOPT=0 USE F FORMAT, OTHERWISE USE E FORMAT
      DO 50  J = 1,NC,8
      JL1 = J - 1
      MAXCOL = 1
      K = NC - JL1
      MAXCOL = MIN0 (K,8)
      MXCLL1 = MAXCOL - 1
      DO 10  L = 1,MAXCOL
   10 NCOL(L) = L + JL1
      WRITE (NPRT,5000)  ( NCOL(N),N=1,MAXCOL )
 5000 FORMAT ('   ROW/COL',I10,7I14)
      DO 40  N = 1,NR
      NL = N + (J-1)*NR
      NH = NL + MXCLL1*NR
      IF ( IOPT )  30,20,30
   20 WRITE (NPRT,5010)  N,( A(I),I=NL,NH,NR )
 5010 FORMAT (I6,8F14.4)
      GO TO 40
   30 WRITE (NPRT,5020)  N,( A(I),I=NL,NH,NR )
 5020 FORMAT (I6,8(1PE14.5))
   40 CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE  RZSIZE (MC,MR,R,Z,IMIN,IMAX,NPTNUM,RMIN,RMAX,
     1                    ZMIN,ZMAX)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     FIND LIMITS OF COORDINATES ALONE EDGES OF MESH
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(MC,MR), Z(MC,MR), NPTNUM(MC,MR), IMIN(MR), IMAX(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      I1 = IMIN(1)
      RMIN = R(I1,1)
      ZMIN = R(I1,1)
      RMAX = RMIN
      ZMAX = ZMIN
C        LEFT AND RIGHT I-J LIMITS
      DO 5  J = 1,MR
        I1 = IMIN(J)
        I2 = IMAX(J)
        RMIN = AMIN1 (RMIN, R(I1,J), R(I2,J) )
        ZMIN = AMIN1 (ZMIN, Z(I1,J), Z(I2,J) )
        RMAX = AMAX1 (RMAX, R(I1,J), R(I2,J) )
    5   ZMAX = AMAX1 (ZMAX, Z(I1,J), Z(I2,J) )
C        BOTTOM AND TOP I-J ROWS
      I1 = IMIN(1)
      I2 = IMAX(1)
      DO 10  I = I1,I2
        IF ( NPTNUM(I,1).EQ.0 )  GO TO 10
        RMIN = AMIN1 (RMIN, R(I,1) )
        ZMIN = AMIN1 (ZMIN, Z(I,1) )
        RMAX = AMAX1 (RMAX, R(I,1) )
        ZMAX = AMAX1 (ZMAX, Z(I,1) )
   10 CONTINUE
      I1 = IMIN(MR)
      I2 = IMAX(MR)
      DO 20  I = I1,I2
        IF ( NPTNUM(I,MR).EQ.0 )  GO TO 20
        RMIN = AMIN1 (RMIN, R(I,MR) )
        ZMIN = AMIN1 (ZMIN, Z(I,MR) )
        RMAX = AMAX1 (RMAX, R(I,MR) )
        ZMAX = AMAX1 (ZMAX, Z(I,MR) )
   20 CONTINUE
      RETURN
      END
      SUBROUTINE  STARTP
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INSTALLATION DEPENDENT PLOT INITIALIZATION
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      PARAMETER ( IPORT = 0 )
      EXTERNAL PLOTS
C     CALCOMP CONTROL
      WRITE (NPRT,*)  ' BEGIN PLOT'
C     IPORT = 0 Landscape, = 1 Portrait
      CALL PLOTS (IPORT, 0, 0, 'CALCOMP')
c     CALL PLOTS (IPORT, 0, 0, '/u4/akin/camd/mesh/m90/CALCOMP')
      RETURN
      END
      SUBROUTINE SWAP (I,J)
C     SWAP VALUES OF I,J
      K = J
      J = I
      I = K
      RETURN
      END
      SUBROUTINE  TOOT
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C                       TOOT YOUR OWN HORN
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      WRITE (NPRT,5000)
 5000 FORMAT ('                              MESH90',/,
     A'              A TWO DIMENSIONAL MESH GENERATOR',//,
     B'       PROGRAMMED BY:',//,
     C'                     PROF. J. E. AKIN',/,
     D'        RICE UNIVERSITY',/,
     E'        DEPT. OF MECHANICAL ENGINEERING AND',/,
     F'              MATERIALS SCIENCE',/,
     G'        HOUSTON, TX  77251-1892')
      RETURN
      END
      SUBROUTINE  TWOPT (MC,MR,R,Z,INNER,IMIN,IMAX,N1)
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     GENERATE STRAIGHT LINES BETWEEN GIVEN POINTS ON J LINE
C     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      DIMENSION  R(MC,MR), Z(MC,MR), INNER(MC,MR), 
     1           IMIN(MR), IMAX(MR)
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
      J2 = MR - 1
      DO 50  J = 2,J2
      I1 = IMIN(J) + 1
      I2 = IMAX(J)
      IBACK = I1 - 1
      DO 40  I = I1,I2
      IF ( INNER(I,J) .EQ. 0 )   GO TO 30
      DELTAI = I - IBACK
      DR = ( R(I,J)-R(IBACK,J) )/DELTAI
      DZ = ( Z(I,J)-Z(IBACK,J) )/DELTAI
      ISTR1 = IBACK + 1
      ISTR2 = I - 1
      IF ( ISTR2 .LT. ISTR1 )  GO TO 30
        DO 10  ISTR = ISTR1,ISTR2
        R(ISTR,J) = R(ISTR-1,J) + DR
   10   Z(ISTR,J) = Z(ISTR-1,J) + DZ
        IBACK = I
   30 CONTINUE
   40 CONTINUE
   50 CONTINUE
      N1 = 1
      RETURN
      END
      SUBROUTINE  VALUE  (MC,MR,ALINE,INNER,IMIN,IMAX,
     1                    NPTNUM,ECHO,NBCCD,IJ2JI)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     I/O BC VALUES FOR LINES IN I-J AND R-Z SPACES
C     I1---I3---I2   ===TO==>>   F1---F3---F2
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      INTEGER ECHO 
      DIMENSION   INNER(MR,MC), NPTNUM(MR,MC),
     1            IMIN(MR),  IMAX(MR)
      DIMENSION   ALINE(MC) 
      COMMON / IOUNIT / NCRD, NPRT, NPUN, NPLT
      COMMON / PROB / NPIO, NPPT, NPEL, NPBC, NPBV, NPLM, NPIP, NPOP
C     INNER(I,J)      > 0 IMPLIES PT (I,J) WAS INPUT
C     IMIN(J),IMAX(J) = MIN, MAX I FOR ANY PT IN ROW J
C       ZERO INNER FOR NEW USE
      DO 2  J = 1,MC
        DO 1  I = 1,MR
    1   INNER(I,J) = 0
    2 CONTINUE
      KOUNT = 0
      IF ( ECHO .EQ. 0 )  WRITE (NPRT,5000)
 5000   FORMAT ( //,' NODAL ESSENTIAL BC INPUT DATA',/,
     1  '    NODE NG    VALUE        (I,  J)')
C       INPUT BC VALUES
      DO 4  N = 1,NBCCD
C==>>  MEANING IF INTEST:
C      =0, SINGLE I,J AND PREVIOUS R,Z POINT
C      =3 or 1, LINE IN I,J TO CONNECT 2 PREVIOUS R,Z ENDS
C      =4 or 2, LINE OR ARC IN I,J TO CONNECT 3 PREVIOUS R,Z PTS
C     NTEST       REQUIRED ?  (Y-YES, N-NO, O-OPTIONAL)
C     VALUE    I1,J1  I2,J2  I3,J3   F1     F2     F3  
C     0          Y      N      N      Y      N      N
C     3 or 1     Y      Y      N      Y      O      N
C     4 or 2     Y      Y      O      Y      Y      O
C     NG  = COMPONENT NUMBER BEING SPECIFIED
        READ (NCRD,5010) INTEST, NG, I1, J1, I2, J2, I3, J3
 5010   FORMAT ( 8I4 )
C       F1,F2,F3 = ENDS & MIDSIDE VALUES BEING INPUT             
        READ (NCRD,5011) F1, F2, F3
 5011   FORMAT ( 7F10.0 )
        IF ( NG .LT. 1 )  NG = 1
        IF ( IJ2JI .EQ. 0 )  THEN    
          CALL SWAP(I1,J1)
          CALL SWAP(I2,J2)
          CALL SWAP(I3,J3)
        ENDIF
C        SINGLE POINT VALUE
        IF ( INTEST .EQ. 0 )  THEN   
          I2 = I1
          J2 = J1
          I3 = I1
          J3 = J1
          F2 = F1
          F3 = F1
        ENDIF
C        LINEAR VALUES BETWEEN 2 POINTS
        IF ( INTEST .EQ. 3 .OR. INTEST .EQ. 1 ) 
     1     F3 = (F2 + F1)*0.5
C       PARABOLIC VALUES BETWEEN 3 POINTS
        IF ( I2 .EQ. 0 )  I2 = I1
        IF ( J2 .EQ. 0 )  J2 = J1
C        FIND BIGGEST LOOP REQUIRED
        LOOP = MAX0 ( IABS(I2-I1), IABS(J2-J1) ) + 1
        DL   = LOOP - 1
        IF ( LOOP .EQ. 1 )  DL = 1.0
C        LOOP OVER POSSIBLE I,J PAIRS
        A1 = I1
        A2 = I2
        A3 = I3
        IF ( I3 .EQ. 0 )  A3 = FLOAT(I1+I2)/2.
        B1 = J1
        B2 = J2
        B3 = J3
        IF ( J3 .EQ. 0 )  B3 = FLOAT(J1+J2)/2.
        DO 3  L = 1, LOOP
C          COMPUTE UNIT COORDINATE
          X = FLOAT(L-1)/DL
C          LOCAL QUADRATIC INTERPOLATIONS
          H1 = (X-1.)*(X+X-1.)
          H2 = X*(X+X-1.)
          H3 = 4.*X*(1.-X)
C          INTERPOLATE IN I,J FOR VALUES
          I = IFIX( A1*H1 + A2*H2 + A3*H3 + 0.1 )
          J = IFIX( B1*H1 + B2*H2 + B3*H3 + 0.1 )
          IF ( I .LT. 1 )   I = 1
          IF ( I .GT. MC )  I = MC
          IF ( J .LT. 1 )   J = 1
          IF ( J .GT. MR )  J = MR
          NODE = NPTNUM(I,J)
C          PARABOLIC ARC
          F = F1*H1 + F2*H2 + F3*H3
          IF ( NODE .GE. 1 )  THEN
            IF ( ECHO .EQ. 0 ) WRITE (NPRT,5080) NODE,NG,F,I,J
 5080       FORMAT ( I8, I3, 1PE13.5, I6, I5 )
            WRITE (NPBV,5080) NODE,NG,F,I,J
C            FLAG VALUES FOR PRINT-PLOT
            INNER(I,J) = 1
            KOUNT      = KOUNT + 1
          ELSE
            WRITE (NPRT,*) ' SKIPPED I,J', I, J
          ENDIF
 3      CONTINUE
        IF ( ECHO .EQ. 0 )  WRITE (NPRT,*) 
     1    '----- END OF BC GENERATION CARD ', N
 4    CONTINUE
      WRITE (NPRT,*) ' NUMBER OF BC VALUES CREATED =', KOUNT
      WRITE (NPIO,*) ' NUMBER OF BC VALUES CREATED =', KOUNT
C      PRINT_PLOT BC NODES
      ISTEP = 4
      CALL  IJMESH (MC,MR,ALINE,INNER,IMAX,IMIN,ISTEP)
      RETURN
      END
      SUBROUTINE WINDO (RMIN,RMAX,ZMIN,ZMAX)
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     INSTALLATION DEPENDENT PLOTTER LIMITS
C    - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      COMMON /IOUNIT/ NCRD, NPRT, NPUN, NPLT
C     CALCOMP CONTROL
      COMMON /CALCOM/ FIRSTR, FIRSTZ, DELTA
      DIMENSION RFAKE(4), ZFAKE(4)
      DATA RSCREN, ZSCREN /10., 7.8/
      EXTERNAL SCALE
C     SCALE DATA
      RFAKE(1) = RMIN
      RFAKE(2) = RMAX
      ZFAKE(1) = ZMIN
      ZFAKE(2) = ZMAX
      CALL SCALE (RFAKE,RSCREN,2,1)
      CALL SCALE (ZFAKE,ZSCREN,2,1)
      FIRSTR = RFAKE(3)
      FIRSTZ = ZFAKE(3)
      DELTA = AMAX1(RFAKE(4),ZFAKE(4))
      write(nprt,*) rfake, zfake,firstr, firstz, delta
      IF ( DELTA .LE. 0.0 )  THEN
        WRITE (NPRT,*) 'BAD DELTA IN WINDO'
        DELTA = 1.0
        FIRSTR = RMIN
        FIRSTZ = ZMIN
      ENDIF
      RETURN
      END
