C     PROGRAM Forearc
C
C     A program to help in preparing thin-plate finite element
C     grids (.feg files) which include a forearc region in 
C     the upper plate of a subduction zone.
C
C     Before running this program, a grid is prepared and
C     filled with elevation, heat-flow, crustal-thickness,
C     and mantle-lithosphere thickness data at nodes.
C    (This is done by assuming local isostasy and an isothermal
C     surface at the base of the lithosphere.)
C
C     Then, this program Forearc is used to replace some
C     values of crustal thickness and mantle lithosphere
C     thickness in the forearc region.
C
C     In addition to the .feg file with nodal data, this program
C     also requires:
C       (1) A grid (.grd) file of the depth to the subduction
C           shear zone, measured from sea level.
C       (2) A grid (.grd) file of the depth to the Moho,
C           measured from sea level.
C     Grid (1) must cover the whole area of the model, but
C     only needs to be accurate in the forearc area,
C     which will be defined as the area where the subduction
C     shear zone is shallower than some limit depth.
C     Grid (2) is only needed in the forearc region.
C
C     The output from this program is a modified .feg file,
C     which will NOT be isostatic or isothermally-limited in
C     the forearc region.
C
C     by Peter Bird, Dept. of Earth & Space Sciences, 
C     University of California, Los Angeles, CA 90095-1567.
C    (For version date, see 1st WRITE below.)
C
C---------------------------------------------------------------------
C                 PARAMETER (ARRAY-SIZE) STATEMENTS
C
C   SET THE FOLLOWING PARAMETERS AT LEAST AS LARGE AS YOUR PROBLEM:
C
C MAXNOD = MAXIMUM NUMBER OF NODES (INCLUDES BOTH "REAL"AND & "FAKE")
       INTEGER, PARAMETER :: MAXNOD = 5000
C
C MAXEL = MAXIMUM NUMBER OF CONTINUUM ELEMENTS (TRIANGLES).
       INTEGER, PARAMETER :: MAXEL = 10000
C
C MAXFEL = MAXIMUM NUMBER OF FAULT ELEMENTS (LINE SEGMENTS);
       INTEGER, PARAMETER :: MAXFEL = 2000
C
C MAXDAT = MAXIMUM NUMBER OF ROWS AND COLUMNS IN 2-D ARRAYS OF
C          ELEVATION, HEAT-FLOW DATA, AND CRUSTAL-THICKNESS DATA
       INTEGER, PARAMETER :: MAXDAT = 900
C
C---------------------------------------------------------------------
C                            TYPE STATEMENTS
C
C  (NOTE: THE IMPLICIT TYPING OF I-N = INTEGER, AND A-H,O-Z = REAL
C   IS OBSERVED IN THIS PROGRAM.  NO TYPES ARE STATED FOR SUCH NAMES.)
C
       CHARACTER*80 TITLE1
C     
       INTEGER CHOICE,I,IROW,IUNITT,
     +         JCOL,LIST,
     +         MXDATA,MXEL,MXFEL,MXNODE,
     +         NMX,NMY,NSX,NSY,
     +           NEMAX,NFAKEN,NFL,NODEF,NODES,
     +           NREALN,NUMEL,NUMNOD,N1000
C
       LOGICAL BRIEF,CHECKE,CHECKF,CHECKN,
     +         INPUT,MODIFY,
     +         more_feg,more_grd,
     +         OUTPUT,PLATES,SHELLS,WORK
C
       REAL ARCLIM,
     +      CRUST,
     +      DIPMAX,DMOHO,DQDTDA,DSHEAR,
     +      ELEV,
     +      FAZ,FDIP,
     +      MANTLE,MARRAY,MDX,MDY,MX1,MX2,MY1,MY2,
     +      OFFMAX,OFFSET,
     +      PLATE,
     +      QLIMIT,
     +      SARRAY,SDX,SDY,SX1,SX2,SY1,SY2,
     +      TLNODE,TOPOG,
     +      X,XNODE,
     +      Y,YNODE,
     +      ZMNODE
C
C---------------------------------------------------------------------
C                        DIMENSION STATMENTS
C
C  DIMENSIONS USING PARAMETER MAXNOD:
       DIMENSION CHECKN(MAXNOD), DQDTDA(MAXNOD),
     +           ELEV  (MAXNOD), 
     +           TLNODE(MAXNOD),
     +           XNODE (MAXNOD), YNODE (MAXNOD), ZMNODE(MAXNOD)
C
C  DIMENSIONS USING PARAMETER MAXEL:
       DIMENSION CHECKE(MAXEL), NODES(6,MAXEL)
C
C  DIMENSIONS USING PARAMETER MAXFEL:
       DIMENSION CHECKF(MAXFEL), FAZ (2,MAXFEL), FDIP  (3,MAXFEL),
     +           NODEF(6,MAXFEL),
     +           OFFSET (MAXFEL)
C
C  DIMENSIONS USING PARAMETER MAXDAT:
       DIMENSION MARRAY (MAXDAT,MAXDAT),
     +           SARRAY (MAXDAT,MAXDAT)
C---------------------------------------------------------------------
C
C                           DATA STATEMENTS
C
C   "IUNITT" = FORTRAN DEVICE NUMBER FOR TEXT MESSAGES:
       DATA IUNITT /6/
C
C---------------------------------------------------------------------
C
C                   BEGINNING OF EXECUTABLE CODE
C
C   *** KLUDGE ALERT ***
C   CONVERSION OF PARAMETERS (CONSTANTS) TO VARIABLES SHOULD LOGICALLY
C   HAVE NO EFFECT, BUT IN FACT HELPS TO SUPPRESS SOME SPURIOUS
C   MESSAGES FROM THE IBM VS-FORTRAN COMPILER.
       MXNODE=MAXNOD
       MXEL  =MAXEL
       MXFEL =MAXFEL
       MXDATA=MAXDAT
