SUBROUTINE GRIDER (N,X,Y,Z,GRDSIZ,IX,IY,ZZ,NX,NY,
     *  MXMN,XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN)
C
C
C    THIS ROUTINE ENABLES THE USER TO GRID IRREGULARY SPACED
C      THREE-DIMENSIONAL DATA INTO EQUALLY SPACED ELEVATIONS
C      (A WEIGHTED MOVING AVERAGE IS USED)
C
C
C    CALLING SEQUENCE:
C
C         CALL GRIDER (N,X,Y,Z,GRDSIZ,IX,IY,ZZ,NX,NY,
C        *  MXMN,XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN)
C
C
C    N       - (INPUT) THE NUMBER OF (X,Y,Z) TRIPLICATES TO BE GRIDDED
C
C    (X,Y,Z) - (INPUT) THE COORDINATES OF THE IRREGULARY SPACED
C                THREE-DIMENSIONAL DATA TO BE GRIDDED
C                (NOTE - THE X ARRAY CONTAINS THE X VALUES
C                        THE Y ARRAY CONTAINS THE Y VALUES
C                        THE Z ARRAY CONTAINS THE Z VALUES)
C
C    GRDSIZ  - (INPUT)^1 THE DESIRED LENGTH OF EACH SQUARE IN THE GRID
C                (NOTE - IF GRDSIZ .LE. 0 CALCULATE THE BEST GRDSIZ)
C              (OUTPUT) THE LENGTH OF EACH SQUARE USED IN THE GRID
C                (WARNING - GRDSIZ MAY BE CHANGED BY GRIDER)
C
C    IX      - (INPUT)^2 THE NUMBER OF GRID LINES WANTED IN THE X
C                DIRECTION (NOTE - IF IX .LE. 1 CALCUATE THE BEST IX)
C              (OUTPUT) THE NUMBER OF GRID LINES USED IN THE X DIRECTION
C                (WARNING - IX MAY BE CHANGED BY GRIDER)
C
C    IY      - (INPUT)^3 THE NUMBER OF GRID LINES WANTED IN THE Y
C                DIRECTION (NOTE - IF IY .LE. 1 CALCUATE THE BEST IY)
C              (OUTPUT) THE NUMBER OF GRID LINES USED IN THE Y DIRECTION
C                (WARNING - IY MAY BE CHANGED BY GRIDER)
C
C    ZZ      - (OUTPUT) THE TWO-DIMENSIONAL ARRAY THAT CONTAINS THE
C                GRIDDED DATA
C
C    NX      - (INPUT) THE X DIMENSION OF THE ZZ ARRAY
C
C    NY      - (INPUT) THE Y DIMENSION OF THE ZZ ARRAY
C
C    MXMN    - (INPUT) FLAG TO CALCULATE X AND Y MAXIMUMS AND MINUMIMS
C                IF MXMN .EQ. 0, CALCULATE XMAX, XMIN, YMAX, AND YMIN
C                IF MXMN .NE. 0, DON'T CALCULTE XMAX, XMIN, YMAX, AND
C                  YMIN USE WHAT THE CALLER SENDS
C
C    XMAX    - (INPUT) THE MAXIMUM X VALUE
C              (OUTPUT) THE MAXIMUM X VALUE (WARNING - XMAX MAY BE
C                CHANGED BY GRIDER)
C
C    XMIN    - (INPUT) THE MINIMUM X VALUE
C              (OUTPUT) THE MINIMUM X VALUE (WARNING - XMIN MAY BE
C                CHANGED BY GRIDER)
C
C    YMAX    - (INPUT) THE MAXIMUM Y VALUE
C              (OUTPUT) THE MAXIMUM Y VALUE (WARNING - YMAX MAY BE
C                CHANGED BY GRIDER)
C
C    YMIN    - (INPUT) THE MINIMUM Y VALUE
C              (OUTPUT) THE MINIMUM Y VALUE (WARNING - YMIN MAY BE
C                CHANGED BY GRIDER)
C
C    ZMAX    - (OUTPUT) THE MAXIMUM Z VALUE
C
C    ZMIN    - (OUTPUT) THE MINIMUM Z VALUE
C
C
C    NOTE:  REQUIRED SUBROUTINES - MAXMIN
C
C
      DIMENSION X(1),Y(1),Z(1),ZZ(NX,NY),D(0:6),N1(6)
      EQUIVALENCE (DX,XP),(DY,YP)
      REAL MINX,MINY,MINZ,MAXZ
      DATA D(0)/1.0E-37/
