
      PROGRAM WRARFI
C
C This program, given a specified map dataset, creates a data file
C containing area identifiers for a large set of positions on the
C globe.  This data file is to be used by the EZMAPB routine MPLNAM
C to determine an area identifier to use in cases in which the current
C mapped area does not include any portion of any boundary line from
C the selected dataset.
C
C The grid points used are generated as follows:  First, the grid is
C divided into two halves - the left half and the right half - each
C of which is a square.  The center points of each half are used.
C Then, each of those halves is subdivided into four square quarters
C and the eight center points of those are used.  Then, each of the
C eight squares is subdivided into four square quarters and the
C thirty-two center points of those are used.  At each subsequent
C step, the squares from the previous step are subdivided into four
C smaller squares and the center points of those squares are used.
C This order tends to minimize the number of points that have to be
C examined in order to find one giving the desired information: one
C looks first at a very coarse grid, then at a finer one, then at a
C still finer one, and so on.
C
C The area identifiers are written in such a way as to keep the file
C small: each entry is just a pair of numbers, "m" and "n", interpreted
C to mean "m repetitions of the area identifier n", and there are no
C carriage returns.
C
C Declare the number of subdivisions to be carried out.  I'm currently
C using NDIV=13 to get a pretty dense set of points.  (The final grid
C of points will have dimensions 16384x8192, and the spacing of these
C points on the globe will be about .022 degrees - about 1.5 miles at
C the equator.)  At this level, the program runs about 15 minutes on
C a 3.06 GHz iMac.  Upping the value of NDIV to 14 would make it run
C about an hour and create an output file about four times as large.
C
        PARAMETER (NDIV=13)
C
C Declare the length of the area map to be used.
C
        PARAMETER (LAMA=10000000)
C
C Declare the area-map array and the arrays in which ARGTAI will return
C area-identifier information.
C
        DIMENSION IAMA(LAMA),IAAI(3),IAGI(3)
C
C Declare a variable in which to form names of output files.
C
        CHARACTER*71 FLNM
C
C Declare an array to use as an output buffer.
C
        CHARACTER*1 CHRS(1024)
C
C Get the name of the file to be used and check for errors.
C
        PRINT * , ' '
        PRINT * , 'Enter the name of the file to be used:'
        READ (*,'(A64)') FLNM
C
        IF (FLNM.EQ.' ') THEN
          PRINT * , 'ERROR - Input file name is blank.'
          STOP
        END IF
C
        LFNM=MPILNB(FLNM)
C
C Open the output data file.
C
        FLNM(LFNM+1:LFNM+6)='.areas'
        FLNM(LFNM+7:LFNM+7)=CHAR(0)
C
        CALL NGOFWO (FLNM(1:LFNM+7),IFDE,ISTA)
C
        IF (ISTA.NE.0) THEN
          PRINT * , 'ERROR - Bad status from NGOFWO: ',ISTA
          STOP
        END IF
C
C Initialize the count of characters sent to the buffer so far.
C
        NCHR=0
C
C Open GKS.
C
        CALL OPNGKS
C
C Tell EZMAP to use a cylindrical equidistant projection and the
C maximal useful area produced by that projection.
C
        CALL MAPROJ ('CE',0.,0.,0.)
        CALL MAPSET ('MA',0.,0.,0.,0.)
C
C Initialize EZMAP.
C
        CALL MAPINT
C
C Initialize the area map.
C
        CALL ARINAM (IAMA,LAMA)
C
C Put the contents of the outline dataset specified by "FLNM" into the
C area map.
C
        CALL MPLNAM (FLNM(1:LFNM),5,IAMA)
C
C Preset the following: IFLG makes ARGTAI call GETSET only the first
C time it's called; IAID is the area identifier for the last point
C examined; NTMS is the number of times in a row the same area
C identifier has turned up; NCHK is used for an error check on the
C number of points processed.
C
        IFLG=1
        IAID=0
        NTMS=0
        NCHK=0
C
C Loop on the number of times the grid is to be subdivided.
C
        DO 104 IDIV=0,NDIV
C
C Compute the number of rows and columns in the grid at this level of
C division and the resulting number of center points of grid boxes for
C which area identifiers will be generated.
C
          NROW=2**IDIV
          NCOL=2*NROW
          NPTS=NROW*NCOL
C
C Reduce NCHK by the number of points that should be processed at this
C level of subdivision.
C
          NCHK=NCHK-NPTS
C
C Loop through the rows of the grid.  Compute associated latitude.
C
          DO 103 J=1,NROW
C
            RLAT=-90.+(REAL(J)-.5)*(180./REAL(NROW))
C
C Loop through the columns of the grid.  Computed associated longitude.
C
            DO 102 I=1,NCOL