C
C---------------------------------------------------------------------
C
C  INTRODUCTION / HEADER SECTION:
C
    1  WRITE (IUNITT,10)
   10  FORMAT (
     + /'  PROGRAM Forearc  (version of 10 April 2000)'
     +//'  A program to help in preparing thin-plate finite element'
     + /'  grids (.feg files) which include a forearc region in '
     + /'  the upper plate of a subduction zone.'
     +//'  Are you currently using:'
     + /'    (1) DrawGrid to create a flat-Earth .feg in (x,y)'
     + /'        coordinates, for use with Plates, or'
     + /'    (2) OrbWeave to create a round-Earth .feg in (lon,lat)'
     + /'        coordinates, for use with Shells?'
     + /'  Please enter the integer that describes your choice: '\)
       READ (*,*) CHOICE
       IF ((CHOICE.LT.1).OR.(CHOICE.GT.2)) THEN
            WRITE (*,"(' ERROR: Illegal selection.')")
            CALL PAUSE()
            GO TO 1
       END IF
       IF (CHOICE.EQ.1) THEN
            PLATES=.TRUE.
            SHELLS=.FALSE.
       ELSE IF (CHOICE.EQ.2) THEN
            PLATES=.FALSE.
            SHELLS=.TRUE.
       END IF
       IF (PLATES) THEN
            WRITE (IUNITT,21)
   21       FORMAT (
     + /'  Before running this program, a grid is prepared '
     + /'  using the interactive editor DrawGrid, and is'
     + /'  filled with elevation, heat-flow, crustal-thickness,'
     + /'  and mantle-lithosphere thickness data at nodes,'
     + /'  using program FillGrid.'
     + /' (This is done by assuming local isostasy and an isothermal'
     + /'  surface at the base of the lithosphere.)'
     +//'  Then, this program Forearc is used to replace some'
     + /'  values of crustal thickness and mantle lithosphere'
     + /'  thickness in the forearc region.')
       ELSE IF (SHELLS) THEN
            WRITE (IUNITT,22)
   22       FORMAT (
     + /'  Before running this program, a grid is prepared'
     + /'  using the interactive editor OrbWeaver, and is'
     + /'  filled with elevation, heat-flow, crustal-thickness,'
     + /'  and mantle-lithosphere thickness data at nodes,'
     + /'  using program OrbData.'
     + /' (This is done by assuming local isostasy and an isothermal'
     + /'  surface at the base of the lithosphere.)'
     +//'  Then, this program Forearc is used to replace some'
     + /'  values of crustal thickness and mantle lithosphere'
     + /'  thickness in the forearc region.')
       END IF
       CALL PAUSE()
       WRITE (IUNITT,30)
   30  FORMAT (
     + /'  In addition to the .feg file with nodal data, this program'
     + /'  also requires:'
     +//'    (1) A grid (.grd) file of the depth to the subduction'
     + /'         shear zone (thrust fault), measured from sea level.'
     + /'    (2) A grid (.grd) file of the depth to the Moho,'
     + /'         measured from sea level.'
     +//'  Grid (1) must cover the whole area of the model, but'
     + /'  only needs to be accurate in the forearc area,'
     + /'  which will be defined as the area where the subduction'
     + /'  shear zone is shallower than some limit depth.'
     + /'  Grid (2) only needs to cover the forearc region;'
     + /'  if values are given outside the forearc, they are not used.'
     + /'  It may happen that the Moho in dataset (2) is deeper than'
     + /'  the shear zone depth in dataset (1), perhaps because the'
     + /'  Moho was found in the subducting plate.  In this case, '
     + /'  the shallower depth from dataset (1) will be used.')
       IF (PLATES) THEN
            WRITE (IUNITT,31)
   31       FORMAT (
     + /' (The regular grid should be defined in (x,y) Cartesian'
     + /'  coordinates, using the same units as you used for the'
     + /'  node locations in your .feg file.)')
       ELSE IF (SHELLS) THEN
            WRITE (IUNITT,32)
   32       FORMAT (
     + /' (The regular grid should be defined in decimal degrees'
     + /'  of East longitude (- for West) and North latitude'
     + /' (- for South), just like the node locations in your'
     + /'  .feg file.)')
       END IF
       CALL PAUSE()
       WRITE (IUNITT,40)
   40  FORMAT (
     + /'  The output from this program is a modified .feg file,'
     + /'  which will NOT be isostatic or isothermally-limited in'
     + /'  the forearc region.')
       CALL PAUSE()
       WRITE (IUNITT,*)
       CALL Prompt_for_Logical(
     +     "Do you want information about .feg files?",.TRUE.,more_feg)
       IF (more_feg) THEN
            IF (PLATES) THEN
                 WRITE (IUNITT,51)
   51            FORMAT (
     +  ' ------------------------------------------------------------'
     + ,'----------'
     + /'                About Finite Element Grid (.feg) Files'
     +//' Flat-Earth finite-element grids for use with F-E program'
     + ,' Plates are'
     + /' created in program DrawGrid, renumbered by program Number'
     + ,' and filled'
     + /' with nodal data by program FillGrid.'
     +//' They are composed of 6-node isoparametric-triangle continuum'
     + ,' elements and'
     + /' 6-node curvilinear fault elements with specified dip'
     + ,' angles.'
     +//' At each node point, nodal data include Cartesian x, y,'
     + /' elevation/depth of the surface with respect to sea level,'
     + ,' heat-flow,'
     + /' crustal thickness, and mantle-lithosphere thickness (not'
     + ,' including'
     + /' crust).'
     +//' Dip angles, in degrees, are given at three points in each'
     + ,' fault element.'
     + /' ------------------------------------------------------------'
     + ,'----------')
            ELSE IF (SHELLS) THEN
                 WRITE (IUNITT,52)
   52            FORMAT (
     +  ' ------------------------------------------------------------'
     + ,'----------'
     + /'                About Finite Element Grid (.feg) Files'
     +//' Spherical-Earth finite-element grids for use with F-E program'
     + ,' Shells are'
     + /' created in program OrbWeaver, renumbered by program'
     + ,' OrbNumber, and filled'
     + /' with nodal data by program OrbData.'
     +//' They are composed of 3-node spherical-triangle continuum'
     + ,' elements and'
     + /' 4-node great-circle-arc fault elements with specified dip'
     + ,' angles.'
     +//' At each node point, nodal data include East longitude,'
     + ,' North latitude,'
     + /' elevation/depth of the surface with respect to sea level,'
     + ,' heat-flow,'
     + /' crustal thickness, and mantle-lithosphere thickness'
     + ,' (not including'
     + /' crust).'
     +//' Dip angles, in degrees, are given for each end of each'
     + ,' fault element.'
     + /' ------------------------------------------------------------'
     + ,'----------')
            END IF
            CALL PAUSE()
       END IF
       CALL Prompt_for_Logical(
     +    "Do you want information about .grd files?",.TRUE.,more_grd)
       IF (more_grd) THEN
            IF (PLATES) THEN
                 WRITE (IUNITT,61)
   61            FORMAT (
     + /' -------------------------------------------------------------'
     +,'---------'
     + /'                    About Gridded Data (.grd) Files'
     +//' These files contain scalar data values on a regular rectangul'
     +,'ar grid,'
     + /'    defined in Cartesian (x, y) space.'
     +//' For flat-planet (Plates) problems:'
     + /' The first line has 3 numbers: X_min, d_X, X_max;'
     + /'   the 2nd also has 3 numbers: Y_min, d_Y, Y_max.')
            ELSE IF (SHELLS) THEN
                 WRITE (IUNITT,62)
   62            FORMAT (
     + /' -------------------------------------------------------------'
     +,'---------'
     + /'                    About Gridded Data (.grd) Files'
     +//' These files contain scalar data values on a regular rectangul'
     +,'ar grid,'
     + /' defined in Cartesian (East longitude, North latitude) space,'
     + /' using units of decimal degrees (+ = E, N; - = W, S).'
     +//' For spherical-planet (Shells) problems:'
     + /' The first line has 3 numbers: Lon_min, d_Lon, Lon_max;'
     + /'   the 2nd also has 3 numbers: Lat_min, d_Lat, Lat_max.')
            END IF
            WRITE (*,"(
     + /' Following lines give the gridded data in text order, i.e., '
     +,'beginning'
     + /'   with the NorthWest corner, going W->E along the top row, th'
     +,'en W->E '
     + /'   along the 2nd row, etc.  The number and position of line br'
     +,'eaks'
     + /'   is NOT important in this part of the file.'
     +//' -------------------------------------------------------------'
     +,'---------')")
            CALL PAUSE()
       END IF