C
C   BOUNDARIES AREN'T DEFINED BY THE CALLER, SO GO GET THEM
C
1000  IF (MXMN .NE. 0) GOTO 1100
      CALL MAXMIN (X,N,XMAX,XMIN)
      CALL MAXMIN (Y,N,YMAX,YMIN)
C
C   COMPUTE THE MAXIMUM AND MINIMUM Z VALUES
C   AND GET SOME DELTA VALUES
C
1100  CALL MAXMIN (Z,N,ZMAX,ZMIN)
      DX = XMAX - XMIN
      DY = YMAX - YMIN
C
C   GRID SPACING SPECIFIED
C
      IF (GRDSIZ .LE. 0.) GOTO 1200
      IX = DX / GRDSIZ + 1.5
      IY = DY / GRDSIZ + 1.5
      GOTO 1700
C
C   NUMBER OF X GRID LINES SPECIFIED
C
1200  IF (IX .LE. 1) GOTO 1400
1300  GRDSIZ = DX / FLOAT(IX - 1)
      IY = DY / GRDSIZ + 1.5
      GOTO 1700
C
C   NUMBER OF Y GRID LINES SPECIFIED
C
1400  IF (IY .LE. 1) GOTO 1600
1500  GRDSIZ = DY / FLOAT(IY - 1)
      IX = DX / GRDSIZ + 1.5
      GOTO 1700
C
C   NO SPECIFICATIONS GIVEN
C   PROGRAM COMPUTES BEST GRID SIZE
C
1600  GRDSIZ = SQRT((DX * DY) / FLOAT(N))
      IX = DX / GRDSIZ + 1.5
      IY = DY / GRDSIZ + 1.5

1700  IF (IX .LE. NX) GOTO 1800
      IX = NX
      GOTO 1200

1800  IF (IY .LE. NY) GOTO 1900
      IY = NY
      GOTO 1500
C
C   PERFORM INTERPOLATION
C
1900  MINX = XMIN + (DX - FLOAT(IX - 1) * GRDSIZ ) / 2.
      MINY = YMIN + (DY - FLOAT(IY - 1) * GRDSIZ ) / 2.
      CLOSE = GRDSIZ / 25.
      DUM = (ZMAX - ZMIN) / 2.
      MAXZ = ZMAX + DUM
      MINZ = ZMIN - DUM
C
C   LOOP THROUGH THE GRID 'X' BY 'X'
C
2000  DO 2800 I = 1,IX
         XP = MINX + (FLOAT(I - 1)) * GRDSIZ
C
C   LOOP THROUGH EACH 'X', 'Y' BY 'Y'
C
         DO 2700 J = 1,IY
            YP = MINY + (FLOAT(J - 1)) * GRDSIZ
            L = 0
C
C   TEST EACH DATA POINT TO FIND THE NEAREST ONES
C
            DO 2500 M = 1,N
               DIST = (XP - X(M))**2 + (YP - Y(M))**2
               IF (DIST .GE. CLOSE) GOTO 2100
C
C   AN ORIGINAL DATA POINT LIES WITHIN A 1/25TH OF A GRID
C   SPACING, INTERPOLATION STOPPED FOR THIS POINT
C
               ZZ(I,J) = Z(M)
               GOTO 2700
C
C   FIND SIX NEAREST POINTS
C
2100           IF (DIST .GE. D(L)) GOTO 2400
               IF (L .LT. 6) L = L + 1
               L1 = L - 1
               KK = L1
               DO 2200 K = L1,1,-1
                  IF (DIST .GE. D(K)) GOTO 2300
                  D(K+1) = D(K)
                  N1(K+1) = N1(K)
2200              KK=K
2300           D(KK) = DIST
               N1(KK) = M
               GOTO 2500
2400           IF(L .GE. 6) GOTO 2500
               L = L + 1
               D(L) = DIST
               N1(L) = M
2500           CONTINUE
C
C   COMPUTE WEIGHTED MOVING AVERAGE
C
            DUM2 = 0.
            PSUM = 0.
            DO 2600 M = 1,6
               DUM = 1. / (SQRT(D(M)))
               PSUM = PSUM + DUM
2600           DUM2 = DUM2 + Z(N1(M)) * DUM
            DUM1 = Z(N1(1))
            DUM = (DUM1 + DUM2 / PSUM) / 2.
            IF (DUM .LE. MAXZ .AND. DUM .GE. MINZ) DUM1 = DUM
            ZZ(I,J) = DUM1
2700        CONTINUE
2800     CONTINUE
      RETURN
      END