
      PROGRAM WRARSZ
C
C This program reads the area file for a specified EZMAP outline dataset
C and writes, to standard output, estimates of the sizes of named areas
C of the globe.
C
C The grid points involved are as follows:  First, the grid is divided
C into two halves - the left half and the right half - each of which is
C a square.  The center points of each half are used.  Then, each of
C those halves is subdivided into four square quarters and the eight
C center points of those are used.  Then, each of the eight squares is
C subdivided into four square quarters and the thirty-two center points
C of those are used.  At each subsequent step, the squares from the
C previous step are subdivided into four smaller squares and the center
C points of those squares are used.  This order tends to minimize the
C number of points that have to be examined in order to find one giving
C the desired information: one looks first at a very coarse grid, then
C at a finer one, then at a 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.  Beware: if this
C number is bigger than the value of NDIV used when the area file was
C created, all the sizes reported will be zeroes.
C
        PARAMETER (NDIV=13)
C
C Declare the largest number of areas the program can handle.
C
        PARAMETER (MNAI=8000)
C
C Declare variables with which to keep track of area estimates.
C
        DOUBLE PRECISION DTOR,PI,HBOX,WBEQ,AGBX,AOTE,AREA(MNAI)
C
C Declare the type of an EZMAPB function used below.
C
        CHARACTER*64 MPNAME
C
C Declare a variable in which to form a file name.
C
        CHARACTER*71 FLNM
C
C Declare an array to use as an input buffer.
C
        CHARACTER*1 CHRS(512)
C
C Define double-precision values of the constants PI/180 and PI.
C
        DATA DTOR / .017453292519943D0 /
        DATA PI   / 3.14159265358979D0 /
C
C Get the name of the file to be used and check for error.
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 Make EZMAPB read the name file for the dataset specified by "FLNM"
C so that we can use calls to retrieve information from it.
C
        CALL MPLNRI (FLNM(1:LFNM))
C
C Open the area file and initialize variables involved in reading it.
C
        FLNM(LFNM+1:LFNM+6)='.areas'
        FLNM(LFNM+7:LFNM+7)=CHAR(0)
C
        CALL NGOFRO (FLNM(1:LFNM+7),IFDE,ISTA)
C
        IF (ISTA.NE.0) THEN
          PRINT * , 'ERROR - Bad status from NGOFRO: ',ISTA
          STOP
        END IF
C
        NTMS=0
        MCHR=0
        NCHR=0
C
C Loop on the number of times the grid is to be subdivided.
C
        DO 105 IDIV=0,NDIV
C
C Compute the number of rows and columns in the grid at this level of
C subdivision.
C
          NROW=2**IDIV
          NCOL=2*NROW
C
C Compute the height of a grid box (in radians) and the width of a grid
C box at the equator (in kilometers), and zero all area estimates.
C
          HBOX=(180.D0/DBLE(NROW))*DTOR
          WBEQ=40000.D0/DBLE(NCOL)

          AOTE=0.D0
C
          DO 101 I=1,MNAI
            AREA(I)=0.D0
  101     CONTINUE
C
C Loop through the rows of the grid.  Compute associated latitude (in
C double precision).
C
          DO 104 J=1,NROW
C
            DLAT=-90.D0+(DBLE(J)-.5D0)*(180.D0/DBLE(NROW))
C
C Compute the area of a grid box at the current latitude.
C
            AGBX=WBEQ*WBEQ*COS(DTOR*DLAT)
C
C Loop through the columns of the grid.
C
            DO 103 I=1,NCOL
C
C Update the estimate of the area of the earth.
C
              AOTE=AOTE+AGBX
C
C If the last area identifier had been used the correct number of times,
C read a new pair of numbers from the area file (a count and a new area
C identifier).
C
              IF (NTMS.EQ.0) THEN
                CALL MDRDNM (IFDE,CHRS,512,MCHR,NCHR,NTMS)
                IF (MCHR.EQ.0) GO TO 106
                CALL MDRDNM (IFDE,CHRS,512,MCHR,NCHR,IAID)
                IF (MCHR.EQ.0) GO TO 106
              END IF
C
C Update the area estimates for the area and for its parents.
C
              ITMP=IAID
C
  102         IF (ITMP.GE.1.AND.ITMP.LE.MNAI) THEN
                AREA(ITMP)=AREA(ITMP)+AGBX
                ITMP=MPIPAR(ITMP)
                GO TO 102
              END IF
C
C Decrement the counter saying how many times to use the current area
C identifier.
C
              NTMS=NTMS-1
C
C End of loop through columns.
C
  103       CONTINUE
C
C End of loop through rows.
C
  104     CONTINUE
C
C End of loop on subdivisions of the grid.
C
  105   CONTINUE
C
C Close the area file.
C
  106   CALL NGCLFI (IFDE)
C
C Print the area estimates.
C
        PRINT * , ' '
        PRINT '(''  Estimated size of the surface of earth:'',21X,F16.5)
     +' , AOTE
        PRINT * , ' '
        PRINT '(''  Estimated size of various named areas:'')'
        PRINT * , ' '
C
        DO 107 I=1,MNAI
          IF (MPIATY(I).NE.0) THEN
            PRINT '(''  '',A60,F16.5)' , MPNAME(I),AREA(I)
          END IF
  107   CONTINUE
C
        PRINT * , ' '
C
C Done.
C
        STOP
C
      END