C
C   READ FINITE ELEMENT GRID ON UNIT 1
C
C   (note: using modified GETNET that reads both types of .feg)
C
       CALL GETNET (INPUT,PLATES,SHELLS,
     +                       1,IUNITT,
     +                    MXEL,MXFEL,MXNODE,
     +             OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP,
     +                    NFAKEN,NFL,NODEF,NODES,NREALN,
     +                    NUMEL,NUMNOD,N1000,OFFMAX,OFFSET,
     +                    TITLE1,TLNODE,XNODE,YNODE,ZMNODE,
     +               WORK,CHECKE,CHECKF,CHECKN)
       WRITE (IUNITT,70)
   70  FORMAT (/' The FE grid has been read.')
C
C   READ IN SUBDUCTION-ZONE-DEPTH-BELOW-SEA-LEVEL GRID (SARRAY):
C
       WRITE(*,81)
   81  FORMAT(/ /' Attempting to read .grd file of'
     +          /' Subduction Zone Depth Below Sea Level'
     +          /' from unit 2...'/)
       READ (2,*) SX1,SDX,SX2
       READ (2,*) SY1,SDY,SY2
       NSX=(SX2-SX1)/SDX + 1.5
       NSY=(SY2-SY1)/SDY + 1.5
       IF ((NSX.GT.MXDATA).OR.(NSY.GT.MXDATA)) THEN
            NEMAX=MAX(NSX,NSY)
            WRITE (IUNITT,85) NEMAX
   85       FORMAT (/' Increase parameter MAXDAT to ',I5,
     +              /' or more and recompile.')
            STOP
       ENDIF
       READ (2,*) ((SARRAY(IROW,JCOL),JCOL=1,NSX),IROW=1,NSY)
       CLOSE(2)
C
C   READ IN MOHO-DEPTH-BELOW-SEA-LEVEL GRID (MARRAY):
C
       WRITE(*,91)
   91  FORMAT(/ /' Attempting to read grd file of'
     +          /' depth to Moho, below sea level,'
     +          /' from unit 3...'/)
       READ (3,*) MX1,MDX,MX2
       READ (3,*) MY1,MDY,MY2
       NMX=(MX2-MX1)/MDX + 1.5
       NMY=(MY2-MY1)/MDY + 1.5
       IF ((NMX.GT.MXDATA).OR.(NMY.GT.MXDATA)) THEN
            NEMAX=MAX(NMX,NMY)
            WRITE (IUNITT,85) NEMAX
            STOP
       ENDIF
       READ (3,*) ((MARRAY(IROW,JCOL),JCOL=1,NMX),IROW=1,NMY)
       CLOSE(3)
C
C   DEFINE THE FOREARC REGION
C
       WRITE(*,100)
  100  FORMAT(/' The forearc region will be defined as the region'
     +        /' in which the depth (below sea level) of the subduction'
     +        /' zone is less than some limit.'
     +        /' Enter this limit, using the same units as the depths'
     +        /' in the .grd file already read: '\)
       READ (*,*) ARCLIM
C
       WRITE(*,200)
  200  FORMAT(/' It is usually desirable to limit the heat-flow of'
     +        /' the forearc region, to prevent unreasonable'
     +        /' temperatures along the subduction shear zone.'
     +        /' Enter the desired heat-flow limit, using the same'
     +        /' units as in the .grd and .feg files already read: '\)
       READ (*,*) QLIMIT
C
C   PROCESS ALL NODES:
C
       WRITE (IUNITT,600)
  600  FORMAT (/' Computing layer thicknesses at forearc nodes...'
     +         /'      0 nodes completed...')
       DO 680 I=1,NUMNOD
C
C         NODE POSITION
C
            X=XNODE(I)
            Y=YNODE(I)
C
C         GET DEPTH (BELOW SEA LEVEL) TO SUBDUCTION ZONE:
C
            CALL INTERP (INPUT,SARRAY,
     +                         MAXDAT,NSX,NSY,
     +                         SX1,SDX,SX2,
     +                         SY1,SDY,SY2,
     +                         X,Y,
     +                  OUTPUT,DSHEAR)
C
            IF (DSHEAR.LE.ARCLIM) THEN
C
C              GET LAYER THICKNESSES BEFORE ANY CHANGE:
C
                 TOPOG=ELEV(I)
                 CRUST=ZMNODE(I)
                 MANTLE=TLNODE(I)
                 PLATE=CRUST+MANTLE
C
C              GET PROPER CRUSTAL THICKNESS IN FOREARC REGION:
C
                 CALL INTERP (INPUT,MARRAY,
     +                              MAXDAT,NMX,NMY,
     +                              MX1,MDX,MX2,
     +                              MY1,MDY,MY2,
     +                              X,Y,
     +                       OUTPUT,DMOHO)
C
C              NOTE THAT THIS MOHO MAY BE IN THE SUBDUCTING PLATE!
C             (ESPECIALLY IF IT IS DEFINED AS THE SEISMIC MOHO)
C              LIMIT THIS DEPTH TO BE NO MORE THAN DSHEAR:
C
                 DMOHO=MIN(DMOHO,DSHEAR)
                 DMOHO=MAX(DMOHO,0.0)
