PROGRAM GEOGRD

***********************************************************************
*                                                                     *
* PURPOSES: 1) To extract a GEOID subgrid from a larger GEOID grid.   *
*                                                                     *
*           2) To translate a GEOID grid file to binary from the      *
*              ASCII transfer format and visa versa.                  *
*                                                                     *
*           3) To translate a GEOID grid file from binary or the      *
*              ASCII transfer format to the ASCII graphics format.    *
*                                                                     *
*           4) To print information about a GEOID grid to a file.     *
*                                                                     *
*              The program GEOID reads a binary grid of geoid heights.*
*              The first record in a grid file consists of header     *
*              information.  All the other records consist of REAL*4  *
*              grid numbers.  The grid files are unformatted and      *
*              direct access.                                         *
*                                                                     *
*              The GEOID grid consist of geoid heights in meters.     *
*              The geoid heights were obtained from a complete        *
*              expansion of the spherical harmonic coefficients       *
*              from the Ohio State University model known as OUS89B.  *
*              The names of the grids are consist of the location of  *
*              the grid (e.g. 'conus' or 'alaska') plus an extension. *
*              The extension for the GEOID grids of geoid height is   *
*              '.geo'.                                                *
*                                                                     *
*              The ASCII transfer file is used for porability.        *
*              The first two records in the ASCII file contain the    *
*              header information.  The rest of the records are       *
*              formatted as 6F12.6.                                   *
*                                                                     *
*              This program is intended to be used with GEOID         *
*              development.                                           *
*                                                                     *
*              All grids start and finish at even degrees.            *
*              NOTE that there must be minimum number of degrees      *
*              between the minimum and maximum longitude.  This is    *
*              explained further in subroutine NLON.                  *
*                                                                     *
* VERSION CODE:  0.90                                                 *
*                                                                     *
* VERSION DATE:  SEPTEMBER 1, 1990                                    *
*                                                                     *
* AUTHOR:   ALICE R. DREW                                             *
*           WARREN T. DEWHURST, PH.D.                                 *
***********************************************************************

***********************************************************************
*                                                                     *
*                  DISCLAIMER                                         *
*                                                                     *
*   THIS PROGRAM AND SUPPORTING INFORMATION IS FURNISHED BY THE       *
* GOVERNMENT OF THE UNITED STATES OF AMERICA, AND IS ACCEPTED AND     *
* USED BY THE RECIPIENT WITH THE UNDERSTANDING THAT THE UNITED STATES *
* GOVERNMENT MAKES NO WARRANTIES, EXPRESS OR IMPLIED, CONCERNING THE  *
* ACCURACY, COMPLETENESS, RELIABILITY, OR SUITABILITY OF THIS         *
* PROGRAM, OF ITS CONSTITUENT PARTS, OR OF ANY SUPPORTING DATA.       *
*                                                                     *
*   THE GOVERNMENT OF THE UNITED STATES OF AMERICA SHALL BE UNDER NO  *
* LIABILITY WHATSOEVER RESULTING FROM ANY USE OF THIS PROGRAM.  THIS  *
* PROGRAM SHOULD NOT BE RELIED UPON AS THE SOLE BASIS FOR SOLVING A   *
* PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO PERSON   *
* OR PROPERTY.                                                        *
*                                                                     *
*   THIS PROGRAM IS PROPERTY OF THE GOVERNMENT OF THE UNITED STATES   *
* OF AMERICA.  THEREFORE, THE RECIPIENT FURTHER AGREES NOT TO ASSERT  *
* PROPRIETARY RIGHTS THEREIN AND NOT TO REPRESENT THIS PROGRAM TO     *
* ANYONE AS BEING OTHER THAN A GOVERNMENT PROGRAM.                    *
*                                                                     *
***********************************************************************

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* NMAX is the maximum size of the buffer.  It is at least as large
* as the largest number of columns in the input grid

      REAL VRSION
      PARAMETER (VRSION = 1.00E0)
      INTEGER NMAX
      PARAMETER (NMAX = 1024)

      REAL XMIN, XMAX, YMIN, YMAX
      INTEGER NCFRST, NCLAST, NRFRST, NRLAST
      INTEGER KIN, KOUT
      CHARACTER*56 RIDENT
      CHARACTER*32 FOUTG
      CHARACTER*28 BASEIN
      CHARACTER*8 PGM
      LOGICAL FLAG

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

* Initialize the input/output unit numbers

      LUIN = 5
      LUOUT = 6
      NING = 11
      NOUTG = 12

* Header

      CALL HEADR (VRSION)

* Open the input and output files

      FLAG = .FALSE.
      CALL IFILES (KIN, FLAG, BASEIN, RIDENT, PGM)
      IF (FLAG) GOTO 9000

      CALL OFILES (KIN, KOUT, FOUTG, BASEIN, RIDENT, PGM)

* Get the variables for the new grid file.

      CALL NPARMS (KOUT, XMIN, XMAX, YMIN, YMAX,
     +             NRFRST, NRLAST, NCFRST, NCLAST)

* Get maximum and minimum Z values in new grid area

      CALL ZEXTRM (KIN, KOUT, NCFRST, NCLAST, NRFRST, NRLAST)

*** Copy from input file to output file

      IF (KOUT .NE. 0) THEN
        CALL TONEW (KIN, KOUT, XMIN, XMAX, YMIN, YMAX,
     +              NCFRST, NCLAST, NRFRST, NRLAST, FOUTG)
      ENDIF

* Close all files

 9000 CONTINUE
      CLOSE (NING,STATUS='KEEP')
      IF (.NOT. FLAG) THEN
        CLOSE (NOUTG,STATUS='KEEP')
      ENDIF
      WRITE (LUOUT,9010)
 9010 FORMAT (' End of program GEOGRD')
      STOP
      END
      SUBROUTINE HEADR (VRSION)