C
              RLON=-180.+(REAL(I)-.5)*(360./REAL(NCOL))
C
C Save the area identifier found for the last grid point.
C
              IAIP=IAID
C
C Find the user coordinates associated with the current grid point.
C
              CALL MAPTRA (RLAT,RLON,XVAL,YVAL)
C
C Get area-identifier information associated with the point.
C
              CALL ARGTAI (IAMA,XVAL,YVAL,IAAI,IAGI,3,NGPS,IFLG)
              IFLG=0
C
C Search the returned information for the new group-1 area identifier.
C
              IAID=0
C
              DO 101 K=1,NGPS
                IF (IAGI(K).EQ.1) IAID=IAAI(K)
  101         CONTINUE
C
C If the new group-1 area identifier is different from the last one,
C dump out information about the last one and reset the counter.
C
              IF (IAID.NE.IAIP) THEN
C
                IF (NTMS.NE.0) THEN
                  CALL WRNUMB (IFDE,CHRS,1024,NCHR,NTMS)
                  CALL WRCHAR (IFDE,CHRS,1024,NCHR, ' ')
                  CALL WRNUMB (IFDE,CHRS,1024,NCHR,IAIP)
                  CALL WRCHAR (IFDE,CHRS,1024,NCHR, ' ')
                END IF
C
                NCHK=NCHK+NTMS
                NTMS=0
C
              END IF
C
C Bump the count of the number of times the current area identifier
C has been seen.
C
              NTMS=NTMS+1
C
C End of loop through columns.
C
  102       CONTINUE
C
C End of loop through rows.
C
  103     CONTINUE
C
C End of loop on subdivisions of the grid.
C
  104   CONTINUE
C
C If the count is non-zero, dump out the last two numbers.
C
        IF (NTMS.NE.0) THEN
          CALL WRNUMB (IFDE,CHRS,1024,NCHR,NTMS)
          CALL WRCHAR (IFDE,CHRS,1024,NCHR, ' ')
          CALL WRNUMB (IFDE,CHRS,1024,NCHR,IAIP)
          CALL WRCHAR (IFDE,CHRS,1024,NCHR, ' ')
        END IF
C
C If there's something left in the output buffer, dump it.
C
        IF (NCHR.GT.0) THEN
          CALL NGWRCH (IFDE,CHRS,NCHR,ISTA)
          IF (ISTA.LE.0) THEN
            PRINT * , 'TESTIT - BAD STATUS FROM NGWRCH = ',ISTA
            STOP
          END IF
        END IF
C
        NCHK=NCHK+NTMS
C
C Check for a discrepancy in the number of points processed.
C
        IF (NCHK.NE.0) PRINT * , 'WRONG NUMBER OF POINTS PROCESSED'
C
C Close the file just written.
C
        CALL NGCLFI (IFDE)
C
C Close GKS.
C
        CALL CLSGKS
C
C Done.
C
        STOP
C
      END



      SUBROUTINE WRNUMB (IFDE,CHRS,MCHR,NCHR,INUM)
C
        CHARACTER*1 CHRS(MCHR)
C
C Translate the integer INUM to a character form and put it in the
C buffer CHRS.
C
        CHARACTER*1 CTMP(20)
C
C Copy the integer to a local variable.
C
        ITMP=INUM
C
C Generate digits from right to left in the temporary buffer CTMP.
C
        NDGT=0
C
  101   NDGT=NDGT+1
C
        IF (NDGT.GT.20) THEN
          PRINT * , 'WRNUMB - INTEGER TOO BIG - STOP'
          STOP
        END IF
C
        CTMP(NDGT)=CHAR(ICHAR('0')+MOD(ITMP,10))
C
        ITMP=ITMP/10
        IF (ITMP.NE.0) GO TO 101
C
C Send the digits to the output buffer in correct order.
C
        DO 102 IDGT=NDGT,1,-1
          CALL WRCHAR (IFDE,CHRS,MCHR,NCHR,CTMP(IDGT))
  102   CONTINUE
C
C Done.
C
        RETURN
C
      END



      SUBROUTINE WRCHAR (IFDE,CHRS,MCHR,NCHR,ICHR)
C
        CHARACTER*1 CHRS(MCHR),ICHR
C
C Put the character ICHR in the buffer CHRS.
C
        NCHR=NCHR+1
        CHRS(NCHR)=ICHR
C
C If the buffer is now full, dump it to the file specified by IFDE.
C
        IF (NCHR.GE.MCHR) THEN
          CALL NGWRCH (IFDE,CHRS,NCHR,ISTA)
          IF (ISTA.LE.0) THEN
            PRINT * , 'WRCHAR - BAD STATUS FROM NGWRCH = ',ISTA
            STOP
          END IF
          NCHR=0
        END IF
C
C Done.
C
        RETURN
C
      END