C
C              MODIFY THE LAYER THICKNESSES:
C
                 CRUST=MAX((TOPOG+DMOHO),0.0)
                 MANTLE=MAX((DSHEAR-DMOHO),0.0)
C
C              STORE THE MODIFIED LAYER THICKNESSES:
C
                 ZMNODE(I)=CRUST
                 TLNODE(I)=MANTLE
C
C              LIMIT THE HEAT-FLOW OF THE FOREARC:
C
                 DQDTDA(I)=MIN(DQDTDA(I),QLIMIT)
C
            END IF
            WRITE (IUNITT,679) I
  679       FORMAT ('+',I6,' nodes completed...')
  680  CONTINUE
C
C   OUTPUT THE MODIFIED .FEG FILE:
C
C  (note: using modified PUTNET that writes both types of .feg)
C
       CALL PUTNET (INPUT,PLATES,SHELLS,
     +                        4,
     +                    BRIEF,DQDTDA,ELEV,FAZ,FDIP,
     +                    MXEL,MXFEL,MXNODE,N1000,
     +                    NFAKEN,NFL,NODEF,NODES,
     +                    NREALN,NUMEL,NUMNOD,OFFSET,
     +                    TITLE1,TLNODE,XNODE,YNODE,ZMNODE)
       WRITE (IUNITT,800)
  800  FORMAT (/' Remember to edit the title line of ',
     +           'this new .feg file!')
       CALL PAUSE()
       STOP
       END
C
C
C
       REAL FUNCTION ATAN2F (Y,X)
C    
C    CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS
C
       REAL X,Y
       IF ((X.NE.0.).OR.(Y.NE.0.)) THEN
            ATAN2F=ATAN2(Y,X)
       ELSE
            ATAN2F=0.
       END IF
       RETURN
       END FUNCTION ATAN2F
C
C
C
       SUBROUTINE GETNET (INPUT,PLATES,SHELLS,
     +                          IUNIT7,IUNIT8,
     +                          MXEL,MXFEL,MXNODE,
     +                   OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP,
     +                          NFAKEN,NFL,NODEF,NODES,NREALN,
     +                          NUMEL,NUMNOD,N1000,OFFMAX,OFFSET,
     +                          TITLE1,TLNODE,XNODE,YNODE,ZMNODE,
     +                     WORK,CHECKE,CHECKF,CHECKN)
C
C   Note: Special version of GetNet; works with both kinds of .feg;
C         flat-Earth (PLATES = T) or round-Earth (SHELLS = T).
C         Also, code is kept simple by assuming BRIEF = T.
C   READ FINITE ELEMENT GRID FROM UNIT "IUNIT7".
C   ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8".
C
       CHARACTER*80 TITLE1
       INTEGER I,IFN,INDEX,IOS,IUNIT7,IUNIT8,
     +         J,K,L,
     +         MXEL,MXFEL,MXNODE,
     +         N,NBASE,NFAKEN,NFL,
     +           NODEF,NODES,NREALN,NTOP,NUMEL,NUMNOD,N1000
       LOGICAL ALLOK,BRIEF,CHECKE,CHECKF,CHECKN,INPUT,OUTPUT,
     +         PLATES,SHELLS,WORK
       REAL AZ1,AZ3,
     +      DIPS,DQDTDA,ELEV,ELEVI,FAZ,FDIP,
     +      OFF,QI,OFFMAX,OFFSET,
     +      TLI,TLNODE,VECTOR,
     +      XI,XNODE,YI,YNODE,ZMI,ZMNODE
C
C   EXTERNAL ARRAYS:
       DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE),
     +           DQDTDA(MXNODE),ELEV(MXNODE),
     +           FAZ(2,MXFEL),FDIP(3,MXFEL),
     +           NODEF(6,MXFEL),NODES(6,MXEL),OFFSET(MXFEL),
     +           TLNODE(MXNODE),
     +           XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE)
C   INTERNAL ARRAYS:
       DIMENSION DIPS(3),IFN(6),VECTOR(7)
C
       WRITE (IUNIT8,1) IUNIT7
    1  FORMAT (/' Attempting to read FINITE ELEMENT GRID from unit',I3/)
       TITLE1='                                        '//
     +        '                                        '
       READ (IUNIT7,2,IOSTAT=IOS) TITLE1
    2  FORMAT (A80)
       WRITE (IUNIT8,3) TITLE1
    3  FORMAT(/' TITLE OF FINITE ELEMENT GRID ='/' ',A80)
C
C   READ NUMBER OF NODES, AND HOW MANY ARE "REAL" VERSUS "FAKE".
C   INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC
C   GRAVITY ANOMALIES (ONE RECORD PER NODE).
C    (OPTION "BRIEF" SUPPRESSES MOST OUTPUT)
C
       READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,N1000,BRIEF
C
       IF (NUMNOD.NE.(NREALN+NFAKEN)) THEN
            WRITE (IUNIT8,5)
    5       FORMAT (/' INCONSISTENT DATA:'/
     +               ' NUMBER OF NODES SHOULD EQUAL TOTAL OF REAL',
     +               ' NODES AND FAKE NODES.')
            CALL PAUSE()
            STOP
       ENDIF
       IF (NUMNOD.GT.MXNODE) THEN
            WRITE (IUNIT8,10) NUMNOD
   10       FORMAT(/' INCREASE PARAMETER MAXNOD TO BE AT LEAST'
     +           /' THE NUMBER OF NODES (',I6,') AND RECOMPILE.')
            CALL PAUSE()
            STOP
       ENDIF
       IF (NREALN.GT.N1000) THEN
            WRITE (IUNIT8,20) NREALN
   20       FORMAT (/' INCREASE THE DATA VALUE N1000 TO BE GREATER'
     +              /' OR EQUAL TO NREALN (',I6,') AND RECOMPILE.')
            CALL PAUSE()
            STOP
       ENDIF