*** This subroutine prints the header information

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

      REAL VRSION
      CHARACTER*1 ANS

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

      WRITE (LUOUT,920)
  920 FORMAT (15X, '                   Welcome', /,
     + 15X,        '                   to  the', /,
     + 15X,        '          National Geodetic Survey', /,
     + 15X,        '      Geoid Grid Manipulation program.', //,
     + 15X,        '          For use when a GEOID grid', /,
     + 15X,        '            need to be translated', /,
     + 15X,        '          between ASCII and binary', /,
     + 15X,        '                     or', /,
     + 15X,        ' when a grid covering a smaller areal extent', /,
     + 15X,        '    need to be extracted from a standard', /,
     + 15X,        '                 GEOID grid.', /,
     + 15X,        '                     or', /,
     + 15X,        '          for obtaining information', /,
     + 15X,        '             about a GEOID grid.', /)

      WRITE (LUOUT,930) VRSION
  930 FORMAT (15X, '               (Version', F5.2, ')', /,
     + 15X,        '              September 1, 1990', /,
     + 15X,        '                Alice R. Drew', /,
     + 15X,        '          Warren T. Dewhurst, Ph.D.', /,
     + 15X,        '         Lieutenant Commander, NOAA', /)

      WRITE (LUOUT,931)
  931 FORMAT (15X,  '          (Hit RETURN to continue.)')

      READ (LUIN,'(A1)') ANS
      WRITE (LUOUT,2)
*   2 FORMAT ('1')
    2 FORMAT ('')

      WRITE (LUOUT,932)
  932 FORMAT ( /, 32X, 'DISCLAIMER' ,//,
     + ' This program and supporting information is furnished by',
     + ' the government of', /,
     + ' the United States of America, and is accepted/used by the',
     + ' recipient with', /,
     + ' the understanding that the U. S. government makes no',
     + ' warranties, express or', /,
     + ' implied, concerning the accuracy, completeness, reliability,',
     + ' or suitability', /,
     + ' of this program, of its constituent parts, or of any',
     + ' supporting data.', //,
     + ' The government of the United States of America shall be',
     + ' under no liability', /,
     + ' whatsoever resulting from any use of this program.',
     + '  This program should', /,
     + ' not be relied upon as the sole basis for solving a problem',
     + ' whose incorrect', /,
     + ' solution could result in injury to person or property.')
        WRITE (LUOUT,933)
  933   FORMAT ( /,
     + ' This program is the property of the government of the',
     + ' United States of', /,
     + ' America. Therefore, the recipient further agrees not to',
     + ' assert proprietary', /,
     + ' rights therein and not to represent this program to anyone as',
     + ' being other', /,
     + ' than a government program.', /)

      WRITE (LUOUT,931)
      READ (LUIN,'(A1)') ANS

      RETURN
      END
      SUBROUTINE IFILES (KIN, FLAG, BASEIN, RIDENT, PGM)

* Interactively get input file basename.
* Open the input file

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

      CHARACTER*32 B32
      PARAMETER (B32 = '                                ')

      REAL ANGLE
      INTEGER NZ, LRECL, IFLAG2, N2
      INTEGER KIN
      CHARACTER*56 RIDENT
      CHARACTER*32 FING
      CHARACTER*28 BASEIN
      CHARACTER*8 PGM
      CHARACTER*1 ANS
      LOGICAL FLAG

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NC, NR
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

      DATA IFLAG2 /2/

* Find out about the input file

      WRITE (LUOUT,2)
*   2 FORMAT ('1')
    2 FORMAT ('')
  220 WRITE (LUOUT,225)
  225 FORMAT (   ' For the input file enter:',
     +        /, '     ''A'' for ASCII transfer format;',
     +        /, '     ''B'' for binary - GEOID data file format);',
     +        /, ' (Default is B).')
      READ (LUIN,'(A1)') ANS
      IF (ANS .EQ. ' '  .OR.  ANS .EQ. 'B'  .OR.  ANS .EQ. 'b') THEN
        KIN = 1
      ELSEIF (ANS .EQ. 'A'  .OR.  ANS .EQ. 'a') THEN
        KIN = -1
      ELSE
        WRITE (LUOUT, 210) ANS
  210   FORMAT (/, ' ERROR - ''', A1, ''' is not a legal answer.')
        GOTO 220
      ENDIF

* Get basename

      IF (KIN .EQ. 1) THEN
        WRITE (LUOUT,110)
  110   FORMAT (   ' Enter the basename for the input GEOID',
     +             ' grid file from which the',
     +          /, ' new grid will be extracted.  The ''.geo''',
     +             ' extension will be',
     +          /, ' added to the basename by GEOGRD.  The default',
     +             ' full name is ''conus.geo''.')
      ELSEIF (KIN .EQ. -1) THEN
        WRITE (LUOUT,115)
  115   FORMAT (   ' Enter the basename for the input GEOID',
     +             ' grid file from which the',
     +          /, ' new grid will be extracted.  The ''.asc''',
     +             ' extension will be',
     +          /, ' added to the basename by GEOGRD.  The default',
     +             ' full name is ''conus.asc''.')
      ENDIF

      READ (LUIN,'(A28)') BASEIN
      CALL NBLANK (BASEIN, IFLAG2, N2)
      IF (N2 .EQ. 0) THEN
        BASEIN = 'conus'
        N2 = 5
      ENDIF

      FING = B32
      FING(1:N2) = BASEIN(1:N2)

      IF (KIN .EQ. 1) THEN

* Open binary input file

        N2 = N2 + 4
        FING(N2-3:N2) = '.geo'
        LRECL = 256
        CLOSE (NING)
        OPEN (NING,FILE=FING,FORM='UNFORMATTED',ACCESS='DIRECT',
     +        RECL=LRECL,STATUS='OLD',ERR= 900)
        READ (NING,REC=1) RIDENT, PGM, NC, NR, NZ, XMIN0, DX,
     +                    YMIN0, DY, ANGLE
      ELSEIF (KIN .EQ. -1) THEN

* Open ASCII input file

        N2 = N2 + 4
        FING(N2-3:N2) = '.asc'
        OPEN (NING,FILE=FING,FORM='FORMATTED',ACCESS='SEQUENTIAL',
     +        STATUS='OLD',ERR=900)
        READ (NING,140) RIDENT, PGM
  140   FORMAT (A56, A8)
        READ (NING, 150,ERR=900)  NC, NR, NZ, XMIN0, DX,
     +                            YMIN0, DY, ANGLE
  150   FORMAT (3I4, 5F12.5)
      ENDIF

* Check for corrupted file

      IF (NZ .NE. 1  .OR.  ANGLE .NE. 0.0) GOTO 960