C
       NBASE=N1000+1
       NTOP=N1000+NFAKEN
       DO 90 K=1,NUMNOD
            CHECKN(K)=.FALSE.
   90  CONTINUE
       DO 100 K=1,NUMNOD
            CALL READN (INPUT,IUNIT7,IUNIT8,7,
     +                 OUTPUT,VECTOR)
            INDEX=VECTOR(1)+0.5
            XI=VECTOR(2)
            YI=VECTOR(3)
            ELEVI=VECTOR(4)
            QI=VECTOR(5)
            ZMI=VECTOR(6)
            TLI=VECTOR(7)
            IF (INDEX.LE.NREALN) THEN
                 I=INDEX
            ELSE
                 I=INDEX-N1000+NREALN
            ENDIF
            CHECKN(I)=.TRUE.
            ELEV(I)=ELEVI
            DQDTDA(I)=QI
            IF (QI.LT.0.) THEN
                 WRITE (IUNIT8,96)
   96            FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.')
                 CALL PAUSE()
                 STOP
            ENDIF
            XNODE(I)=XI
            YNODE(I)=YI
            IF (ZMI.LT.0.) THEN
                 WRITE (IUNIT8,97)
   97            FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.')
                 CALL PAUSE()
                 STOP
            ENDIF
            ZMNODE(I)=ZMI
            IF (TLI.LT.0.) THEN
                 WRITE (IUNIT8,98)
   98            FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS IS',
     +                  ' NON-PHYSICAL.')
                 CALL PAUSE()
                 STOP
            ENDIF
            TLNODE(I)=TLI
  100  CONTINUE
       ALLOK=.TRUE.
       DO 101 I=1,NUMNOD
            ALLOK=ALLOK.AND.CHECKN(I)
  101  CONTINUE
       IF (.NOT.ALLOK) THEN
            WRITE (IUNIT8,102)
  102       FORMAT(' THE FOLLOWING NODES WERE NEVER READ:')
            DO 104 I=1,NUMNOD
                 IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)I
  103            FORMAT (' ',36X,I6)
  104       CONTINUE
            CALL PAUSE()
            STOP
       ENDIF
C
C  READ TRIANGULAR ELEMENTS
C
       READ (IUNIT7,*) NUMEL
       IF (NUMEL.GT.MXEL) THEN
            WRITE (IUNIT8,108) NUMEL
  108       FORMAT(/' INCREASE PARAMETER MAXEL TO BE AT LEAST EQUAL'
     +           /' TO THE NUMBER OF ELEMENTS (',I6,') AND RECOMPILE.')
            CALL PAUSE()
            STOP
       ENDIF
       DO 109 K=1,NUMEL
            CHECKE(K)=.FALSE.
  109  CONTINUE
       DO 200 K=1,NUMEL
C   (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.)
            IF (PLATES) THEN
                 READ (IUNIT7,*) I,(IFN(J),J=1,6)
                 CHECKE(I)=.TRUE.
                 DO 122 J=1,6
                      N=IFN(J)
                      IF ((N.LE.0).OR.(N.GT.NTOP).OR.
     +                     ((N.GT.NREALN).AND.(N.LE.N1000))) THEN
                           WRITE (IUNIT8,121) N
  121                      FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.')
                           CALL PAUSE()
                           STOP
                      END IF
                      IF (N.GT.NREALN) N=N-N1000+NREALN
                      NODES(J,I)=N
  122            CONTINUE
            ELSE IF (SHELLS) THEN
                 READ (IUNIT7,*) I,(IFN(J),J=1,3)
                 CHECKE(I)=.TRUE.
                 DO 127 J=1,3
                      N=IFN(J)
                      IF ((N.LE.0).OR.(N.GT.NTOP).OR.
     +                     ((N.GT.NREALN).AND.(N.LE.N1000))) THEN
                           WRITE (IUNIT8,121) N
                           CALL PAUSE()
                           STOP
                      END IF
                      IF (N.GT.NREALN) N=N-N1000+NREALN
                      NODES(J,I)=N
  127            CONTINUE
            END IF
  200  CONTINUE
       ALLOK=.TRUE.
       DO 201 I=1,NUMEL
            ALLOK=ALLOK.AND.CHECKE(I)
  201  CONTINUE
       IF (.NOT.ALLOK) THEN
            WRITE (IUNIT8,202)
  202       FORMAT (' THE FOLLOWING ELEMENTS WERE NEVER READ:')
            DO 204 I=1,NUMEL
                 IF (.NOT.CHECKE(I)) WRITE(IUNIT8,203)I
  203            FORMAT (' ',39X,I6)
  204       CONTINUE
            CALL PAUSE()
            STOP
       ENDIF
C
C   READ FAULT ELEMENTS
C
       READ (IUNIT7,*) NFL
       IF (NFL.GT.MXFEL) THEN
            WRITE (IUNIT8,220)NFL
  220       FORMAT (/' INCREASE PARAMETER MAXFEL TO BE AT LEAST EQUAL'
     +           /' TO THE NUMBER OF FAULTS (',I6,') AND RECOMPILE.')
            CALL PAUSE()
            STOP
       ENDIF
       OFFMAX=0.
       DO 222 I=1,NFL
            CHECKF(I)=.FALSE.
  222  CONTINUE
       DO 300 K=1,NFL
C
         IF (PLATES) THEN
C
            OFF=0.
            READ (IUNIT7,*,ERR=211) I,(IFN(J),J=1,6),(DIPS(L),L=1,3),
     +                      AZ1,AZ3,OFF
  211       CHECKF(I)=.TRUE.
            DO 213 J=1,6
                 N=IFN(J)
                 IF ((N.LE.0).OR.(N.GT.NTOP).OR.
     +                ((N.GT.NREALN).AND.(N.LE.N1000))) THEN
                      WRITE (IUNIT8,121) N
                      CALL PAUSE()
                      STOP
                 ENDIF
                 IF (N.GT.NREALN) N=N-N1000+NREALN
                 NODEF(J,I)=N
  213       CONTINUE
            DO 217 L=1,3
                 IF (ABS(DIPS(L)).GT.90.) THEN
                      WRITE(IUNIT8,215) DIPS(L)
  215                 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN',
     +                     ' RANGE OF -90. TO +90. DEGREES.'/
     +                     ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE',
     +                     ' HORIZONAL;'/
     +                     '  A + PREFIX (OR NONE) INDICATES A DIP',
     +                      ' TOWARD THE N1-N2-N3 SIDE;'/
     +                     '  A - PREFIX INDICATES A DIP TOWARD',
     +                      ' THE N6-N5-N4 SIDE.)')
                      CALL PAUSE()
                      STOP
                 ENDIF
                 IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L)
                 FDIP(L,I)=DIPS(L)*0.017453293
  217       CONTINUE
            IF ((ABS(AZ1).GT.361.).OR.(ABS(AZ3).GT.361.)) THEN
                 WRITE (IUNIT8,218) AZ1,AZ3
  218            FORMAT (' ILLEGAL ARGUMENT OF ',F10.4,' OR ',F10.4,
     +                '; SHOULD BE IN RANGE -360. TO +360. DEGREES.')
                 CALL PAUSE()
                 STOP
            ENDIF
            FAZ(1,I)=AZ1*0.017453293
            FAZ(2,I)=AZ3*0.017453293
            IF (OFF.LT.0.) THEN
                 WRITE (IUNIT8,219) OFF
  219            FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2,
     +                   ' FOR FAULT ELEMENT',I6/
     +                   ' OFFSETS MAY NOT BE NEGATIVE.')
                 CALL PAUSE()
                 STOP
            ENDIF
            OFFSET(I)=OFF
            OFFMAX=MAX(OFFMAX,OFF)
C
         ELSE IF (SHELLS) THEN
C
            OFF=0.
            READ (IUNIT7,*,ERR=221) I,(IFN(J),J=1,4),(DIPS(L),L=1,2),
     +                              OFF
  221       CHECKF(I)=.TRUE.
            DO 223 J=1,4
                 N=IFN(J)
                 IF ((N.LE.0).OR.(N.GT.NTOP).OR.
     +                ((N.GT.NREALN).AND.(N.LE.N1000))) THEN
                      WRITE (IUNIT8,121) N
                      CALL PAUSE()
                      STOP
                 ENDIF
                 IF (N.GT.NREALN) N=N-N1000+NREALN
                 NODEF(J,I)=N
  223       CONTINUE
            DO 227 L=1,2
                 IF (ABS(DIPS(L)).GT.90.) THEN
                      WRITE(IUNIT8,225) DIPS(L)
  225                 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN',
     +                     ' RANGE OF -90. TO +90. DEGREES.'/
     +                     ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE',
     +                     ' HORIZONAL;'/
     +                     '  A + PREFIX (OR NONE) INDICATES A DIP',
     +                      ' TOWARD THE N1-N2 SIDE;'/
     +                     '  A - PREFIX INDICATES A DIP TOWARD',
     +                      ' THE N3-N4 SIDE.)')
                      CALL PAUSE()
                      STOP
                 ENDIF
                 IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L)
                 FDIP(L,I)=DIPS(L)*0.017453293
  227       CONTINUE
            IF (OFF.LT.0.) THEN
                 WRITE (IUNIT8,229) OFF
  229            FORMAT (' ILLEGAL FAULT OFFSET OF ',1P,E10.2,
     +                   ' FOR FAULT ELEMENT',I6/
     +                   ' OFFSETS MAY NOT BE NEGATIVE.')
                 CALL PAUSE()
                 STOP
            ENDIF
            OFFSET(I)=OFF
            OFFMAX=MAX(OFFMAX,OFF)
C
         END IF
  300  CONTINUE
       ALLOK=.TRUE.
       DO 301 I=1,NFL
            ALLOK=ALLOK.AND.CHECKF(I)
  301  CONTINUE
       IF (.NOT.ALLOK) THEN
            WRITE(IUNIT8,302)
  302       FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:')
            DO 304 I=1,NFL
                 IF (.NOT.CHECKF(I)) WRITE(IUNIT8,303)I
  303            FORMAT(' ',36X,I6)
  304       CONTINUE
            CALL PAUSE()
            STOP
       ENDIF
       RETURN
       END SUBROUTINE GETNET
C
C
C
       SUBROUTINE INTERP (INPUT,ARRAY,
     +                          MAXDAT,NX,NY,
     +                          X1,DX,X2,
     +                          Y1,DY,Y2,
     +                          X,Y,
     +                   OUTPUT,VALUE)
C
C      OBTAINS A SINGLE VALUE FROM A .GRD FILE DATA SET
C
       INTEGER IC1,IC2,INPUT,IR1,IR2,MAXDAT,NX,NY,OUTPUT
       LOGICAL OUTSID,WAYOUT
       REAL ARRAY,BOT,DX,DY,FC,FR,TOP,VALUE,X,X1,X2,Y,Y1,Y2
       DIMENSION ARRAY(MAXDAT,MAXDAT)
C
       OUTSID=.FALSE.
       WAYOUT=.FALSE.
       IR1=((Y2-Y)/DY)+1.00001
       IR1=MAX(IR1,1)
       IR1=MIN(IR1,NY-1)
       IR2=IR1+1
       FR=((Y2-DY*(IR1-1))-Y)/DY
       IC1=((X-X1)/DX)+1.00001
       IC1=MAX(IC1,1)
       IC1=MIN(IC1,NX-1)
       IC2=IC1+1
       FC=(X-(X1+DX*(IC1-1)))/DX
       OUTSID=OUTSID.OR.(FR.LT.-0.01).OR.(FR.GT.1.01).OR.
     +                  (FC.LT.-0.01).OR.(FC.GT.1.01)
       WAYOUT=WAYOUT.OR.(FR.LT.-1.01).OR.(FR.GT.2.01).OR.
     +                  (FC.LT.-1.01).OR.(FC.GT.2.01)
       FR=MIN(1.,MAX(0.,FR))
       FC=MIN(1.,MAX(0.,FC))
       TOP=ARRAY(IR1,IC1)+FC*(ARRAY(IR1,IC2)-ARRAY(IR1,IC1))
       BOT=ARRAY(IR2,IC1)+FC*(ARRAY(IR2,IC2)-ARRAY(IR2,IC1))
       VALUE=TOP+FR*(BOT-TOP)
C
       IF (WAYOUT) THEN
            WRITE (*,10) X,Y
   10       FORMAT(' WARNING: Point (',ES10.3,',',ES10.3,') is '
     +             'slightly outside the .grd area.')
       ELSE IF (OUTSID) THEN
            WRITE (*,20) X,Y
   20       FORMAT(' ERROR: Point (',ES10.3,',',ES10.3,') is '
     +             'far outside the .grd area.')
            CALL PAUSE()
            STOP
       END IF
C     
       END SUBROUTINE INTERP
C
C
C
       SUBROUTINE PAUSE()
       WRITE(*,"(/' Press [Enter]...'\)")
       READ(*,*)
       RETURN
       END
C
C
C
      SUBROUTINE Prompt_for_Logical (prompt_text, default, answer)