* Reopen binary file with correct record length

      IF (KIN .EQ. 1) THEN
        CLOSE (NING)
        LRECL = 4*(NC + 1)
        OPEN (NING,FILE=FING,FORM='UNFORMATTED',ACCESS='DIRECT',
     +        RECL=LRECL,STATUS='OLD')
      ENDIF

* Write the (user-significant) input file variables to the screen

      XMAX0 = DX * REAL( NC-1 )  +  XMIN0
      YMAX0 = DY * REAL( NR-1 )  +  YMIN0

      IF (KIN .EQ. 1) THEN
        WRITE (LUOUT,130) FING(1:N2)
  130   FORMAT (/, ' Binary file ''', A, ''' has been opened.')
      ELSEIF (KIN .EQ. -1) THEN
        WRITE (LUOUT,135) FING(1:N2)
  135   FORMAT (/, ' ASCII file ''', A, ''' has been opened.')
      ENDIF
      WRITE (LUOUT,180) YMIN0, YMAX0, -XMAX0, -XMIN0
  180 FORMAT (/, ' The minimum and maximum latitudes  (+N) are: ',
     +           2F8.0,
     +        /, ' The minimum and maximum longitudes (+W) are: ',
     +           2F8.0, /)
      WRITE (LUOUT,200)
  200 FORMAT (14X, '(Hit RETURN to continue.)')
      READ (LUIN,'(A1)') ANS
      WRITE (LUOUT,2)

  999 RETURN

* Error messages

  900 WRITE (LUOUT,910) FING(1:N2)
  910 FORMAT (/, ' *** ERROR *** ''', A, ''' cannot be opened!',
     +        /, ' Do you wish to try again (Y/N)?',
     +        /, ' (Default is Y)')
      READ (LUIN,'(A1)') ANS
      IF (ANS .NE. 'N'  .AND.  ANS .NE. 'n') GOTO 220
      FLAG = .TRUE.
      GOTO 999

* Corrupted input file

  960 WRITE (LUOUT, 965) FING(1:N2)
  965 FORMAT (/, 5X, ' *********** ERROR ***********', /,
     +           ' The header information in grid file', /,
     +           ' ''', A, '''', /,
     +           ' is incorrect!  This file is corrupted or the wrong',
     +                                                       ' format.')
      FLAG = .TRUE.
      GOTO 999
      END
      SUBROUTINE NBLANK (A, IFLAG, NBLK)

*** Return position of last non-blank in string (IFLAG = 2)
*** or position of first non-blank in string (IFLAG = 1)

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

      INTEGER LENG, IFLAG, IBLK, NBLK
      CHARACTER*(*) A

      LENG = LEN(A)

      IF (IFLAG .EQ. 2) THEN
        DO 1 IBLK = LENG, 1, -1
          IF ( A(IBLK:IBLK) .NE. ' ' ) THEN
            NBLK = IBLK
            RETURN
          ENDIF
    1   CONTINUE
      ELSEIF (IFLAG .EQ. 1) THEN
        DO 2 IBLK = 1, LENG, +1
          IF ( A(IBLK:IBLK) .NE. ' ' ) THEN
            NBLK = IBLK
            RETURN
          ENDIF
    2   CONTINUE
      ENDIF

* String contains all blanks

      NBLK = 0

      RETURN
      END
      SUBROUTINE NLAT (KOUT, YMIN, YMAX)

* Get the new latitude variables

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* SMALL is system dependent and should be a very small real number

      REAL SMALL
      PARAMETER (SMALL = 1.E-5)

      REAL RCARD
      REAL YMIN, YMAX
      INTEGER LENG, IERR, KOUT
      CHARACTER*10 ANS

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NR, NC
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

  110 WRITE (LUOUT,115) YMIN0
  115 FORMAT (/, ' Enter the minimum latitude, north positive.',
     +        /, ' The default value is:', F8.0)
      READ (LUIN,'(A10)') ANS
      IF (ANS .EQ. '          ') THEN
        YMIN = YMIN0
      ELSE
        LENG = 10
        YMIN = RCARD (ANS, LENG, IERR)
      ENDIF

      WRITE (LUOUT,120) YMAX0
  120 FORMAT (/, ' Enter the maximum latitude, north positive.',
     +        /, ' The default value is:', F8.0)
      READ (LUIN,'(A10)') ANS
      IF (ANS .EQ. '          ') THEN
        YMAX = YMAX0
      ELSE
        LENG = 10
        YMAX = RCARD (ANS, LENG, IERR)
      ENDIF

* check values

      IF (YMIN .LT. YMIN0  .OR. YMIN .GT. YMAX0  .OR.
     +    YMAX .LT. YMIN0  .OR. YMAX .GT. YMAX0) THEN
        WRITE (LUOUT, 130) YMIN0, YMAX0
  130   FORMAT (/, ' *** ERROR *** ',
     +          /, ' One or both of the values that you gave for the',
     +             ' minimum and maximum',
     +          /, ' latitudes are out of bounds!  For the input',
     +             ' GEOID grid',
     +          /, ' the minimum and maximum latitudes (+N)  are: ',
     +             2F8.0)
        GOTO 110
      ELSEIF ( ( YMIN .GT. YMAX )  .OR.
     +         ( KOUT .EQ. -2  .AND.  YMIN .EQ. YMAX) ) THEN
        WRITE (LUOUT, 140) YMIN0, YMAX0
  140   FORMAT (/, ' *** ERROR *** ',
     +          /, ' The minimum value that you gave for the latitude',
     +             ' is less than the maximum!',
     +          /, ' For the input GEOID grid',
     +          /, ' the minimum and maximum latitudes (+N)  are: ',
     +             2F8.0)
        GOTO 110
      ENDIF

* The subgrid must be buffered by at least one whole degree of
* the original grid - BUT must be wholy within the original grid.

* For ASCII graphics format, buffer only to the nearest whole degree

      IF (KOUT .NE. -2) THEN
        YMIN = AINT(YMIN - 1.E0 + SMALL)
        YMAX = AINT(YMAX + 2.E0 - SMALL)
      ELSE
        YMIN = AINT(YMIN        + SMALL)
        YMAX = AINT(YMAX + 1.E0 - SMALL)
      ENDIF

      IF (YMIN .LT. YMIN0) YMIN = YMIN0
      IF (YMAX .GT. YMAX0) YMAX = YMAX0

      RETURN
      END
      SUBROUTINE NLON (KOUT, XMIN, XMAX)

* Get the new longitude variables.

* The variables LMIN is the minimum number of degrees between the
* minimum longitude and the maximum longitude.  The reason that
* there must be this minimum is that in order to be read by GEOID,
* the extracted grid must have a logical record length of at least 96.
* This is because the header record length is 96 bytes.
* This variables is not used for the graphics output (KOUT=-2)

* What happens is that the output file is opened with a record length
* of 4*(NC2+1) and the first record written.  The second record
* is then written into the file starting at character 4*(NC2+1)+1.
* where NC2 is the number of columns in the output file and is
* determined by the longitude difference and the longitude increment.
* Since the geoid heights are written as REAL*4 numbers, there must be
* at least 24 numbers in a row.

* For example, the conus file has DX=.08333.  GEOGRD requires that the
* grid boundarys be even degrees.  In order to have at least 24
* number in a row, the (maximum-minimum) longitude must be at least 6
* degrees.  However, given buffering, the inputted (maximum-minimum)
* longitude must be at least 4 degrees.

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* SMALL is system dependent and should be a very small real number

      REAL SMALL
      PARAMETER (SMALL = 1.E-5)

      REAL RCARD
      REAL XMIN, XMAX
      INTEGER LENG, IERR, KOUT, LMIN
      CHARACTER*10 ANS

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NR, NC
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

  160 WRITE (LUOUT,150) -XMAX0
  150 FORMAT (/, ' Enter the minimum longitude,',
     +           ' west longitude positive.',
     +        /, ' The default value is:', F8.0)
      READ (LUIN,'(A10)') ANS
      IF (ANS .EQ. '          ') THEN
        XMAX = XMAX0
      ELSE
        LENG = 10
        XMAX = RCARD (ANS, LENG, IERR)
        XMAX = -XMAX
      ENDIF
      IF (KOUT .NE. -2) THEN
        LMIN = INT(24.E0*DX + SMALL)
        IF (LMIN .LT. 1) LMIN = 1
        WRITE (LUOUT,155) LMIN, -XMIN0
  155   FORMAT (/, ' Enter the maximum longitude,',
     +             ' west longitude positive.',
     +          /, ' The (buffered) difference between the minimum and',
     +             ' maximum longitudes',
     +          /, ' must be at least', I3, ' degrees.  The default',
     +             ' value is:', F8.0)
      ELSE
        WRITE (LUOUT,165) -XMIN0
  165   FORMAT (/, ' Enter the maximum longitude,',
     +             ' west longitude positive.',
     +          /, ' The default value is:', F8.0)
      ENDIF
      READ (LUIN,'(A10)') ANS
      IF (ANS .EQ. '          ') THEN
        XMIN = XMIN0
      ELSE
        LENG = 10
        XMIN = RCARD (ANS, LENG, IERR)
        XMIN = -XMIN
      ENDIF

* check values

      IF (XMIN .LT. XMIN0  .OR. XMIN .GT. XMAX0  .OR.
     +    XMAX .LT. XMIN0  .OR. XMAX .GT. XMAX0) THEN
        WRITE (LUOUT, 170) -XMAX0, -XMIN0
  170   FORMAT (/, ' *** ERROR *** ',
     +          /, ' One or both of the values that you gave for the',
     +             ' minimum and maximum',
     +          /, ' longitudes are out of bounds!  For the input',
     +             ' GEOID grid',
     +          /, ' the minimum and maximum longitudes (+W)  are: ',
     +             2F8.0)
        GOTO 160
      ELSEIF ( ( XMIN .GT. XMAX )  .OR.
     +         ( KOUT .EQ. -2  .AND.  XMIN .EQ. XMAX) ) THEN
        WRITE (LUOUT, 180) -XMAX0, -XMIN0
  180   FORMAT (/, ' *** ERROR *** ',
     +          /, ' The minimum value that you gave for the longitude',
     +             ' is less than the maximum!',
     +          /, ' For the input GEOID grid',
     +          /, ' the minimum and maximum longitudes (+W)  are: ',
     +             2F8.0)
        GOTO 160
      ENDIF

* The subgrid must be buffered by at least one whole degree of
* the original grid - BUT must be wholy within the original grid.

* For ASCII graphics format, buffer only to the nearest whole degree

      IF (KOUT .NE. -2) THEN
        XMIN = AINT(XMIN - 2.E0 + SMALL)
        XMAX = AINT(XMAX + 1.E0 - SMALL)
      ELSE
        XMIN = AINT(XMIN - 1.E0 + SMALL)
        XMAX = AINT(XMAX        - SMALL)
      ENDIF

      IF (XMIN .LT. XMIN0) XMIN = XMIN0
      IF (XMAX .GT. XMAX0) XMAX = XMAX0

* check range

      IF ( KOUT .NE. -2  .AND.  (XMAX-XMIN) .LT. LMIN ) THEN
        WRITE (LUOUT, 190) LMIN, -XMAX0, -XMIN0
  190   FORMAT (/, ' *** ERROR *** ',
     +          /, ' The difference between the minimum and maximum',
     +             ' longitudes MUST be at',
     +          /, ' least', I3, ' degrees. ',
     +             ' For the input GEOID grid',
     +          /, ' the minimum and maximum longitudes (+W)  are: ',
     +             2F8.0)
        GOTO 160
      ENDIF

      RETURN
      END
      SUBROUTINE NPARMS (KOUT, XMIN, XMAX, YMIN, YMAX,
     +                   NRFRST, NRLAST, NCFRST, NCLAST)

* Get the variables for the new grid file.

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* SMALL is system dependent and should be a very small real number

      REAL SMALL
      PARAMETER (SMALL = 1.E-5)

      REAL XMIN, XMAX, YMIN, YMAX
      INTEGER NCFRST, NCLAST, NRFRST, NRLAST
      INTEGER KOUT
      CHARACTER*1 ANS

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NR, NC
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

* If only printing grid information then set defaults and return

      IF (KOUT .EQ. 0) THEN
        YMIN = YMIN0
        YMAX = YMAX0
        XMAX = XMAX0
        XMIN = XMIN0
        NRFRST = 1
        NRLAST = NR
        NCFRST = 1
        NCLAST = NC
        RETURN
      ENDIF

*** Get the area variables of the grid

      WRITE (LUOUT,2)
*   2 FORMAT ('1')
    2 FORMAT ('')
      IF (KOUT .NE. -2) THEN
        WRITE (LUOUT,100)
  100   FORMAT (   ' The extracted GEOID grid will be automatically',
     +             ' buffered to the next',
     +          /, ' whole degree plus one by GEOGRD.  Please enter',
     +             ' the latitude and longitude',
     +          /, ' extremes for your area of interest.  Enter the',
     +             ' values in decimal degrees.',
     +          /, ' If you enter blanks, the extremes from the input',
     +             ' file will be used.')
      ELSE
        WRITE (LUOUT,105)
  105   FORMAT (   ' The extracted GEOID grid will be automatically',
     +             ' rounded to the next',
     +          /, ' whole degree by GEOGRD.  Please enter the',
     +             ' latitude and longitude',
     +          /, ' extremes for your area of interest.  Enter the',
     +             ' values in decimal degrees.',
     +          /, ' If you enter blanks, the extremes from the input',
     +             ' file will be used.')
      ENDIF

* latitude

  200 CALL NLAT (KOUT, YMIN, YMAX)

* longitude, change to east positive

      CALL NLON (KOUT, XMIN, XMAX)

* Ask if the new maximums and minimums are OK

      WRITE (LUOUT,2)
      IF ( KOUT .NE. -2  .AND.
     +     ( XMIN .NE. XMIN0  .OR.  XMAX .NE. XMAX0  .OR.
     +       YMIN .NE. YMIN0  .OR.  YMAX .NE. YMAX0 )  ) THEN
        WRITE (LUOUT, 210) YMIN, YMAX, -XMAX, -XMIN
  210   FORMAT (   ' For the output grid;',
     +          /, ' the (buffered) minimum and maximum latitudes ',
     +             ' (+N) are: ', 2F8.0,
     +          /, ' the (buffered) minimum and maximum longitudes',
     +             ' (+W) are: ', 2F8.0)
      ELSE
        WRITE (LUOUT, 215) YMIN, YMAX, -XMAX, -XMIN
  215   FORMAT (   ' For the output grid;',
     +          /, ' the minimum and maximum latitudes ',
     +             ' (+N) are: ', 2F8.0,
     +          /, ' the minimum and maximum longitudes',
     +             ' (+W) are: ', 2F8.0)
      ENDIF

      WRITE (LUOUT,220)
  220 FORMAT (/, ' If these values are OK, hit RETURN to continue.',
     +           '  Enter any other character',
     +        /, ' to re-enter the minimum and maximum latitudes',
     +           ' and longitudes.')
      READ (LUIN,'(A1)') ANS
      IF (ANS .NE. ' ') GOTO 200

* Calculate the new origin for the subgrid area,
* calculate the number of columns and rows

      NRFRST = NINT( (YMIN - YMIN0)/DY ) + 1
      NRLAST = NINT( (YMAX - YMIN0)/DY ) + 1
      NCFRST = NINT( (XMIN - XMIN0)/DX ) + 1
      NCLAST = NINT( (XMAX - XMIN0)/DX ) + 1

      RETURN
      END
      SUBROUTINE OFILES (KIN, KOUT, FOUTG, BASEIN, RIDENT, PGM)

* Interactively get output file basename.

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

      CHARACTER*28 B28
      PARAMETER (B28 = '                            ')
      CHARACTER*32 B32
      PARAMETER (B32 = '                                ')

      INTEGER KIN, KOUT, IFLAG2, N2
      CHARACTER*56 RIDENT
      CHARACTER*32 FOUTG
      CHARACTER*28 BASEIN, BASOUT, DBNAME
      CHARACTER*8 PGM
      CHARACTER*1 ANS

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NC, NR
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

      DATA IFLAG2 /2/

* Find out about the output file

  220 WRITE (LUOUT,115)
  115 FORMAT (   ' For the output file enter:',
     +        /, '     ''A'' for ASCII transfer format;',
     +        /, '     ''B'' for binary - GEOID data file format);',
     +        /, '     ''G'' for ASCII graphics format -',
     +           ' There are five header information records,',
     +        /, '         and the record lengths',
     +           ' can be VERY large;'
     +        /, '     ''I'' for input file information written to',
     +           ' the output file;',
     +        /, ' (Default is B).')
      READ (LUIN,'(A1)') ANS
      IF (ANS .EQ. ' '  .OR.  ANS .EQ. 'B'  .OR.  ANS .EQ. 'b') THEN
        KOUT = 1
      ELSEIF (ANS .EQ. 'A'  .OR.  ANS .EQ. 'a') THEN
        KOUT = -1
      ELSEIF (ANS .EQ. 'G'  .OR.  ANS .EQ. 'g') THEN
        KOUT = -2
      ELSEIF (ANS .EQ. 'I'  .OR.  ANS .EQ. 'i') THEN
        KOUT = 0
      ELSE
        WRITE (LUOUT, 210) ANS
  210   FORMAT (/, ' ERROR - ''', A1, ''' is not a legal answer.')
        GOTO 220
      ENDIF

* default basename is input basename unless the same format

      IF (KOUT .NE. KIN) THEN
        DBNAME = BASEIN
        CALL NBLANK (DBNAME, IFLAG2, N2)
      ELSE
        DBNAME = 'geogrd'
        N2 = 6
      ENDIF

  180 IF (KOUT .EQ. 1) THEN

* The extension for the binary output files is .geo

        WRITE (LUOUT,160) DBNAME(1:N2)
  160   FORMAT (/, ' Enter the file name for the new',
     +             ' GEOID grid.  Enter only the',
     +          /, ' basename.  The ''.geo'' extension',
     +             ' will be added to the',
     +          /, ' basename by GEOGRD.  The default file name',
     +             ' is ''', A, '.geo''.')
        READ (LUIN,'(A28)') BASOUT
        IF (BASOUT .EQ. B28) BASOUT = DBNAME
        IF (KIN .EQ. 1  .AND.  BASOUT .EQ. BASEIN) THEN
          WRITE (LUOUT,165) BASOUT
  165     FORMAT (/, ' ***** ERROR ****',
     +            /, ' File ''', A, ''' is the name of the input file!',
     +            /, ' You must choose another name.')
          GOTO 180
        ENDIF

        CALL NBLANK (BASOUT, IFLAG2, N2)
        FOUTG = B32
        FOUTG(1:N2) = BASOUT(1:N2)
        FOUTG(N2+1:N2+4) = '.geo'

      ELSEIF (KOUT .EQ. -1) THEN

* The extension for the ASCII output files is .asc

        WRITE (LUOUT,190) DBNAME(1:N2)
  190   FORMAT (/, ' Enter the file name for the new',
     +             ' GEOID grid.  Enter only the',
     +          /, ' basename.  The ''.asc'' extension',
     +             ' will be added to the',
     +          /, ' basename by GEOGRD.  The default file name',
     +             ' is ''', A, '.asc''.')
        READ (LUIN,'(A28)') BASOUT
        IF (BASOUT .EQ. B28) BASOUT = DBNAME
        IF (KIN .EQ. -1  .AND.  BASOUT .EQ. BASEIN) THEN
          WRITE (LUOUT,185) BASOUT
  185     FORMAT (/, ' ***** ERROR ****',
     +            /, ' File ''', A, ''' is the name of the input file!',
     +            /, ' You must choose another name.')
          GOTO 180
        ENDIF

        CALL NBLANK (BASOUT, IFLAG2, N2)
        FOUTG = B32
        FOUTG(1:N2) = BASOUT(1:N2)
        FOUTG(N2+1:N2+4) = '.asc'

      ELSEIF (KOUT .EQ. -2) THEN

* The extensions for the graphics output file is .grd

        WRITE (LUOUT,170) DBNAME(1:N2)
  170   FORMAT (/, ' Enter the file name for the new',
     +             ' GEOID grid.  Enter only the',
     +          /, ' basename.  The ''.grd'' extension',
     +             ' will be added to the',
     +          /, ' basename by GEOGRD.  The default file name',
     +             ' is ''', A, '.grd''.')
        READ (LUIN,'(A28)') BASOUT
        IF (BASOUT .EQ. B28) BASOUT = DBNAME

        CALL NBLANK (BASOUT, IFLAG2, N2)
        FOUTG = B32
        FOUTG(1:N2) = BASOUT(1:N2)
        FOUTG(N2+1:N2+4) = '.grd'

      ELSEIF (KOUT .EQ. 0) THEN