C       Writes a line to the default (*) unit, with:
C         "prompt_text" ["default"]: 
C       and accepts an answer with a logical value.
C       If [Enter] is pressed with nothing preceding it,
C          then "answer" takes the value "default".
C       Note that prompt_text should usually end with '?'.
C       It can be more than 70 bytes long if necessary, but cannot
C          contain line breaks, tabs, or other formatting.
C       Trailing blanks in prompt_text are ignored.
        IMPLICIT NONE
        CHARACTER*(*), INTENT(IN) :: prompt_text
        LOGICAL, INTENT(IN) :: default
        LOGICAL, INTENT(OUT) :: answer
        CHARACTER*1 :: inbyte
        CHARACTER*3 :: yesno
        INTEGER :: blank_at, bytes, written
        LOGICAL :: finished
        bytes = LEN_TRIM(prompt_text)
        finished = .FALSE.
        DO WHILE (.NOT. finished)
            IF (default) THEN
                yesno = 'Yes'
            ELSE
                yesno = 'No'
            END IF
            written = 0
            DO WHILE ((bytes - written) > 70)
                blank_at = written + 
     +           INDEX(prompt_text((written+1):(written+70)),' ',.TRUE.)
                IF (blank_at < (written + 2)) blank_at = written + 70
                WRITE (*,"(' ',A)") prompt_text((written+1):blank_at)
                written = blank_at
            END DO
            WRITE (*,"(' ',A,' [',A']: '\)") 
     +          prompt_text((written+1):bytes), TRIM(yesno)
            finished = .TRUE.
            READ (*,"(A)") inbyte
            IF (LEN_TRIM(inbyte) == 0) THEN
                answer = default
            ELSE
                SELECT CASE (inbyte)
                CASE ('Y')
                    answer = .TRUE.
                CASE ('y')
                    answer = .TRUE.
                CASE ('T')
                    answer = .TRUE.
                CASE ('t')
                    answer = .TRUE.
                CASE ('R')
                    answer = .TRUE.
                CASE ('r')
                    answer = .TRUE.
                CASE ('O')
                    answer = .TRUE.
                CASE ('o')
                    answer = .TRUE.
                CASE ('N')
                    answer = .FALSE.
                CASE ('n')
                    answer = .FALSE.
                CASE ('F')
                    answer = .FALSE.
                CASE ('f')
                    answer = .FALSE.
                CASE ('W')
                    answer = .FALSE.
                CASE ('w')
                    answer = .FALSE.
                CASE DEFAULT
      WRITE(*,
     +"(' ERROR: Your response (',A,') was not recognized.')") inbyte
      WRITE(*,
     +"('       (Only the first letter of your answer is used.)')")
      WRITE(*,"('        To agree, enter Y, y, T, t, O, o, R, or r.')")
      WRITE(*,"('        To disagree, enter N, n, F, f, W, or w.')")
      WRITE(*,"('        Please try again:')")
                    finished = .FALSE.
                END SELECT
            END IF ! a byte was entered
        END DO ! until finished
      END SUBROUTINE Prompt_for_Logical
C
C
C
       SUBROUTINE PUTNET (INPUT,PLATES,SHELLS,
     +                          IUNITO,
     +                          BRIEF,DQDTDA,ELEV,FAZ,FDIP,
     +                          MXEL,MXFEL,MXNODE,N1000,
     +                          NFAKEN,NFL,NODEF,NODES,
     +                          NREALN,NUMEL,NUMNOD,OFFSET,
     +                          TITLE1,TLNODE,XNODE,YNODE,ZMNODE)
C
C   Note: Special version of PutNet; works with both types of .feg:
C   flat-Earth (PLATES = T) or round-Earth (SHELLS = T),
C   but does NOT work with FAULTS like the original PutNet!
C
C   WRITES FINITE ELEMENT GRID TO UNIT "IUNITO".
C
       CHARACTER*80 TITLE1
       INTEGER I,IPRINT,IUNITO,J,K,KM,KP,MXEL,MXFEL,MXNODE,
     +         N1000,NFAKEN,NFL,NODEF,NODES,NP,NREALN,NUMEL,NUMNOD
       LOGICAL BRIEF,INPUT,PLATES,SHELLS
       REAL AZ1,AZ3,DIPS,DQDTDA,ELEV,FAZ,FDIP,
     +      OFF2,OFFSET,SIDE2,TLNODE,
     +      XC,XM,XNODE,XP,YC,YM,YNODE,YP,ZMNODE
C
       DIMENSION DQDTDA(MXNODE),ELEV(MXNODE),
     +           FAZ(2,MXFEL),FDIP(3,MXFEL),NP(6),
     +           NODEF(6,MXFEL),NODES(6,MXEL),
     +           OFFSET(MXFEL),TLNODE(MXNODE),
     +           XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE)
       DIMENSION DIPS(3)
C
C      RETURN COORDINATES OF STRAIGHT MID-POINT NODES TO ZERO
C
       IF (PLATES) THEN
            DO 100 I=1,NUMNOD
                 DO 90 J=1,NUMEL
                      DO 80 K=4,6
                           IF (I.EQ.NODES(K,J)) THEN
                                KM=K-3
                                XM=XNODE(NODES(KM,J))
                                YM=YNODE(NODES(KM,J))
                                KP=MOD(K-3,3)+1
                                XP=XNODE(NODES(KP,J))
                                YP=YNODE(NODES(KP,J))
                                SIDE2=(XP-XM)**2+(YP-YM)**2
                                XC=(XM+XP)*0.5
                                YC=(YM+YP)*0.5
                                OFF2=(XNODE(I)-XC)**2+(YNODE(I)-YC)**2
                                IF (OFF2/SIDE2.LE.0.0001) THEN
                                     XNODE(I)=0.0
                                     YNODE(I)=0.0
                                END IF
                                GO TO 99
                           END IF
   80                 CONTINUE
   90            CONTINUE
   99            CONTINUE
  100       CONTINUE
       END IF
C
       WRITE(*,"(/ /' READY TO CREATE OUTPUT .FEG FILE ON UNIT ',I3)")
       WRITE(*,"(' Please enter a (new) filename with .feg ending: ')")
       WRITE (IUNITO,101) TITLE1
  101  FORMAT (A80)
C
       WRITE (IUNITO,102) NUMNOD, NREALN,NFAKEN,N1000,BRIEF
  102  FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)')
C
       DO 108 I=1,NUMNOD
            IF (I.LE.NREALN) THEN
                 IPRINT=I
            ELSE
                 IPRINT=N1000+(I-NREALN)
            ENDIF
            WRITE (IUNITO,106) I,XNODE(I),YNODE(I),
     +                         ELEV(I),DQDTDA(I),ZMNODE(I),
     +                         TLNODE(I)
  106       FORMAT (I5,1P,2E13.5,4E10.2)
  108  CONTINUE
C
       WRITE (IUNITO,110) NUMEL
  110  FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)')