* The extension for the information file is .inf

        WRITE (LUOUT,195) DBNAME(1:N2)
  195   FORMAT (/, ' Enter the file name for the grid information',
     +             ' file.  Enter only the basename.',
     +          /, ' The ''.inf'' extension will be added to the',
     +             ' basename by GEOGRD.  The default',
     +          /, ' base file name is ''', A, '''.')
        READ (LUIN,'(A28)') BASOUT
        IF (BASOUT .EQ. B28) BASOUT = DBNAME
        CALL NBLANK (BASOUT, IFLAG2, N2)
        FOUTG = B32
        FOUTG(1:N2) = BASOUT(1:N2)
        FOUTG(N2+1:N2+4) = '.inf'

* Write input header information to the information file

        OPEN (NOUTG,FILE=FOUTG,FORM='FORMATTED',ACCESS='SEQUENTIAL',
     +        STATUS='UNKNOWN')
        IF (KIN .EQ. 1) THEN
          CALL NBLANK (BASEIN, IFLAG2, N2)
          WRITE (NOUTG,230) BASEIN(1:N2)//'.geo'
  230     FORMAT (/, ' File ''', A, ''' is binary.')
        ELSEIF (KIN .EQ. -1) THEN
          CALL NBLANK (BASEIN, IFLAG2, N2)
          WRITE (NOUTG,235) BASEIN(1:N2)//'.asc'
  235     FORMAT (/, ' File ''', A, ''' is ASCII.')
        ENDIF
        WRITE (NOUTG,250) RIDENT, PGM
  250   FORMAT (/, ' From the header record(s):',
     +          /, ' Data Identification =''', A56, '''',
     +          /, ' Originating software =''', A8, '''')
        WRITE (NOUTG,260) NC, NR
  260   FORMAT (/, ' Number of columns =', I5,
     +             '    Number of rows =', I5)
        WRITE (NOUTG,270) YMIN0, YMAX0, DY
  270   FORMAT (/, ' Latitude:   minimum =', F8.0,
     +             '   maximum =', F8.0,
     +             '  increment =', F9.5)
        WRITE (NOUTG,280) -XMAX0, -XMIN0, DX
  280   FORMAT (   ' Longitude:  minimum =', F8.0,
     +             '   maximum =', F8.0,
     +             '  increment =', F9.5)
      ENDIF

      RETURN
      END
      SUBROUTINE TONEW (KIN, KOUT, XMIN, XMAX, YMIN, YMAX,
     +                  NCFRST, NCLAST, NRFRST, NRLAST, FOUTG)

*** Copy from input file to output file

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* NMAX is the maximum size of the buffer.  It is at least as large
* as the largest number of columns in the input grid

      INTEGER NMAX
      PARAMETER (NMAX = 1024)

* Variables ending in an 'G' are associated with geoid height files.

      REAL XMIN, XMAX, YMIN, YMAX
      REAL ANGLE, SCR
      REAL BUFG(NMAX)
      INTEGER NCFRST, NCLAST, NRFRST, NRLAST
      INTEGER NZ, LRECL, KIN, KOUT, NRPR1, NRPR2
      INTEGER IFLAG2, N2
      INTEGER I, I2, IDUM, IREC1, IREC2, J, JJ, J1, J2, NR2, NC2
      CHARACTER*56 RIDENT
      CHARACTER*32 FOUTG
      CHARACTER*8 PGM
      CHARACTER*1 ASCR

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NC, NR
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

      DOUBLE PRECISION ZSUMG, ZSQRG
      REAL ZMING, ZMAXG
      COMMON /ZSTATS/ ZMING, ZMAXG, ZSUMG, ZSQRG

      DATA IFLAG2 /2/

* Open output file, write header information

      ANGLE = 0.0
      RIDENT = 'GEOID EXTRACTED REGION'
      PGM = 'GEOGRD'
      NZ = 1
      NR2 = NRLAST - NRFRST + 1
      NC2 = NCLAST - NCFRST + 1

      IF (KOUT .EQ. 1) THEN

* Binary output file

        LRECL = 4*(NC2 + 1)
        OPEN (NOUTG,FILE=FOUTG,FORM='UNFORMATTED',ACCESS='DIRECT',
     +        RECL=LRECL,STATUS='UNKNOWN')
        WRITE (NOUTG,REC=1) RIDENT, PGM, NC2, NR2, NZ, XMIN, DX,
     +                      YMIN, DY, ANGLE

      ELSEIF (KOUT .EQ. -1) THEN

* ASCII transfer format output file

        OPEN (NOUTG,FILE=FOUTG,FORM='FORMATTED',ACCESS='SEQUENTIAL',
     +        STATUS='UNKNOWN')
        WRITE (NOUTG,140) RIDENT, PGM
  140   FORMAT (A56, A8)
        WRITE (NOUTG,150) NC2, NR2, NZ, XMIN, DX,
     +                    YMIN, DY, ANGLE
  150   FORMAT (3I4, 5F12.5)

      ELSEIF (KOUT .EQ. -2) THEN

* ASCII graphics format output file (SURFER header information)

        OPEN (NOUTG,FILE=FOUTG,FORM='FORMATTED',ACCESS='SEQUENTIAL',
     +        STATUS='UNKNOWN')
        WRITE (NOUTG,160)
  160   FORMAT ('DSAA')
        WRITE (NOUTG,165) NC2, NR2
  165   FORMAT (2I12)
        WRITE (NOUTG,170) XMIN, XMAX
  170   FORMAT (2F12.4)
        WRITE (NOUTG,170) YMIN, YMAX
        WRITE (NOUTG,170) ZMING, ZMAXG
      ENDIF

      CALL NBLANK (FOUTG, IFLAG2, N2)
      WRITE (LUOUT,80) FOUTG(1:N2)
   80 FORMAT (/, ' GEOGRD is now creating the ', A, ' file.',
     +        /, ' This process may take several minutes.')

* number of ASCII rows per binary row (note that there are always
* an odd number of grid columns)

        IF (KIN .EQ. -1) THEN
          NRPR1 = NC/6 + 1

* For ASCII input file, skip up to the NRFRST equivalent row

          REWIND NING
          READ (NING, '(A1)') ASCR
          READ (NING, '(A1)') ASCR
          IF (NRFRST .GT. 1) THEN
            I2 = (NRFRST - 1)*NRPR1
            DO 210 I = 1, I2
              READ (NING, 200) SCR
  200         FORMAT (F12.6)
  210       CONTINUE
          ENDIF

        ENDIF

        IF (KOUT .EQ. -1) THEN
          NRPR2 = (NCLAST - NCFRST + 1)/6 + 1
        ELSEIF (KOUT .EQ. 1) THEN
          IREC2 = 1
        ENDIF

* Copy input file to output file

      DO 300 I = NRFRST, NRLAST

* read input file

        IF (KIN .EQ. 1) THEN
          IREC1 = I + 1
          READ (NING,REC=IREC1) IDUM, (BUFG(J), J = 1, NC)
        ELSE
          DO 320 JJ = 1, NRPR1-1
            J1 = (JJ-1)*6 + 1
            J2 = J1 + 5
            READ (NING, 310) (BUFG(J), J = J1, J2)
  310       FORMAT (6F12.6)
  320     CONTINUE
          J1 = J2 + 1
          READ (NING, 310) (BUFG(J), J = J1, NC)
        ENDIF

* write output file

        IF (KOUT .EQ. 1) THEN
          IDUM = 0
          IREC2 = IREC2 + 1
          WRITE (NOUTG,REC=IREC2) IDUM, (BUFG(J), J = NCFRST, NCLAST)
        ELSEIF (KOUT .EQ. -1) THEN
          DO 420 JJ = 1, NRPR2-1
            J1 = (JJ-1)*6 + NCFRST
            J2 = J1 + 5
            WRITE (NOUTG, 310) (BUFG(J), J = J1, J2)
  420     CONTINUE
          J1 = J2 + 1
          WRITE (NOUTG, 310) (BUFG(J), J = J1, NCLAST)
        ELSEIF (KOUT .EQ. -2) THEN
          WRITE (NOUTG,440) (BUFG(J), J = NCFRST, NCLAST)
  440     FORMAT (1024F12.6)
        ENDIF

  300 CONTINUE

      IF (KOUT .EQ. 1) THEN
        CALL NBLANK (FOUTG, IFLAG2, N2)
        WRITE (LUOUT,130) FOUTG(1:N2)
  130   FORMAT (/, ' Binary file ''', A, '''',
     +             ' has been created.')
      ELSEIF (KOUT .EQ. -1  .OR.  KOUT .EQ. -2) THEN
        CALL NBLANK (FOUTG, IFLAG2, N2)
        WRITE (LUOUT,135) FOUTG(1:N2)
  135   FORMAT (/, ' ASCII file ''', A, '''',
     +             ' has been created.')
      ENDIF

      RETURN
      END
      SUBROUTINE ZEXTRM (KIN, KOUT, NCFRST, NCLAST, NRFRST, NRLAST)

* Read from the input file to obtain the minimum and maximum Z
* values for the new grid, and other statistical values

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* NMAX is the maximum size of the buffer.  It is at least as large
* as the largest number of columns in the input grid

      INTEGER NMAX
      PARAMETER (NMAX = 1024)

* Variables ending in an 'G' are associated with the geoid height
* grid file.

      DOUBLE PRECISION AVEG, SIGG
      REAL BUFG(NMAX)
      REAL SCR
      INTEGER I2, NRPR1, JJ, J1, J2, NTOT
      INTEGER IREC, I, J, IDUM, KIN, KOUT, ILAST
      INTEGER NCFRST, NCLAST, NRFRST, NRLAST
      CHARACTER*1 ASCR

      REAL XMIN0, XMAX0, YMIN0, YMAX0, DX, DY
      INTEGER NC, NR
      COMMON /GRID0/ XMIN0, XMAX0, YMIN0, YMAX0, DX, DY, NR, NC

      INTEGER LUIN, LUOUT, NING, NOUTG
      COMMON /INOUT/ LUIN, LUOUT, NING, NOUTG

      DOUBLE PRECISION ZSUMG, ZSQRG
      REAL ZMING, ZMAXG
      COMMON /ZSTATS/ ZMING, ZMAXG, ZSUMG, ZSQRG

* Initialize

      ZMING = 9.9E10
      ZMAXG = -9.9E10
      ZSUMG = 0.D0
      ZSQRG = 0.D0

      NTOT = 0

      IF (KOUT .NE. 0) THEN
        WRITE (LUOUT,80)
   80   FORMAT (/, ' GEOGRD is now finding the maximum and minimum',
     +             ' geoid height values in the',
     +          /, ' new grid.  This process may take several minutes.')
      ELSE
        WRITE (LUOUT,85)
   85   FORMAT (/, ' GEOGRD is now finding the maximum and minimum',
     +             ' geoid height values in the',
     +          /, ' input grid.  This process may take several',
     +             ' minutes.')
      ENDIF

      IF (KIN .EQ. 1) THEN
        DO 300 I = NRFRST, NRLAST
          IREC = I + 1

          READ (NING,REC=IREC) IDUM, (BUFG(J), J = 1, NC)

          CALL ZSUMS (NCFRST, NCLAST, BUFG)
          NTOT = NTOT + NCLAST - NCFRST + 1

  300   CONTINUE

      ELSEIF (KIN .EQ. -1) THEN

        NRPR1 = NC/6 + 1

* For ASCII input file, skip up to the NRFRST equivalent row

        REWIND NING
        READ (NING, '(A1)') ASCR
        READ (NING, '(A1)') ASCR
        IF (NRFRST .GT. 1) THEN
          I2 = (NRFRST - 1)*NRPR1
          DO 210 I = 1, I2
            READ (NING, 200) SCR
  200       FORMAT (F12.6)
  210     CONTINUE
        ENDIF

        ILAST = NRLAST - NRFRST + 1
        DO 340 I = 1, ILAST
          DO 320 JJ = 1, NRPR1-1
            J1 = (JJ-1)*6 + 1
            J2 = J1 + 5
            READ (NING, 310) (BUFG(J), J = J1, J2)
  310       FORMAT (6F12.6)
  320     CONTINUE
          J1 = J2 + 1
          READ (NING, 310) (BUFG(J), J = J1, NC)
          CALL ZSUMS (NCFRST, NCLAST, BUFG)
          NTOT = NTOT + NCLAST - NCFRST + 1
  340   CONTINUE

      ENDIF

      AVEG = ZSUMG / NTOT
      SIGG = DSQRT(  ( DBLE(NTOT)*ZSQRG - ZSUMG*ZSUMG )/
     +       ( DBLE(NTOT)*DBLE(NTOT-1) )  )

      WRITE (LUOUT, 500)
  500 FORMAT (/, 7X, 'GEOID HEIGHT STATISTICS',
     +        /, 1X, 36('-'),
     +        /, 25X, 'Geoid height',
     +        /, 27X, '(meters)')
      WRITE (LUOUT, 510) ZMING
  510 FORMAT (' Minimum', 17X, F9.3)
      WRITE (LUOUT, 520) ZMAXG
  520 FORMAT (' Maximum', 17X, F9.3)
      WRITE (LUOUT, 530) AVEG
  530 FORMAT (/, ' Average', 17X, F9.3)
      WRITE (LUOUT, 540) SIGG
  540 FORMAT (' Standard Deviation', 6X, F9.3)

      IF (KOUT .EQ. 0) THEN
        WRITE (NOUTG, 500)
        WRITE (NOUTG, 510) ZMING
        WRITE (NOUTG, 520) ZMAXG
        WRITE (NOUTG, 530) AVEG
        WRITE (NOUTG, 540) SIGG
        WRITE (LUOUT,*)
      ENDIF

      RETURN
      END
      SUBROUTINE ZSUMS (NCFRST, NCLAST, BUFG)

* Do the summations for the subroutine ZEXTRM

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

* NMAX is the maximum size of the buffer.  It is at least as large
* as the largest number of columns in the input grid

      INTEGER NMAX
      PARAMETER (NMAX = 1024)

* Variables ending in an 'G' are associated with geoid height files.

      REAL BUFG(NMAX)
      REAL ZG
      INTEGER K
      INTEGER NCFRST, NCLAST

      DOUBLE PRECISION ZSUMG, ZSQRG
      REAL ZMING, ZMAXG
      COMMON /ZSTATS/ ZMING, ZMAXG, ZSUMG, ZSQRG

      DO 400 K = NCFRST, NCLAST

        ZG = BUFG(K)
        ZMING = MIN(ZMING, ZG)
        ZMAXG = MAX(ZMAXG, ZG)
        ZSUMG = ZSUMG + DBLE(ZG)
        ZSQRG = ZSQRG + DBLE(ZG*ZG)

  400 CONTINUE

      RETURN
      END
      REAL FUNCTION RCARD (CHLINE, LENG, IERR)

*** Read a real number from a line of card image.
*** LENG is the length of the card
*** blanks are the delimiters of the REAL*4 variable

*     IMPLICIT REAL (A-H, O-Z)
*     IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

      INTEGER LENG, IERR, I, J, ILENG
      REAL VAR
      CHARACTER*80 CHLINE

      IERR = 0

* Find first non-blank character

* DO WHILE line character is blank, I is first non-blank character

      I = 1
   10 IF ( CHLINE(I:I) .EQ. ' '  .OR.  CHLINE(I:I) .EQ. ',' ) THEN
        I = I + 1

* Check for totally blank card

        IF ( I .GE. LENG) THEN
          RCARD = 0.0E0
          LENG = 0
          RETURN
        ENDIF

      GOTO 10
      ENDIF

* Find first blank character (or end of line)

* DO WHILE line character is not a blank

      J = I + 1
   20 IF ( CHLINE(J:J) .NE. ' '  .AND.  CHLINE(J:J) .NE. ',' ) THEN
        J = J + 1

* Check for totally filed card

        IF ( J .GT. LENG) THEN
          GOTO 40
        ENDIF

      GOTO 20
      ENDIF

* J is now 1 more than the position of the last non-blank character

   40 J = J - 1

* ILENG is the length of the real string, it cannot be greater
* than 15 characters

      ILENG = J - I + 1

      IF (ILENG .GT. 15) THEN
        STOP 'RCARD'
      ENDIF

* Read the real variable from the line, and set the return VAR to it

      READ (CHLINE(I:J), 55, ERR=9999) VAR
   55 FORMAT (F15.0)
      RCARD = VAR

* Now reset the values of LENG and CHLINE to the rest of the card

      CHLINE( 1 : LENG ) = CHLINE( (J+1) : LENG )
      LENG = LENG - J

      RETURN

* Read error

 9999 IERR = 1
      RETURN
      END