C
       DO 200 I=1,NUMEL
            IF (PLATES) THEN
                 DO 150 K=1,6
                      IF (NODES(K,I).LE.NREALN) THEN
                           NP(K)=NODES(K,I)
                      ELSE
                           NP(K)=N1000+(NODES(K,I)-NREALN)
                      END IF
  150            CONTINUE
                 WRITE (IUNITO,160) I,(NP(K),K=1,6)
  160            FORMAT (I5,6I5)
            ELSE IF (SHELLS) THEN
                 DO 170 K=1,3
                      IF (NODES(K,I).LE.NREALN) THEN
                           NP(K)=NODES(K,I)
                      ELSE
                           NP(K)=N1000+(NODES(K,I)-NREALN)
                      END IF
  170            CONTINUE
                 WRITE (IUNITO,180) I,(NP(K),K=1,3)
  180            FORMAT (I5,3I5)
            END IF
  200  CONTINUE
C
       WRITE (IUNITO,210) NFL
  210  FORMAT (I10,' (NFL =  NUMBER OF FAULT ELEMENTS)')
C
       DO 300 I=1,NFL
            IF (PLATES) THEN
                 DO 220 K=1,6
                      IF (NODEF(K,I).LE.NREALN) THEN
                           NP(K)=NODEF(K,I)
                      ELSE
                           NP(K)=N1000+(NODEF(K,I)-NREALN)
                      END IF
  220            CONTINUE
                 DO 230 K=1,3
                      DIPS(K)=FDIP(K,I)
                      DIPS(K)=DIPS(K)*57.2958
                      IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180.
  230            CONTINUE
                 AZ1=FAZ(1,I)*57.2958
                 AZ3=FAZ(2,I)*57.2958
                 WRITE (IUNITO,250) I,(NP(K),K=1,6),
     +               (DIPS(K),K=1,3),AZ1,AZ3,OFFSET(I)
  250            FORMAT (I5,6I6,3F6.1,2F7.1,1P,E10.2)
            ELSE IF (SHELLS) THEN
                 DO 260 K=1,4
                      IF (NODEF(K,I).LE.NREALN) THEN
                           NP(K)=NODEF(K,I)
                      ELSE
                           NP(K)=N1000+(NODEF(K,I)-NREALN)
                      END IF
  260            CONTINUE
                 DO 270 K=1,2
                      DIPS(K)=FDIP(K,I)
                      DIPS(K)=DIPS(K)*57.2958
                      IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180.
  270            CONTINUE
                 WRITE (IUNITO,280) I,(NP(K),K=1,4),
     +               (DIPS(K),K=1,2),OFFSET(I)
  280            FORMAT (I5,4I6,2F6.1,1P,E10.2)
            END IF
  300  CONTINUE
C
       RETURN
       END SUBROUTINE PUTNET
C
C
C
       SUBROUTINE READN (INPUT,IUNITP,IUNITT,N,
     +                  OUTPUT,VECTOR)
C
C   A UTILITY ROUTINE DESIGNED TO PERMIT -FAULTS- INPUT FILES
C      TO ALSO BE USED BY -PLATES-, WHICH EXPECTS MORE NUMBERS
C      IN SOME RECORDS.
C   THIS ROUTINE ATTEMPTS TO READ 'N' FLOATING-POINT VALUES
C      (USING * FORMAT) FROM THE NEXT RECORD ON DEVICE 'IUNITP'.
C   IF ANYTHING GOES WRONG, THE MISSING VALUES ARE SET TO ZERO.
C
       CHARACTER*1 C
       CHARACTER*80 LINE
       INTEGER I,IUNITP,IUNITT,N,NUMBER
       LOGICAL ANYIN,DOTTED,EXPON,INPUT,OUTPUT,SIGNED
       REAL VECTOR
       DIMENSION VECTOR(N)
C
       READ (IUNITP,1) LINE
    1  FORMAT (A80)
C
       NUMBER=0
       ANYIN=.FALSE.
       EXPON=.FALSE.
       SIGNED=.FALSE.
       DOTTED=.FALSE.
       DO 10 I=1,80
            C=LINE(I:I)
            IF ((C.EQ.' ').OR.(C.EQ.',').OR.(C.EQ.'/')) THEN
                 SIGNED=.FALSE.
                 EXPON=.FALSE.
                 DOTTED=.FALSE.
                 IF (ANYIN) THEN
                      NUMBER=NUMBER+1
                      ANYIN=.FALSE.
                 ENDIF
            ELSE IF ((C.EQ.'+').OR.(C.EQ.'-')) THEN
                 IF (SIGNED) THEN
                      GO TO 50
                 ELSE
                      SIGNED=.TRUE.
                 ENDIF
            ELSE IF ((C.EQ.'E').OR.(C.EQ.'D').OR.
     +               (C.EQ.'e').OR.(C.EQ.'d')) THEN
                 IF (EXPON) THEN
                      GO TO 50
                 ELSE
                      EXPON=.TRUE.
                      SIGNED=.FALSE.
                      DOTTED=.TRUE.
                 ENDIF
            ELSE IF (C.EQ.'.') THEN
                 IF (DOTTED) THEN
                      GO TO 50
                 ELSE
                      DOTTED=.TRUE.
                 ENDIF
            ELSE IF ((C.EQ.'0').OR.(C.EQ.'1').OR.(C.EQ.'2').OR.
     +               (C.EQ.'3').OR.(C.EQ.'4').OR.(C.EQ.'5').OR.
     +               (C.EQ.'6').OR.(C.EQ.'7').OR.(C.EQ.'8').OR.
     +               (C.EQ.'9')) THEN
                 SIGNED=.TRUE.
                 ANYIN=.TRUE.
            ELSE
                 GO TO 50
            ENDIF
   10  CONTINUE
       IF (ANYIN) NUMBER=NUMBER+1
C
   50  IF (NUMBER.EQ.0) THEN
            WRITE (IUNITT,91) N,LINE
   91       FORMAT (/' ERR0R: A LINE OF PARAMETER INPUT WHICH',
     +               ' WAS SUPPOSED TO CONTAIN 1-',I2,' NUMBERS'/
     +               ' COULD NOT BE INTERPRETED.  LINE FOLLOWS:'/
     +               ' ',A80)
            CALL PAUSE()
            STOP
       ELSE
            NUMBER=MIN(NUMBER,N)
            BACKSPACE IUNITP
            READ (IUNITP,*) (VECTOR(I),I=1,NUMBER)
            IF (NUMBER.LT.N) THEN
                 DO 99 I=NUMBER+1,N
                      VECTOR(I)=0.
   99            CONTINUE
            ENDIF
       ENDIF
       RETURN
       END SUBROUTINE READN