      PROGRAM GEOID

***********************************************************************
*                                                                     *
* PROGRAM :   GEOID                                                   *
*                                                                     *
* PURPOSE:    COMPUTATION PROGRAM TO INTERPOLOATE GEOID HEIGHT        *
*             VALUES FROM A GRID OF REGULARLY-SPACED ESTIMATES.       *
*             THE INPUT REQUIRES LATITUDE AND LONGITUDE VALUES        *
*             REFERRENCED TO THE NORTH AMERICAN DATUM OF 1983         *
*             (NAD 83) AND THE WORLD GEODETIC SYSTEM OF 1980 (WGS 80).*
*             INPUT DATA REFERRENCED TO THE NORTH AMERICAN DATUM OF   *
*             1927 (NAD27) SHOULD BE CONVERTED TO NAD 83 VIA NADCON   *
*             PRIOR TO USING THIS PROGRAM.                            *
*                                                                     *
*             THE ACTUAL COMPUTATION IS PERFORMED AS AN INTERPOLATION *
*             FROM A REGULARLY-SPACED GRID OF POINTS OBTAINED FROM    *
*             TWO DIFFERENT SOURCES.  ONE SOURCE IS THE GEOID93 MODEL *
*             COMPUTED BY DR. DENNIS G. MILBERT, NGS, NOAA.  THE OTHER*
*             IS AN EVALUATION OF A GEOPOTENTIAL MODEL BASED ON A     *
*             COMPLETE EXPANSION OF SPHERICAL HARMONIC COEFFICIENTS   *
*             OBTAINED FROM OHIO STATE UNIVERSITY AND KNOWN AS OSU89B.*
*             THESE COEFFICIENTS WERE OBTAINED FROM DR. RICHARD RAPP  *
*             OF THE GEODETIC SCIENCE AND SURVEYING DEPARTMENT.       *
*                                                                     *
*             THE INTERPOLATION IS ACCOMPLISHED BY A LOCALLY FIT      *
*             BIQUADRATIC FUNCTION.                                   *
*                                                                     *
*             THE POLYNOMIAL SURFACE IS FIT TO NINE DATA POINTS       *
*             DEFINING THE AREA SURROUNDING THE (X,Y) PAIR            *
*             WHERE THE INTERPOLATION IS TO TAKE PLACE.               *
*                                                                     *
*             THE PROGRAM REQUIRES THAT THE USER SPECIFY:             *
*                                                                     *
*             1)  THE NAME OF AN OUTPUT FILE                          *
*                                                                     *
*             2)  THE NAME OF AN INPUT FILE (IF AVAILABLE).           *
*                                                                     *
*                                                                     *
*             THIS PROGRAM ALLOWS FOR EITHER NGS STANDARD HORIZONTAL  *
*             DATA FORMATS AS SPECIFIED IN THE FGCC PUBLICATION,      *
*             COMMONLY KNOWN AS THE 'HORIZONTAL BLUE BOOK' (SEE       *
*             SUBROUTINES TYPE3, TYPE4, TYPE5), OR IN A GENERIC FILE  *
*             FORMAT (SEE SUBROUTINE TYPE1 OR SUBROUTINE TYPE2).      *
*                                                                     *
*             THE CODE CAN BE EASILY MODIFIED TO ACCOMMODATE CUSTOM   *
*             FILE SPECIFICATIONS BY MODIFYING SUBROUTINES: ENDREP,   *
*             GETPT, IPARMS, WRTPT, AND (OPTIONALLY) FHELP.           *
*                                                                     *
*                                                                     *
* VERSION CODE:  1.00   (20-dec-90)                                   *
* VERSION CODE:  1.01   (27-feb-91)                                   *
*   increaded MXAREA to 10                                            *
* VERSION CODE:  2.00   (10-mar-93)                                   *
*                                                                     *
* VERSION DATE:  march 10, 1993                                       *
*                                                                     *
*                DENNIS MILBERT, Ph.D.        (GEOID93 model)         *
*                DONALD SCHULTZ               (GEOID93 model)         *
*                ADVANCED GEODETIC SCIENCE BRANCH                     *
*                NATIONAL GEODETIC SURVEY, NOS, NOAA                  *
*                ROCKVILLE, MD   20852                                *
***********************************************************************

***********************************************************************
*                                                                     *
*                  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)

      REAL VRSION
      INTEGER MXAREA
      PARAMETER (VRSION = 2.00E0, MXAREA =10)

      DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG
      DOUBLE PRECISION AGHT, VGHT
      DOUBLE PRECISION SDGHT, SDGHT2
      REAL GHT, SMGHT, BGGHT
      INTEGER IPAGE, ITYPE, NCONV, IFILE
      CHARACTER*80 NAME
      LOGICAL PAGE, NODATA, SCREEN

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

**********************
* INITIALIZE VARIABLES
**********************

      CALL INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE,
     +            SMGHT, BGGHT, AGHT, VGHT, SDGHT, SDGHT2,
     +            XSMALL, XBIG, YSMALL, YBIG)

**************************
* PRINT HEADER INFORMATION
**************************

      CALL HEADR (VRSION)

***********************
* OPEN GEOID DATA FILES
***********************

      CALL NGRIDS (NODATA)
      IF (NODATA) GOTO 9999

******************************************************
* REQUEST FOR THE NEEDED VARIABLE VALUES FROM THE USER
******************************************************

      CALL IPARMS (ITYPE, SCREEN)

********************************
* LOOP ONCE FOR EACH COMPUTATION
********************************

      CALL MLOOP (NCONV, IPAGE, ITYPE, VRSION, GHT,
     +            SDGHT, SDGHT2, SMGHT, BGGHT,
     +            XSMALL, XBIG, YSMALL, YBIG,
     +            PAGE, SCREEN, NAME)

***********************************************
* FINISHED WITH ALL COMPUTATIONS - WRITE REPORT
***********************************************

      CALL ENDREP (IPAGE, NCONV, ITYPE, GHT, AGHT, VGHT,
     +             SDGHT, SDGHT2, SMGHT, BGGHT,
     +             XSMALL, XBIG, YSMALL, YBIG)

*****************
* CLOSE ALL FILES
*****************

      DO 1010 IFILE = 1, NAREA
        CLOSE ( LUAREA(IFILE) )
 1010 CONTINUE
      CLOSE (NIN, STATUS='KEEP')
      CLOSE (NOUT, STATUS='KEEP')
      CLOSE (NAPAR, STATUS='KEEP')

 9999 STOP
      END

      SUBROUTINE ASKPT (NCONV, NAME, IDLA, IMLA, SLA,
     +                  IDLO, IMLO, SLO, XPT, YPT, EOF, NOPT)

* Interactively ask for the name and location of a point

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*40 B40
      PARAMETER (B40 = '                                        ')

      DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD
      REAL SLA, SLO, RMLA, RMLO, RCARD
      INTEGER NCONV
      INTEGER IDLA, IMLA, IDLO, IMLO
      INTEGER IFLAG1, IFLAG2, N1, N2
      INTEGER LENG, IERR, IOS
      CHARACTER*80 NAME
      CHARACTER*40 ANS, DUM
      LOGICAL EOF, NOPT

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA IFLAG1 /1/, IFLAG2 /2/

      NAME = '    '
      WRITE (LUOUT,*) ' What is the NAME for this station or',
     +                ' point?'
      READ (LUIN,'(A80)') NAME

**********
* LATITUDE
**********

      IF (NCONV .EQ. 0) THEN
        WRITE (LUOUT,110)
  110   FORMAT (/, ' Latitudes and Longitudes may be entered',
     +          ' in three formats:', /,
     +          '   (1) degrees, minutes, and decimal seconds, OR', /,
     +          '   (2) degrees, and decimal minutes OR', /,
     +          '   (3) decimal degrees.', /,
     +          ' Degrees, minutes and seconds may be separated',
     +          ' by blanks or commas.', /,
     +          ' A latitude or longitude of 0 will end data',
     +          ' entry.', /)
      ENDIF

        WRITE (LUOUT,*) ' What is its NAD 83 latitude?'
        WRITE (LUOUT,*) ' '

      READ (LUIN,170,ERR=9930,IOSTAT=IOS) ANS
  170 FORMAT (A40)
      IF (ANS .EQ. B40) GOTO 9999

      DUM = ANS
      CALL NBLANK (DUM, IFLAG2, N2)
      LENG = N2
      RDLA = DCARD( DUM(1:N2), LENG, IERR )
      IF (IERR .NE. 0) GOTO 9950
      IF (LENG .GT. 0) THEN

        RMLA = RCARD( DUM, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9950

        IF (LENG .GT. 0) THEN
          SLA  = RCARD( DUM, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9950
        ELSE
          SLA = 0.E0
        ENDIF

      ELSE
        RMLA = 0.E0
        SLA = 0.E0
      ENDIF

      IF ( (RDLA .EQ. 0.D0)  .AND.  (RMLA .EQ. 0.E0)  .AND.
     +     (SLA .EQ. 0.E0) ) GOTO 9999

* Check for illogical values

      IF (RDLA .LT.   0.D0) GOTO 9940
      IF (RDLA .GT.  90.D0) GOTO 9950
      IF (RMLA .LT. 0.E0  .OR.  RMLA .GT. 60.E0) GOTO 9950
      IF ( SLA .LT. 0.E0  .OR.   SLA .GT. 60.E0) GOTO 9950

***********
* LONGITUDE
***********

      WRITE (LUOUT,*) ' What is its NAD 83 longitude?',
     +                '  (Longitude is positive west.)'
      WRITE (LUOUT,*) ' '

      READ (LUIN,170,ERR=9930,IOSTAT=IOS) ANS
      IF (ANS .EQ. B40) GOTO 9999

      DUM = ANS
      CALL NBLANK (DUM, IFLAG2, N2)
      LENG = N2
      RDLO = DCARD( DUM(1:N2), LENG, IERR )
      IF (IERR .NE. 0) GOTO 9960
      IF (LENG .GT. 0) THEN

        RMLO = RCARD( DUM, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9960

        IF (LENG .GT. 0) THEN
          SLO  = RCARD( DUM, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9960
        ELSE
          SLO = 0.E0
        ENDIF

      ELSE
        RMLO = 0.E0
        SLO = 0.E0
      ENDIF

      IF ( (RDLO .EQ. 0.D0)  .AND.  (RMLO .EQ. 0.E0)  .AND.
     +     (SLO .EQ. 0.E0) ) GOTO 9999

* Check for illogical values

      IF (RDLO .LT.   0.D0) GOTO 9940
      IF (RDLO .GT. 360.D0) GOTO 9960
      IF (RMLO .LT. 0.E0  .OR.  RMLO .GT. 60.E0) GOTO 9960
      IF ( SLO .LT. 0.E0  .OR.   SLO .GT. 60.E0) GOTO 9960

* Calculate decimal degrees

      YPT = RDLA + DBLE(RMLA)/60.D0 + DBLE(SLA)/3600.D0
      XPT = RDLO + DBLE(RMLO)/60.D0 + DBLE(SLO)/3600.D0

* Get degrees, minutes, seconds

      CALL HMS (YPT, IDLA, IMLA, SLA)
      CALL HMS (XPT, IDLO, IMLO, SLO)

 9000 RETURN

* Error messages

 9930 CONTINUE
      CALL NBLANK (ANS, IFLAG1, N1)
      CALL NBLANK (ANS, IFLAG2, N2)
      WRITE (LUOUT,9935) ANS(N1:N2)
 9935 FORMAT (' ERROR - in the answer:', /,
     +        9X, '''', A, '''', /,
     +        '         Must enter number in prescribed format!', /)
      NOPT = .TRUE.
      GOTO 9000

 9940 WRITE (LUOUT,9945)
      CALL NBLANK (ANS, IFLAG1, N1)
      CALL NBLANK (ANS, IFLAG2, N2)
      WRITE (LUOUT,9945) ANS(N1:N2)
 9945 FORMAT (' ERROR - in the answer:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude and Longitudes must be positive!', /,
     +        '         Longitude is positive west.', /)
      NOPT = .TRUE.
      GOTO 9000

 9950 CONTINUE
      CALL NBLANK (ANS, IFLAG1, N1)
      CALL NBLANK (ANS, IFLAG2, N2)
      WRITE (LUOUT,9955) ANS(N1:N2)
 9955 FORMAT (' ERROR - Illogical value for latitude in the answer:', /,
     +        '         ''', A, '''', /,
     +        '         Latitude must be between 0 and 90 degrees.', /,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9960 CONTINUE
      CALL NBLANK (ANS, IFLAG1, N1)
      CALL NBLANK (ANS, IFLAG2, N2)
      WRITE (LUOUT,9965) ANS(N1:N2)
 9965 FORMAT (' ERROR - Illogical value for longitude in the answer:',/,
     +        '         ''', A, '''', /,
     +        '         Longitude must be between 0 and 360 degrees.',/,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9999 CONTINUE
      EOF = .TRUE.
      GOTO 9000
      END

      SUBROUTINE DGRIDS

* This subroutine opens the geoid grid(s) using the default grid
* names and locations.  The default names of the grid areas are
* given in DAREAS and the default base file locations are in DFILES

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

      INTEGER MXAREA, MXDEF
      PARAMETER (MXAREA =10, MXDEF = MXAREA)
      CHARACTER*80 B80
      CHARACTER*20 B20
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)
      CHARACTER*65 B65
      PARAMETER (B65 = B20//B20//B20//'     ')

      REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1
      INTEGER IDEF, ITEMP, NR1, NC1
      CHARACTER*80 DUM
      CHARACTER*65 AFILE, DFILES(MXAREA)
      CHARACTER*15 DAREAS(MXDEF)
      LOGICAL NOGO, GFLAG

      CHARACTER*15 AREAS
      COMMON /AREAS/ AREAS(MXAREA)

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      DATA DUM / B80 /

* DFILES contains the default locations (pathname) of the grid files
* without the .geo extension. (For example 'conus' would
* indicate that the conus.geo grid file is in the
* current working directory.)  The length of each entry in DFILES may
* be up to 65 characters.  DAREAS contains the default names of these
* areas.  The names are used internally in the program and in the
* program output.  They may be no longer than 15 characters.  They
* must correspond on a one-to-one basis with the file locations in
* the DFILES array.  That is, the first area name in DAREAS must
* be the name that you wish for the first data file set in the
* DFILES array.  You may, of course, have the arrays the same if
* the location of the data file is no longer than 15 characters.
* The locations of the grid files may be differ for each
* installation.  If the pathnames are not correct DFILES (and, possibly,
* DAREAS) may be changed and the program recompiled.

      DATA DFILES /'conus', 'hawaii', 'prvi', 'alaska',
     +              ' ', ' ', ' ', ' ', ' ', ' '/
      DATA DAREAS /'Conus OSU89B', 'Hawaii OSU89B',
     +             'PR & VI OSU89B', 'Alaska OSU89B',
     +             ' ', ' ', ' ', ' ', ' ', ' '/

      GFLAG = .FALSE.
      WRITE (LUOUT, 80)
   80 FORMAT (/, '      Default Data Grids', /,
     +           '   #  AREA NAME', /, 1X, 79('=') )

      DO 140 IDEF = 1, MXDEF
        AFILE = DFILES(IDEF)
        IF (AFILE .EQ. B65) GOTO 999

* Try to open a set of default files.
* Do not print error messages for non-existing files.

        ITEMP = NAREA + 1
        CALL OPENFL (AFILE, ITEMP, NOGO, DX1, DY1,
     +               XMAX1, XMIN1, YMAX1, YMIN1, NR1, NC1, DUM, GFLAG)

        IF (.NOT. NOGO) THEN

* Set of files opened OK and variables read

          NAREA = ITEMP
          AREAS(NAREA) = DAREAS(IDEF)
          DX(NAREA) = DX1
          DY(NAREA) = DY1
          XMAX(NAREA) = XMAX1
          XMIN(NAREA) = XMIN1
          YMAX(NAREA) = YMAX1
          YMIN(NAREA) = YMIN1
          NR(NAREA) = NR1
          NC(NAREA) = NC1

          WRITE (LUOUT,120) NAREA, AREAS(NAREA)
  120     FORMAT (2X, I2, 2X, A15)

        ENDIF

  140 CONTINUE

  999 RETURN
      END

      SUBROUTINE DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG)

* This subroutine prints out a small diagram showing the area
* that was computed.

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

      DOUBLE PRECISION XTEMP, XSMALL, XBIG, YSMALL, YBIG
      REAL SLOMIN, SLOMAX, SLAMIN, SLAMAX
      INTEGER LODMIN, LOMMIN, LODMAX, LOMMAX
      INTEGER LADMIN, LAMMIN, LADMAX, LAMMAX
      INTEGER LU, NCONV

      XTEMP = XSMALL
      XSMALL = XBIG
      XBIG = XTEMP
      CALL HMS (XBIG,   LODMIN, LOMMIN, SLOMIN)
      CALL HMS (XSMALL, LODMAX, LOMMAX, SLOMAX)
      CALL HMS (YSMALL, LADMIN, LAMMIN, SLAMIN)
      CALL HMS (YBIG,   LADMAX, LAMMAX, SLAMAX)

      WRITE (LU,90) NCONV
   90 FORMAT (//, ' The total number of geoid estimates: ', I8)
      WRITE (LU,100)
  100 FORMAT (//, 30X, 'Region of Estimation')

      WRITE (LU,1000) LODMAX, LOMMAX, SLOMAX, LODMIN, LOMMIN, SLOMIN,
     +                LADMAX, LAMMAX, SLAMAX, LADMAX, LAMMAX, SLAMAX,
     +                LADMIN, LAMMIN, SLAMIN, LADMIN, LAMMIN, SLAMIN,
     +                LODMAX, LOMMAX, SLOMAX, LODMIN, LOMMIN, SLOMIN

 1000 FORMAT (5(/), T4, 'NAD 83', /,
     +       T4, 'Longitude:', 8X, I4, 1X, I2.2, 1X, F6.3, 13X, I4, 1X,
     +       I2.2, 1X, F6.3, /,
     +       T4, 'Latitude:', 9X, I4, 1X, I2.2, 1X, F6.3,
     +       ' ************', I4, 1X, I2.2, 1X, F6.3, /,
     +       5(T27, 3X, '*', T57, '*', /),
     +       T27, 3X, '*', 9X, ' NAD 83', T57, '*', /,
     +       T27, 3X, '*', 9X, '  data ', T57, '*', /,
     +       T27, 3X, '*', 9X, ' points', T57, '*', /,
     +       5(T27, 3X, '*', T57, '*', /),
     +       T4, 'Latitude:', 9X, I4, 1X, I2.2, 1X, F6.3,
     +       ' ************', I4, 1X, I2.2, 1X, F6.3, /,
     +       T4, 'Longitude:', 8X, I4, 1X, I2.2, 1X, F6.3, 13X, I4, 1X,
     +       I2.2, 1X, F6.3, //)

      RETURN
      END

      SUBROUTINE ENDREP (IPAGE, NCONV, ITYPE, GHT, AGHT, VGHT,
     +                   SDGHT, SDGHT2, SMGHT, BGGHT,
     +                   XSMALL, XBIG, YSMALL, YBIG)

*** Gather statistics and write the end-of-program report

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG
      DOUBLE PRECISION AGHT, VGHT
      DOUBLE PRECISION SDGHT, SDGHT2
      REAL GHT, SMGHT, BGGHT
      INTEGER IPAGE, NCONV, ITYPE
      INTEGER LU, I
      LOGICAL PAGE

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      PAGE = .TRUE.
      IPAGE = IPAGE + 1

*******************
* DO THE STATISTICS
*******************

      IF (NCONV .GE. 2) THEN

* calculate mean, variance, standard deviation for geoid height
* estimates (meters).

*******************
* GEOID HT (METERS)
*******************

        AGHT = SDGHT/DBLE(NCONV)
        VGHT = SDGHT2/DBLE(NCONV-1) - SDGHT**2/DBLE( NCONV*(NCONV-1 ))

        IF (VGHT .GT. 1.0D-6) THEN
          SDGHT = DSQRT(VGHT)
        ELSE
          VGHT = 0.0D0
          SDGHT = 0.0D0
        ENDIF
      ENDIF

      IF (NCONV .LT. 2) THEN
        AGHT = GHT
        VGHT = 0.0D0
        SDGHT = 0.0D0
      ENDIF

***************************************
* PRINT OUT THE STATISTICS FOR THIS JOB
***************************************

      LU = LUOUT
      IF (NCONV .GT. 0) THEN

*************************************************
* FIRST REPORT THE FINAL STATISTICS TO THE SCREEN
*************************************************

        CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT, IPAGE, PAGE)

        CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG)

****************************************************
* NOW REPORT THE FINAL STATISTICS TO THE OUTPUT FILE
****************************************************

        IF (ITYPE .EQ. 0) THEN

**************************************
* INTERACTIVE USE ONLY - NO INPUT FILE
**************************************

          LU = NOUT

          CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT,
     +                 IPAGE, PAGE)

          CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG)

          IF (NCONV .EQ. 0) THEN
            DO 1007 I = 1, 2
              WRITE (LU,*) ' All of your NAD 83 LOCATIONS are out of',
     +                     ' bounds.'
              LU = NOUT

 1007       CONTINUE
          ENDIF

        ELSEIF (ITYPE .EQ. 1) THEN

* For file format types 1 and 2

          LU = NOUT
          CALL REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT,
     +                 IPAGE, PAGE)

          CALL DIAGRM (LU, NCONV, XSMALL, XBIG, YSMALL, YBIG)

        ELSEIF (ITYPE .EQ. 2) THEN

* ITYPE = 2, (free format input, free format output) does not have
* anything written to the output file in this subroutine

        ELSEIF (ITYPE .EQ. 3  .OR.  ITYPE .EQ. 4  .OR.
     +          ITYPE .EQ. 5 .or.  ITYPE .EQ. 6) THEN
      
* ITYPE=3, ITYPE=4, and ITYPE=5 & 6 the NGS Horizontal Blue Book
* file formats do not have a * report written to the output file

        ENDIF

      ENDIF

      RETURN
      END

      SUBROUTINE FGRID (XPT, YPT, DX, DY, XMAX, XMIN,
     +                  YMAX, YMIN, XGRID, YGRID, IROW, JCOL, NOGO)

**********************************************************************
** SUBROUTINE FGRID: IDENTIFIES THE LOCAL GRID SQUARE FOR INTERP.    *
**********************************************************************

* This subroutine is designed to identify the grid square in which a
* particular point is located and get the corner coordinates
* converted into the index coordinate system.

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

      DOUBLE PRECISION XPT, YPT
      REAL DX, DY, XMAX, XMIN, YMAX, YMIN, XGRID, YGRID
      INTEGER IROW, JCOL
      LOGICAL NOGO

      NOGO = .FALSE.

* Check to see it the point is outside the area of the gridded data

      IF (XPT .GE. DBLE(XMAX)  .OR.  XPT .LE. DBLE(XMIN)   .OR.
     +    YPT .GE. DBLE(YMAX)  .OR.  YPT .LE. DBLE(YMIN) ) THEN
        NOGO = .TRUE.
*       WRITE (*,*) '***THE POINT IS OUT OF BOUNDS***'
        GOTO 200
      ENDIF

* Calculate the coordinate values for the point to be interpolated
* in terms of grid indices

      XGRID = SNGL(XPT - DBLE(XMIN) )/DX + 1.E0
      YGRID = SNGL(YPT - DBLE(YMIN) )/DY + 1.E0

* Find the I,J values for the SW corner of the local square

      IROW = IFIX(YGRID)
      JCOL = IFIX(XGRID)

  200 RETURN
      END

      SUBROUTINE FHELP

*** Print information about the formats of the input data
*** file types used by GEOID.

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

      INTEGER MXAREA, ITYPE, IOS
      PARAMETER (MXAREA =10)

      CHARACTER*1 ANS

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

*****************************
* CHOSE THE INPUT FILE FORMAT
*****************************

 9001 WRITE (LUOUT,*) ' '
      WRITE (LUOUT,*) ' What format do you want information about?'
      WRITE (LUOUT,*) '  1) Free Format Type 1'
      WRITE (LUOUT,*) '  2) Free Format Type 2'
      WRITE (LUOUT,*) '  3) NGS Blue Book (old) format (default)'
      WRITE (LUOUT,*) '  4) NGS Blue Book (new *84*) format'
      WRITE (LUOUT,*) '  5) NGS Blue Book (new *86*) format'
      WRITE (LUOUT,*) '  6) NGS Blue Book (extended) format'

      READ (LUIN,'(A1)') ANS
      IF (ANS .EQ. ' ') THEN

        ITYPE = 3

      ELSE
        READ (ANS,90,ERR=9940,IOSTAT=IOS) ITYPE
   90   FORMAT (I1)

        IF (ITYPE .GT. 6  .OR.  ITYPE .LT. 1) THEN

* Not a good answer - Try again.

          WRITE (LUOUT,*) ' Gotta pick one of these -',
     +                    ' sorry try again.'
          GOTO 9001
        ENDIF
      ENDIF

* Print information

      IF (ITYPE .EQ. 1) THEN

*******************************
* FOR FILE FORMAT ITYPE = 1
* FREE FORMAT TYPE1 INPUT FILE
* PRETTY OUTPUT FORMAT
*******************************

        WRITE (LUOUT,2)
*   2   FORMAT ('1')
    2   FORMAT ('')
        WRITE (LUOUT, 110)
  110   FORMAT (' Free Format Type 1 - The first 40 characters of',
     +                                  ' the input data record may', /,
     +          ' contain the station name or be blank.  The rest of',
     +                                 ' the record (columns 41-80)', /,
     +          ' must contain the latitude and longitude.  They may',
     +                                    ' be given in (1) decimal', /,
     +          ' degrees; (2) integer degrees and decimal minutes,',
     +                                    ' or (3) integer degrees,', /,
     +          ' integer minutes, and decimal seconds.  The decimal',
     +                                    ' portion of the latitude', /,
     +          ' MUST contain a decimal point as it is used to',
     +                                ' determine which is the last', /,
     +          ' number forming part of the latitude.  The output',
     +                                ' will be in "pretty" format.', /)
        WRITE (LUOUT, 120)
  120   FORMAT (' The following three records are examples of valid',
     +                                            ' input records:', /,
     +          ' <------------ Columns 1-40 ------------>',
     +                     '<------------ Columns 41-80----------->', /,
     +          ' AAA                                     34.',
     +                                     '4444444      98.8888888', /,
     +          ' BBB                                     25',
     +                                   ' 55.55555     76 56.66666', /,
     +          ' CCC                                     45 45',
     +                                     ' 45.555   111 11 11.111', /)

        WRITE (LUOUT,931)
  931   FORMAT (15X,    '     (Hit RETURN to continue.)')
        READ (LUIN,'(A1)') ANS
        WRITE (LUOUT,2)

        WRITE (LUOUT, 130)
  130   FORMAT (' The following is an example of the output.  Note',
     +                               ' that with Free Format Type 1', /,
     +          ' data, the input latitude and',
     +                                 ' longitude are expressed in', /,
     +          ' degrees, minutes, and seconds regardless of the',
     +                                           ' method of input.', /)
        WRITE (LUOUT, 140)
  140   FORMAT ('                           Computation #:    1',
     +                               '        Region: Conus OSU89B', //,
     +          '  Station Name:  AAA', //,
     +          '                      Latitude                 ',
     +                               'Longitude       Geoid Height', /,
     +          '                                                  ',
     +                                    '               (meters)', //,
     +          '                   34 26 39.99984           98 53',
     +                                     ' 19.99968       -26.868', /)

      ELSEIF (ITYPE .EQ. 2) THEN

*******************************
* FOR FILE FORMAT ITYPE = 2
* FREE FORMAT TYPE2 INPUT FILE
* FREE FORMAT OUTPUT
*******************************

        WRITE (LUOUT,2)
        WRITE (LUOUT, 210)
  210   FORMAT (' Free Format Type 2 - The first 32 characters of',
     +                                 ' the input data record must', /,
     +          ' contain the latitude and longitude.  They may be',
     +                              ' given in (1) decimal degrees;', /,
     +          ' (2) integer degrees and decimal minutes, or (3)',
     +                                   ' integer degrees, integer', /,
     +          ' minutes, and decimal seconds.  The decimal portion',
     +                               ' of the latitude MUST contain', /,
     +          ' a decimal point as it is used to determine which',
     +                                 ' is the last number forming', /,
     +          ' part of the latitude.  The next eight characters',
     +                              ' (columns 33-40) will have the', /,
     +          ' calculated geoid height while the rest of the input',
     +                                 ' record (columns 41-80) may', /,
     +          ' contain the station name or be blank.  The output',
     +                              ' will be in the same format as', /,
     +          ' the input but will contain the geoid height.', /)

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

        WRITE (LUOUT, 220)
  220   FORMAT (' The following three records are examples of',
     +                                      ' valid input records:', //,
     +          ' <------------ Columns 1-40 ------------>',
     +                     '<------------ Columns 41-80----------->', /,
     +          ' 45 45 45.55555 111 11 11.11111          ',
     +                                                         'one', /,
     +          ' 25 55.5555555   76 56.6666666           ',
     +                                                         'two', /,
     +          ' 34.444444444    98.888888888            ',
     +                                                          'three')

        WRITE (LUOUT, 230)
  230   FORMAT (/, ' The following is an example of the output.', //,
     +             ' GEOID      Version 2.00 - NGS',
     +                                              ' GEOID93 Model', /,
     +             '   45 45 45.55555 111 11 11.11111 -10.859one', /,
     +             '   25 55.5555555   76 56.6666666  -39.988two', /,
     +             '   34.444444444    98.888888888   -26.868three', /)

      ELSEIF (ITYPE .EQ. 3) THEN

****************************************
* FOR INPUT FILE ITYPE = 3
* THE HORIZONTAL BLUE BOOK SPECIFICATION
****************************************

* The following two format statements logically belong together
* but are split in order to have fewer that 19 continuation lines

        WRITE (LUOUT,2)
        WRITE (LUOUT, 310)
  310   FORMAT (' NGS Horizontal ''Old'' Blue Book format - *80*',
     +                            ' (Control Point) and *84* (Geoid', /,
     +          ' Height) Records.  GEOID uses information from the',
     +                              ' *80* and *84* records in Blue', /,
     +          ' Book files.  The Station Serial Number (SSN), the',
     +                                 ' latitude and the longitude', /,
     +          ' are read from the *80* records, the rest of the',
     +                               ' record is ignored.  Only the', /,
     +          ' *84* record is modified by GEOID - the *80*',
     +                          ' records and all other records are')
        WRITE (LUOUT, 311)
  311   FORMAT (' passed through without change to the output file.',
     +                              '  On the *84* records, the SSN', /,
     +          ' is used to determine which is the corresponding',
     +                                      ' the *80* record.  The', /,
     +          ' information to the right of the SSN (columns',
     +                               ' 15-80) in the *84* record is', /,
     +          ' created by GEOID.  On the *80* records, the',
     +                             ' direction of the latitude must', /,
     +          ' be north positive (''N'') and the direction of the',
     +                                     ' longitude must be west', /,
     +          ' positive (''W'').', /)

        WRITE (LUOUT, 320)
  320   FORMAT (' This format is out of date and should no longer be',
     +                         ' used for submitting data', /,
     +          ' to the National Geodetic Survey.', //,
     +          ' This format was defined in:', /,
     +          '   ''Input Formats and Specifications of the',
     +                       ' National Geodetic Survey Data Base''', /,
     +          '   ''Volume 1. Horizontal Control Data''.', /,
     +          ' Published by the Federal Geodetic Control Committee',
     +                                      ' in 1980 and available', /,
     +          ' from: the National Geodetic Survey, NOAA, Rockville,',
     +                                                   'MD 20852.', /)

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

        WRITE (LUOUT, 330)
  330   FORMAT (' The following input example is *80* and *84* records',
     +                                  ' from an ''old'' Blue Book', /,
     +          ' file:', //,
     +          '004560*80*096 KNOXVILLE CT HSE',
     +          '              411906578  N0930548534  W 277  MIA33', /,
     +          '004565*84*096 RAPP - OSU86G                        ',
     +                                  '                     -319  10')

        WRITE (LUOUT, 340)
  340   FORMAT (/, ' The following example is of the output *84*',
     +                          ' record with the new, interpolated', /,
     +          ' geoid height.', //,
     +          '004565*84*096 NGS   PROGRAM GEOID      Conus OSU89B',
     +                               '                     -316   5', /)

      ELSEIF (ITYPE .EQ. 4) THEN

* The following two format statements logically belong together
* but are split in order to have fewer that 19 continuation lines

        WRITE (LUOUT,2)
        WRITE (LUOUT,410)
  410   FORMAT (' NGS Horizontal ''New'' Blue Book format - *80*',
     +                            ' (Control Point) and *84* (Geoid', /,
     +          ' Height) Records.  GEOID uses information from the',
     +                              ' *80* and *84* records in Blue', /,
     +          ' Book files.  The Station Serial Number (SSN), the',
     +                                 ' latitude and the longitude', /,
     +          ' are read from the *80* records, the rest of the',
     +                               ' record is ignored.  Only the', /,
     +          ' *84* record is modified by GEOID - the *80*',
     +                          ' records and all other records are')
        WRITE (LUOUT, 411)
  411   FORMAT (' passed through without change to the output file.',
     +                              '  On the *84* records, the SSN', /,
     +          ' is used to determine which is the corresponding',
     +                                      ' the *80* record.  The', /,
     +          ' information to the right of the SSN (columns',
     +                               ' 15-80) in the *84* record is', /,
     +          ' created by GEOID.  On the *80* records, the',
     +                             ' direction of the latitude must', /,
     +          ' be north positive (''N'') and the direction of the',
     +                                     ' longitude must be west', /,
     +          ' positive (''W'').', /)

        WRITE (LUOUT, 420)
  420   FORMAT (' For more information on this format,',
     +                                           ' please refer to:', /,
     +          '   ''Input Formats and Specifications of the',
     +                       ' National Geodetic Survey Data Base''', /,
     +          '   ''Volume 1. Horizontal Control Data''.', /,
     +          ' Published by the Federal Geodetic Control Committee',
     +                                       ' in January, 1989 and', /,
     +          ' available from: the National Geodetic Survey, NOAA,',
     +                                       ' Rockville, MD 20852.', /)

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

        WRITE (LUOUT, 430)
  430   FORMAT (' The following input example is *80* and *84* records',
     +                                   ' from a ''new'' Blue Book', /,
     +          ' file:', //,
     +          '004560*80*0096KNOXVILLE CT HSE              ',
     +                        '411906578  N0930548534  W 277  MIA33', /,
     +          '004565*84*0096RAPP                          ',
     +                           '                          -3190  100')

        WRITE (LUOUT, 440)
  440   FORMAT (/, ' The following example is of the output *84*',
     +                          ' record with the new, interpolated', /,
     +          ' geoid height.', //,
     +          '004565*84*0096NGS   PROGRAM GEOID      Conus OSU89B',
     +                               '                   -3164   50', /)

*********************************************************************
***   *86* records
*********************************************************************
      ELSEIF (ITYPE .EQ. 5) THEN

* The following two format statements logically belong together
* but are split in order to have fewer that 19 continuation lines

        WRITE (LUOUT,2)
        WRITE (LUOUT,510)
  510   FORMAT (' NGS Horizontal ''New'' Blue Book format - *80*',
     +                            ' (Control Point) and *86* (Combo', /,
     +          ' Height) Records.  GEOID uses information from the',
     +                              ' *80* and *86* records in Blue', /,
     +          ' Book files.  The Station Serial Number (SSN), the',
     +                                 ' latitude and the longitude', /,
     +          ' are read from the *80* records, the rest of the',
     +                               ' record is ignored.  Only the', /,
     +          ' *86* record is modified by GEOID - the *80*',
     +                          ' records and all other records are')
        WRITE (LUOUT, 511)
  511   FORMAT (' passed through without change to the output file.',
     +                              '  Any *84* records are ', /,
     +          ' deleted.  Information in columns',
     +                               ' 36-43) in the *86* record is', /,
     +          ' created by GEOID.  On the *80* records, the',
     +                             ' direction of the latitude must', /,
     +          ' be north positive (''N'') and the direction of the',
     +                                     ' longitude must be west', /,
     +          ' positive (''W'').', /)

        WRITE (LUOUT, 520)
  520   FORMAT (' For more information on this format,',
     +                                           ' please refer to:', /,
     +          '   ''Input Formats and Specifications of the',
     +                       ' National Geodetic Survey Data Base''', /,
     +          '   ''Volume 1. Horizontal Control Data''.', /,
     +          ' Published by the Federal Geodetic Control Committee',
     +                                       ' in January, 1989 and', /,
     +          ' available from: the National Geodetic Survey, NOAA,',
     +                                       ' Rockville, MD 20852.', /)

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

        WRITE (LUOUT, 530)
  530   FORMAT (' The following input example is *80* and *86* records',
     +                                   ' from a ''new'' Blue Book', /,
     +          ' file:', //,
     +          '004560*80*0096KNOXVILLE CT HSE              ',
     +                        '411906578  N0930548534  W 277  MIA33', /,
     +          '004565*86*0096 any values in here           ',
     +                        'will be reproduced if ssn matches   ')

        WRITE (LUOUT, 540)
  540   FORMAT (/, ' The following example is of the output *86*',
     +                          ' record with the new, interpolated', /,
     +          ' geoid height.', //,
     +          '004565*86*0096 any values in here   -31768D ',
     +                        'will be reproduced if ssn matches   ',/)

*********************************************************************
*********************************************************************
      ELSEIF (ITYPE .EQ. 6) THEN

        WRITE (LUOUT,2)
        WRITE (LUOUT,610)
  610   FORMAT (' NGS Horizontal ''Extended'' Blue Book format',
     +                         ' (BIGADJUST format) - *80* (Control', /,
     +          ' Point) and *84* (Geoid Height) Records.  Only the',
     +                                 ' *80* records are used from', /,
     +          ' the input file.  The output consists only of *84*'
     +                               ' records.  The SSN''s for the', /,
     +          ' *84* records are taken from the *80* records - that',
     +                                       ' is, there is an *84*', /,
     +          ' record for each *80* record, regardless of any',
     +                           ' *84* records in the input file.', //,
     +          ' This format is for internal NGS use only!!!', /)

      ENDIF

      WRITE (LUOUT,*) ' Do you want more information (Y/N)?'
      WRITE (LUOUT,*) ' (Default is Y)'
      READ (LUIN,'(A1)') ANS
      IF (ANS .NE. 'N'  .AND.  ANS .NE. 'n') GOTO 9001

      RETURN

* Error message

 9940 WRITE (LUOUT,*) ' Gotta pick ''1'' through ''6'' -',
     +                ' sorry try again.'
      GOTO 9001
      END

      SUBROUTINE GETPT (NCONV, ITYPE, NAME, IDLA, IMLA, SLA, IDLO,
     +                  IMLO, SLO, XPT, YPT, EOF, NOPT)

* Get the name, latitude, and longitude of a point either interactively
* or from an input data file

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      DOUBLE PRECISION XPT, YPT
      REAL SLA, SLO
      INTEGER NCONV, ITYPE
      INTEGER IDLA, IMLA, IDLO, IMLO
      CHARACTER*80 NAME
      CHARACTER*1 ANS
      LOGICAL EOF, NOPT

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      EOF = .FALSE.
      NOPT = .FALSE.

      IF (ITYPE .EQ. 0) THEN

*************************************
* FOR INTERACTIVE USE - NO INPUT FILE
*************************************

        IF (NCONV .GE. 1) THEN
          WRITE (LUOUT,*) ' Do you want to do another geoid',
     +                    ' computation (Y/N)?'
          WRITE (LUOUT,*) ' (Default is Y)'
          READ (LUIN,'(A1)') ANS
          IF (ANS .EQ. 'n'  .OR.  ANS .EQ. 'N') GOTO 9999
        ENDIF

* Get a point (X,Y) to compute

        CALL ASKPT (NCONV, NAME, IDLA, IMLA, SLA,
     +              IDLO, IMLO, SLO, XPT, YPT, EOF, NOPT)
        IF (NOPT) GOTO 9000

      ELSEIF (ITYPE .EQ. 1) THEN

* Free format type 1

        CALL TYPE1 (NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +              XPT, YPT, EOF, NOPT)

      ELSEIF (ITYPE .EQ. 2) THEN

* Free format type 2

        CALL TYPE2 (IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +              XPT, YPT, EOF, NOPT)

      ELSEIF (ITYPE .EQ. 3  .OR.  ITYPE .EQ. 4  .OR.
     +        ITYPE .EQ. 5  .or. ITYPE .EQ. 6 ) THEN

* ITYPE=3 is NGS old Horizontal Blue Book
* ITYPE=4 is NGS new Horizontal Blue Book with *84*
* ITYPE=5 is NGS new Horizontal Blue Book with *86*
* ITYPE=6 is NGS extended Horizontal Blue Book

        CALL TYPE3 (ITYPE, NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +              XPT, YPT, EOF, NOPT)

      ENDIF

*********************************************************
* CHANGE THE LONGITUDE TO POSITIVE EAST FOR INTERPOLATION
*********************************************************

      XPT = -XPT

 9000 RETURN

* End of file

 9999 CONTINUE
      EOF = .TRUE.
      GOTO 9000
      END

      SUBROUTINE GRIDS

* This subroutine opens the geoid grid(s) that are requested in
* a file named 'AREA.PAR' (if it exists).
* AREA.PAR (or area.par) will be read for the names
* and locations of the gridded geoid height data set files.

* The data in the AREA.PAR file has the following format:
* Columns 1-15 contain the name of the area.  This name may
*   contain blanks or any other characters.
* Columns 16-80 (the rest of the record) contain the base name of the
*   grid files.  That is the '.geo' extension is not
*   included - it is added by this subroutine before each file is
*   opened.

* Comments are indicated by an '*' in column 1, blank records are
* also treated as comments.  Comment records are printed to the
* output file but otherwise ignored.

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*80 B80
      CHARACTER*20 B20
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1
      INTEGER NR1, NC1
      INTEGER IOS, I, IFLAG1, IFLAG2, N1, N2, LENG, IERR, ITEMP
      CHARACTER*80 CARD, CCARD, DUM
      CHARACTER*65 AFILE, GFILE
      CHARACTER*15 AAREA
      LOGICAL NOGO, GFLAG

      CHARACTER*15 AREAS
      COMMON /AREAS/ AREAS(MXAREA)

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      DATA DUM / B80 /
      DATA IFLAG1 /1/, IFLAG2 /2/

***********************************************************************
* TRY AND OPEN THE 'AREA.PAR' FILE, IF NOT THEN TRY AND OPEN 'area.par'
***********************************************************************

      GFILE = 'AREA.PAR'
      OPEN (NAPAR,FILE=GFILE,FORM='FORMATTED',STATUS='OLD',
     +      ACCESS='SEQUENTIAL',ERR=9100,IOSTAT=IOS)

* File containing grid names exits

   10 GFLAG = .TRUE.
      WRITE (LUOUT, 80)
   80 FORMAT (/, '      Data Grids named in the AREA.PAR file', /,
     +           '   #  AREA NAME      LOCATION', /, 1X, 79('=') )

      DO 120 I = 1, MXAREA
  100   READ (NAPAR,110,ERR=9000,END=9000,IOSTAT=IOS) CARD
  110   FORMAT (A80)

* Check for comment records and blank records

        IF ( CARD(1:1) .EQ. '*' ) THEN
          CALL NBLANK (CARD, IFLAG2, N2)
          WRITE (LUOUT,'(5X, A)') CARD(1:N2)
          GOTO 100
        ELSEIF ( CARD .EQ. B80 ) THEN
          WRITE (LUOUT,*) ' '
          GOTO 100
        ENDIF

* Get area name and basename of file (i.e. location without extensions)

        DUM = CARD
        AAREA = DUM(1:15)
        CALL NBLANK (CARD(16:80), IFLAG1, N1)
        DUM(1:65) = CARD(15+N1:80)
        LENG = 65
        AFILE = CCARD(DUM, LENG, IERR)
        IF (IERR .NE. 0) STOP 'Coding error in GRIDS -- AFILE'
        IF (AFILE .EQ. '        ') GOTO 9000

* Try and open files

        ITEMP = NAREA + 1
        CALL OPENFL (AFILE, ITEMP, NOGO, DX1, DY1,
     +               XMAX1, XMIN1, YMAX1, YMIN1, NR1, NC1, CARD, GFLAG)

        IF (.NOT. NOGO) THEN

* Files opened OK and variables read

          NAREA = ITEMP
          AREAS(NAREA) = AAREA
          DX(NAREA) = DX1
          DY(NAREA) = DY1
          XMAX(NAREA) = XMAX1
          XMIN(NAREA) = XMIN1
          YMAX(NAREA) = YMAX1
          YMIN(NAREA) = YMIN1
            NR(NAREA) = NR1
            NC(NAREA) = NC1

          CALL NBLANK (CARD, IFLAG2, N2)
          WRITE (LUOUT,140) NAREA, CARD(1:N2)
  140     FORMAT (2X, I2, 2X, A)
        ENDIF

  120 CONTINUE

 9000 RETURN

* 'AREA.PAR' does not exist, try the name 'area.par'
* If it exists, open it and continue, if it does not exist, return.

 9100 CONTINUE
      GFILE = 'area.par'
      OPEN (NAPAR,FILE=GFILE,FORM='FORMATTED',STATUS='OLD',
     +      ACCESS='SEQUENTIAL',ERR=9000,IOSTAT=IOS)
      GOTO 10
      END

      SUBROUTINE HEADR (VRSION)

*** This subroutine prints the header information and the disclaimer

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      REAL VRSION
      CHARACTER*1 ANS

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      WRITE (LUOUT,920)
  920 FORMAT (10X, '                   Welcome', /,
     +        10X, '                   to  the', /,
     +        10X, '           National Geodetic Survey', /,
     +        10X, '           Geoid Computation Program', //,
     +        10X, '               NGS GEOID93 Model', /)

      WRITE (LUOUT,930) VRSION

  930 FORMAT (//,10X, '               (Version', F5.2, ')', /,
     +       10X,     '               March 10, 1993', /,
     +       10X,     '            Dennis Milbert, Ph.D.',/,
     +       10X,     '               Donald Schultz',/,
     +       10X,     '       Advanced Geodetic Science Branch', /)

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

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

      WRITE (LUOUT,932)
  932 FORMAT (/, '                           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
      WRITE (LUOUT,2)

      RETURN
      END

      SUBROUTINE HMS (DD, ID, IM, S)

* Use this to change from decimal degrees (double precision)
* to integer degrees, integer minutes, and decimal seconds (real)
* Seconds are assumed to have no more than 5 decimal places

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

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

      DOUBLE PRECISION DD, TEMP
      REAL S
      INTEGER ID, IM

      ID = IDINT(DD)
      TEMP = ( DD - DBLE(ID) )*60.0D0
      IM = IDINT(TEMP)
      S = SNGL( ( TEMP - DBLE(IM) )*60.0D0 )

      IF (IM .EQ. 60) THEN
        IM = 0
        ID = ID + 1
      ENDIF

      IF (S .LT. SMALL) S = 0.E0

      IF (S .GT. (60.E0-SMALL)  ) THEN
        S = 0.E0
        IM = IM + 1
      ENDIF

      RETURN
      END

      SUBROUTINE INITL (SCREEN, PAGE, NAME, IPAGE, ITYPE,
     +                  SMGHT, BGGHT, AGHT, VGHT, SDGHT, SDGHT2,
     +                  XSMALL, XBIG, YSMALL, YBIG)

*** This subroutine initializes all the variables needed in GEOID

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*20 B20
      CHARACTER*80 B80
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG
      DOUBLE PRECISION AGHT, VGHT
      DOUBLE PRECISION SDGHT, SDGHT2
      REAL SMGHT, BGGHT
      INTEGER IPAGE, ITYPE
      CHARACTER*80 NAME
      LOGICAL PAGE, SCREEN

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

* Initialize card variable in common CURNT to blank

      CARD = B80

* Set the logical units for input/output common INOUT
* Note that the unit numbers for the data grids (LUAREA) are defined
* in subroutine OPENFL as the numbers 11 through MXAREA

      LUIN  = 5
      LUOUT = 6
      NOUT  = 101
      NIN   = 102
      NAPAR = 103

******************************************************************
*                             INITIALIZE
******************************************************************

* Defaults: SCREEN = .TRUE. => send results to screen
*           PAGE = .FALSE.  => don't start a new page in the output file
*           NAME = '     '  => no station name
*           IPAGE = 0       => current output file page number is 0
*           ITYPE = 0       => interactive input of points

      SCREEN = .TRUE.
      PAGE = .FALSE.
      NAME = '    '
      IPAGE = 0
      ITYPE = 0

      SMGHT =  1.0E10
      BGGHT = -1.0E10

      AGHT = 0.0D0
      VGHT = 0.0D0
      SDGHT = 0.0D0
      SDGHT2 = 0.0D0

      XSMALL =  1.0D10
      XBIG   = -1.0D10
      YSMALL =  1.0D10
      YBIG   = -1.0D10

      RETURN
      END

      SUBROUTINE INTERP (NOGO, RESP, XPT, YPT, GHT, ITYPE)

*** COMPUTE BIQUADRATIC INTERPOLATION, FLOATING POINT VERSION
*** LOGIC NOT TESTED FOR '0/360 MERIDIAN WRAPAROUND'

* THIS CODE IS ADAPTED FROM THAT OF DENNIS MILBERT - 6/5/90
* DENNIS(XPT,YPT,DX,DY,XMAX,XMIN,YMAX,YMIN,
*     &             ZEE,NR,NC,LUAREA)

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

      INTEGER MNCBIN, MXAREA
      PARAMETER (MNCBIN = 1000, MXAREA =10)
      CHARACTER*20 B20
      CHARACTER*80 B80
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      DOUBLE PRECISION XPT, YPT
      REAL BUF, GHT
      REAL X, Y, FX0, FX1, FX2, QTERP2
      REAL DX0, DY0, XMAX0, YMAX0, XMIN0, YMIN0
      REAL XGRID, YGRID, VAL
      REAL TEE1, TEE2, TEE3, TEE4, TEE5, TEE6, TEE7, TEE8, TEE9
      INTEGER NR0, NC0
      INTEGER IAREA, IROW, JCOL, IFILE, IDUM, I, J, IDIM
      INTEGER IX, IX1, IX2, JY, JY1, JY2
      INTEGER IFLAG1, IFLAG2, N1, N2, ITYPE
      CHARACTER*15 RESP
      LOGICAL NOGO, FLAG
      DIMENSION BUF(MNCBIN)

      CHARACTER*15 AREAS
      COMMON /AREAS/ AREAS(MXAREA)

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      DATA IFLAG1 /1/, IFLAG2 /2/, FLAG /.FALSE./

******************************************************************
*                             INITIALIZE
******************************************************************

      NOGO  =  .FALSE.

****************************************************
* READ WHERE TO GET THE DATA AND HOW IT IS ORGANIZED
****************************************************

* CHECK TO SEE WHICH SET OF GRIDDED FILES XPT,YPT is in.

      DO 100 IAREA = 1, NAREA

        DX0 = DX(IAREA)
        DY0 = DY(IAREA)
        XMAX0 = XMAX(IAREA)
        XMIN0 = XMIN(IAREA)
        YMAX0 = YMAX(IAREA)
        YMIN0 = YMIN(IAREA)
        NR0 = NR(IAREA)
        NC0 = NC(IAREA)

        CALL FGRID (XPT, YPT, DX0, DY0, XMAX0, XMIN0,
     +              YMAX0, YMIN0, XGRID, YGRID, IROW, JCOL, NOGO)
        IF (.NOT. NOGO) GOTO 200

  100 CONTINUE

* Not in any of the grid areas

      NOGO = .TRUE.
      GOTO 950

  200 CONTINUE

* Point in area number IAREA and named AREAS(IAREA)

      IFILE = LUAREA(IAREA)
      RESP = AREAS(IAREA)

      IF (YPT .LT. YMIN0  .OR.  YPT .GT. YMAX0  .OR.
     +    XPT .LT. XMIN0  .OR.  XPT .GT. XMAX0) THEN
        VAL=1.D30
      ELSE

*** Within grid boundaries, get indices in the grid

        IX  = (XPT - XMIN0)/DX0 + 1.D0
        IX1 = IX + 1
        IX2 = IX + 2

*** Check if edge collision

    1   IF (IX2 .GT. NC0) THEN
          IX  = IX -  1
          IX1 = IX1 - 1
          IX2 = IX2 - 1
          GO TO 1
        ENDIF

        X = (XPT - DX0*(IX - 1) - XMIN0)/DX0

*** move grid to get nearer center of 3 X 3

        IF (X .LT. 0.5  .AND.  IX .GT. 1) THEN
          IX  = IX  - 1
          IX1 = IX1 - 1
          IX2 = IX2 - 1
          X = X + 1.
        ENDIF
        IF ( X .LT. 0.  .OR.  X .GT. 2.) STOP 'INTERP - 1'

*** Now do it for Y

        JY = (YPT - YMIN0)/DY0 + 1.D0
        JY1 = JY + 1
        JY2 = JY + 2
    2   IF (JY2 .GT. NR0) THEN
          JY  = JY  - 1
          JY1 = JY1 - 1
          JY2 = JY2 - 1
          GO TO 2
        ENDIF
        Y = (YPT - DY0*(JY - 1) - YMIN0)/DY0
        IF (Y .LT. 0.5  .AND.  JY .GT. 1) THEN
          JY  = JY  - 1
          JY1 = JY1 - 1
          JY2 = JY2 - 1
          Y = Y + 1.
        ENDIF

        IF (Y .LT. 0.  .OR.  Y .GT. 2.) STOP 'INTERP - 2'

*** Get the nine data points needed from gridded geoid height file !!!

*** lower boundary

                READ (IFILE,REC=JY+1) IDUM, (BUF(J), J=1,NC0)
                TEE1 = BUF(IX)
                TEE2 = BUF(IX+1)
                TEE3 = BUF(IX+2)

*** middle row

                READ (IFILE,REC=JY1+1) IDUM, (BUF(J), J=1,NC0)
                TEE4 = BUF(IX)
                TEE5 = BUF(IX+1)
                TEE6 = BUF(IX+2)

*** upper boundary

                READ (IFILE,REC=JY2+1) IDIM, (BUF(J), J=1,NC0)
                TEE7 = BUF(IX)
                TEE8 = BUF(IX+1)
                TEE9 = BUF(IX+2)

*** Compute the value

        FX0 = QTERP2(X, TEE1, TEE2, TEE3)
        FX1 = QTERP2(X, TEE4, TEE5, TEE6)
        FX2 = QTERP2(X, TEE7, TEE8, TEE9)

        GHT = QTERP2(Y, FX0, FX1, FX2)
      ENDIF

9999  RETURN

* Error Messages

  950 CONTINUE
      IF (ITYPE .NE. 0) THEN
        CALL NBLANK (CARD, IFLAG1, N1)
        CALL NBLANK (CARD, IFLAG2, N2)
        WRITE (LUOUT,955) CARD(N1:N2)
  955   FORMAT (' *** THIS POINT IS OUT OF BOUNDS ***', /,
     +          1X, '''', A, '''')
      ELSE
        WRITE (LUOUT,960)
  960   FORMAT (' *** THE POINT IS OUT OF BOUNDS ***')
      ENDIF

* Write out grid areas for the first out-of-bounds error message

      IF (.NOT.FLAG  .OR.  ITYPE .EQ. 0) THEN
        WRITE (LUOUT,*) ' It must be within one of the following grid',
     +                  ' areas;'
        WRITE (LUOUT,975)
  975   FORMAT (18X, 7X, 'Latitude', 7X, 'Longitude', /,
     +          5X, 'Area Name', 5X, 2(5X, 'MIN', 4X, 'MAX', 1X),
     +          '(degrees)' )
        DO 970 I = 1, NAREA
          WRITE (LUOUT,965) AREAS(I), NINT( YMIN(I) ), NINT( YMAX(I) ),
     +                      NINT( -XMAX(I) ), NINT( -XMIN(I) )
  965     FORMAT (1X, '''', A15, '''', 2(2X, 2I7) )
  970   CONTINUE
        FLAG = .TRUE.
      ENDIF

      WRITE (LUOUT,*) ' '
      GOTO 9999

      END

      SUBROUTINE IPARMS (ITYPE, SCREEN)

*** This subroutine interactively requests for information
*** needed by GEOID.

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*80 B80
      CHARACTER*20 B20
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      CHARACTER*80 INFILE, OFILE
      CHARACTER*1 ANS
      INTEGER ITYPE
      INTEGER IFLAG1, IFLAG2, N1, N2, IOS
      LOGICAL SCREEN, EFLAG

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA IFLAG1 /1/, IFLAG2 /2/

**********************************
* GET THE NAME FOR THE OUTPUT FILE
**********************************

   14 WRITE (LUOUT,*) ' What is the name of the file which will',
     +                ' contain the program results?'
      WRITE (LUOUT,*) ' The default name is ''geoid.out''.'
      READ (LUIN,'(A80)') OFILE
      IF (OFILE .EQ. B80) OFILE = 'geoid.out'
      INQUIRE (FILE=OFILE, EXIST=EFLAG)
      IF (EFLAG) THEN
        CALL NBLANK (OFILE, IFLAG1, N1)
        CALL NBLANK (OFILE, IFLAG2, N2)
        WRITE (LUOUT,*) ' The file ''', OFILE(N1:N2), ''''
        WRITE (LUOUT,*) ' already exists.  Do you want to overwrite',
     +                                 ' it (Y/N)?'
        WRITE (LUOUT,*) ' (Default is Y)'
        READ (LUIN, '(A1)') ANS
        IF (ANS .EQ. 'n'  .OR.  ANS .EQ. 'N') GOTO 14
      ENDIF
      OPEN (NOUT,FILE=OFILE,FORM='FORMATTED',STATUS='UNKNOWN',
     +      ACCESS='SEQUENTIAL',ERR=9910,IOSTAT=IOS)

**********************************
* GET THE NAME FOR THE INPUT FILE
**********************************

      WRITE (LUOUT,*) ' '
   13 WRITE (LUOUT,*) ' Do you have an input data file (Y/N)?'
      WRITE (LUOUT,*) ' (Default is N)'

      READ (LUIN,'(A1)') ANS
      IF (ANS .EQ. 'Y'  .OR.  ANS .EQ. 'y') THEN
        WRITE (LUOUT,*) ' '
        WRITE (LUOUT,*) ' What is the name of the input data file?'
        WRITE (LUOUT,*) ' The default name is ''BBOOK''.'
        READ (LUIN,'(A80)') INFILE
        IF (INFILE .EQ. B80) INFILE = 'BBOOK'
        OPEN (NIN,FILE=INFILE,FORM='FORMATTED',STATUS='OLD',
     +        ACCESS='SEQUENTIAL',ERR=9920,IOSTAT=IOS)

*****************************
* CHOSE THE INPUT FILE FORMAT
*****************************

 9001   CONTINUE
        CALL NBLANK (INFILE, IFLAG2, N2)
        WRITE (LUOUT,*) ' '
        WRITE (LUOUT,*) ' What is your file format?'
        WRITE (LUOUT,*) '  0) Help - File format information'
        WRITE (LUOUT,*) '  1) Free Format Type 1'
        WRITE (LUOUT,*) '  2) Free Format Type 2'
        WRITE (LUOUT,*) '  3) NGS Blue Book (old) format (default)'
        WRITE (LUOUT,*) '  4) NGS Blue Book (new) format (make *84*)'
        WRITE (LUOUT,*) '  5) NGS Blue Book (new) format (make *86*)'
        WRITE (LUOUT,*) '  6) NGS Blue Book (extended) format'

        READ (LUIN,'(A1)') ANS
        IF (ANS .EQ. ' ') THEN

* NGS Horizontal Blue Book format

          ITYPE = 3

        ELSE
          READ (ANS,347,ERR=9940,IOSTAT=IOS) ITYPE
  347     FORMAT (I1)

          IF (ITYPE .GT. 6  .OR.  ITYPE .LT. 0) THEN

* Not a good answer - Try again.

            WRITE (LUOUT,*) ' Gotta pick one of these -',
     +                      ' sorry try again.'
            GOTO 9001

          ENDIF

        ENDIF

* Get help information

        IF (ITYPE .EQ. 0) THEN
          CALL FHELP
          GOTO 9001
        ENDIF

********************************
* CHECK FOR A SCREEN OUTPUT ALSO
********************************

        WRITE (LUOUT,*) ' '
        WRITE (LUOUT,*) ' Do you want the results written to the',
     +                  ' terminal screen as well as to'
        WRITE (LUOUT,*) ' the output file (Y/N)?'
        WRITE (LUOUT,*) ' (Default is N)'
        READ (LUIN,'(A1)') ANS
        IF (ANS .NE. 'y'  .AND.  ANS .NE. 'Y') SCREEN = .FALSE.

        GOTO 9002

* Error message

 9940   WRITE (LUOUT,*) ' Gotta pick ''1'' or ''2'' or ''3'' -',
     +                  ' sorry try again.'
        GOTO 9001

 9002 ENDIF

      RETURN

* Error message

 9910 CONTINUE
      CALL NBLANK (OFILE, IFLAG1, N1)
      CALL NBLANK (OFILE, IFLAG2, N2)
      WRITE (LUOUT,9915) IOS, OFILE(N1:N2)
 9915 FORMAT (' ERROR (', I5, ') - The operating system could not',
     +        ' open the file ', /,
     +        ' ''', A, '''', /,
     +        ' Try again.', /)
      GOTO 14

 9920 CONTINUE
      CALL NBLANK (INFILE, IFLAG1, N1)
      CALL NBLANK (INFILE, IFLAG2, N2)
      WRITE (LUOUT,9915) IOS, INFILE(N1:N2)
      GOTO 13

      END

      SUBROUTINE MLOOP (NCONV, IPAGE, ITYPE, VRSION, GHT,
     +                  SDGHT, SDGHT2, SMGHT, BGGHT,
     +                  XSMALL, XBIG, YSMALL, YBIG,
     +                  PAGE, SCREEN, NAME)

**********************************************************************
* THIS SUBROUTINE LOOPS THROUGH THE INPUT DATA (EITHER AN INPUT DATA *
* FILE OR INTERACTIVELY), CALCULATES THE GEOID HEIGHT VALUES,        *
* UPDATES THE MINIMUM, MAXIMUM, AND STATISTICAL SUMMATIONS, AND THEN *
* PRINTS THE RESULTS TO THE OUTPUT FILE AND/OR THE SCREEN.           *
**********************************************************************

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      DOUBLE PRECISION XSMALL, XBIG, YSMALL, YBIG, XPT, YPT
      DOUBLE PRECISION SDGHT, SDGHT2
      REAL VRSION, GHT, SMGHT, BGGHT
      REAL SLA, SLO, GHT2
      INTEGER NCONV, IPAGE, ITYPE
      INTEGER IDLA, IMLA, IDLO, IMLO

      CHARACTER*80 NAME
      CHARACTER*15 RESP
      LOGICAL PAGE, NOGO, SCREEN, NOPT, EOF

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

* set defaults for those variables not used by every format type

********************************************************************
* BEGIN THE COMPUTATION LOOP FOR EACH INTERPOLATION
* DO UNTIL END OF FILE OR NO MORE COMPUTATIONS REQUESTED
********************************************************************

      NCONV = 0
  160 CONTINUE

        PAGE = .FALSE.

********************************************
* GET THE NAME AND LOCATION OF ANOTHER POINT
********************************************

        CALL GETPT (NCONV, ITYPE, NAME, IDLA, IMLA, SLA, IDLO,
     +              IMLO, SLO, XPT, YPT, EOF, NOPT)
        IF (NOPT) GOTO 155
        IF (EOF) GOTO 9999

**********************
* DO THE INTERPOLATION
**********************

        NOGO = .FALSE.
        CALL INTERP (NOGO, RESP, XPT, YPT, GHT, ITYPE)

****************************************************
* CHECK TO SEE IF THIS POINT IS VALID
* IF NOGO IS TRUE THEN GET ANOTHER POINT AND DON'T
* DO THE COMPUTATION - POINT IS OUT OF BOUNDS
* IF NOGO IS NOT TRUE THEN PROCEED - ESTIMATE MADE
****************************************************

        IF (NOGO) GOTO 155
        NCONV = NCONV + 1

********************************************
* CHANGE THE LONGITUDE BACK TO POSITIVE WEST
* AND CHANGE BACK TO D.M.S FORMAT
********************************************

        XPT = -XPT

**************************
* DO THE LITTLE STATISTICS
**************************

* First, the basics
* meters....

        GHT2   = ( DBLE(GHT) )**2
        SDGHT2 = SDGHT2 + GHT2
        SDGHT  = SDGHT  + DBLE(GHT)

* then the ranges

        XSMALL = DMIN1( XSMALL, XPT)
        XBIG   = DMAX1( XBIG  , XPT)
        YSMALL = DMIN1( YSMALL, YPT)
        YBIG   = DMAX1( YBIG  , YPT)

        SMGHT = MIN( SMGHT, GHT)
        BGGHT = MAX( BGGHT, GHT)

**********************************
** WRITE TO OUTPUT FILE AND SCREEN
**********************************

        CALL WRTPT (ITYPE, NCONV, VRSION, NAME,
     +              IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +              GHT, RESP, IPAGE, PAGE, SCREEN)

**********************
* START THE LOOP AGAIN
**********************

  155 GOTO 160

 9999 RETURN
      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 IFLAG, NBLK, LENG, IBLK
      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 NGRIDS (NODATA)

* This subroutine opens the GEOID grids.
* A total of one file is necessary for each area.

* If a file named AREA.PAR exists it will be read for the names and
* locations of the gridded data.  The format of the data in
* the AREA.PAR file is given in the GRIDS subroutine.

* If the AREA.PAR file does not exist, or there is still room in the
* arrays in the GDINFO common then the default area names used.

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      CHARACTER*1 ANS
      LOGICAL NODATA

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      REAL DX, DY, XMAX, XMIN, YMAX, YMIN
      INTEGER NR, NC, NAREA
      COMMON /GDINFO/ DX(MXAREA), DY(MXAREA), XMAX(MXAREA),
     +                XMIN(MXAREA), YMAX(MXAREA), YMIN(MXAREA),
     +                NR(MXAREA), NC(MXAREA), NAREA

* Initialize

      NODATA = .FALSE.
      NAREA = 0
      WRITE (LUOUT,100)
  100 FORMAT (' GEOID is now opening the files containing the',
     +        ' gridded data.', /,
     +        ' The areas listed below may be used for GEOID',
     +        ' ESTIMATES.')

* Try to open the 'AREA.PAR' file in the subroutine GRIDS

      CALL GRIDS

* If NAREA>=MXAREA, then skip the section that opens the default files.
* If NAREA<MXAREA or no 'AREA.PAR' file exists, then open default names
* in the subroutine DGRIDS.

  130 IF (NAREA .LT. MXAREA) THEN

        CALL DGRIDS

      ENDIF

      WRITE (LUOUT,975)
  975 FORMAT (/, '  (Hit RETURN to continue.)')
      READ (LUIN,'(A1)') ANS
      WRITE (LUOUT,2)
*   2 FORMAT ('1')
    2 FORMAT ('')

      IF (NAREA .EQ. 0) THEN
        NODATA = .TRUE.
        WRITE (LUOUT, 970)
  970   FORMAT (/, ' ******* ERROR *********', /,
     +          ' No grid files were opened -- program ending!')
      ENDIF

      RETURN
      END

      SUBROUTINE OPENFL (AFILE, ITEMP, NOGO, DX1, DY1, XMAX1,
     +                   XMIN1, YMAX1, YMIN1, NR1, NC1, CARD, GFLAG)

*** Given base name of gridded data files, open the data file

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*23 B23
      PARAMETER (B23 = '                       ')
      CHARACTER*69 B69
      PARAMETER (B69 = B23//B23//B23)

      REAL DX1, DY1, XMAX1, XMIN1, YMAX1, YMIN1
      REAL X01, Y01, ANGLE1
      INTEGER LRECL, IFILE, IOS, NZ1
      INTEGER IFLAG1, IFLAG2
      INTEGER ITEMP, NR1, NC1
      INTEGER N1, N2, N3, N4
      CHARACTER*80 CARD
      CHARACTER*69 AGEO
      CHARACTER*65 AFILE
      CHARACTER*56 RIDENT
      CHARACTER*8 PGM
      LOGICAL NOGO, OFLAG, EFLAG1, GFLAG

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA OFLAG /.FALSE./, EFLAG1 /.FALSE./
      DATA IFLAG1 /1/, IFLAG2 /2/

* Initialize

      NOGO = .FALSE.

* Form complete names of grid files

      CALL NBLANK (AFILE, IFLAG2, N2)
      IF (N2 .EQ. 0) STOP 'Logical Coding Error in OPENF'
      AGEO = B69
      AGEO(1:N2) = AFILE
      AGEO(N2+1:N2+4) = '.geo'

*******************************************************
* DIRECT ACCESS GRID FILES
* Each file is opened once to get the grid variables.
* The file is then closed and reopened to ensure that
* the record length is correct
*******************************************************

      LRECL = 256
      IFILE = 10 + ITEMP
      LUAREA(ITEMP) = IFILE
      INQUIRE (FILE=AGEO, EXIST=EFLAG1, OPENED=OFLAG)
      IF (.NOT. EFLAG1) GOTO 910
      IF (OFLAG) GOTO 980
      OPEN (IFILE,FILE=AGEO,FORM='UNFORMATTED',STATUS='OLD',
     +       ACCESS='DIRECT',RECL=LRECL,ERR=910,IOSTAT=IOS)
      READ (IFILE,REC=1) RIDENT, PGM, NC1, NR1, NZ1, X01, DX1,
     +                   Y01, DY1, ANGLE1
      CLOSE (IFILE)

      LRECL = 4*(NC1+1)
      OPEN (IFILE,FILE=AGEO,FORM='UNFORMATTED',STATUS='OLD',
     +      ACCESS='DIRECT',RECL=LRECL,ERR=910,IOSTAT=IOS)

* Calculate values used in this program

      XMIN1 = X01
      YMIN1 = Y01
      XMAX1 = X01 + (NC1-1)*DX1
      YMAX1 = Y01 + (NR1-1)*DY1
      DX1 = ABS(DX1)
      DY1 = ABS(DY1)

*****************************************
* REPORT SOMETHING ABOUT THE GRIDDED DATA
*****************************************
*     WRITE (LUOUT,4050) RIDENT, PGM, NC1, NR
*4050 FORMAT (1X, A56, /, 1X, A8, /, I5, I5)
*     WRITE (LUOUT,*) 'DX,DY,NR,NC', DX1, DY1, NR1, NC1
*     WRITE (LUOUT,4055) -XMAX1, -XMIN1, YMIN1, YMAX1
*4055 FORMAT (' MIN Longitude = ', F10.4, ' MAX Longitude = ', F10.4, /,
*    +        ' MIN Latitude  = ', F10.4, ' MAX Latitude  = ', F10.4, /)
*****************************************

 9999 RETURN

****************************
* WARNING AND ERROR MESSAGES
****************************

* Grid file does not exist

  910 CONTINUE
      NOGO = .TRUE.

      IF (GFLAG) THEN
        CALL NBLANK (AGEO, IFLAG1, N1)
        CALL NBLANK (CARD, IFLAG1, N3)
        CALL NBLANK (AGEO, IFLAG2, N2)
        CALL NBLANK (CARD, IFLAG2, N4)

* .geo does not exist

        WRITE (LUOUT, 915) AGEO(N1:N2), CARD(N3:N4)
  915   FORMAT (/, 5X, ' *********** WARNING ***********', /,
     +          5X, ' The grid file', /,
     +          6X, '''', A, '''', /,
     +          5X, ' from record:', /,
     +          6X, '''', A, '''', /,
     +          5X, ' does not exist!', /,
     +          5X, ' *******************************', /)
      ENDIF

      CLOSE (IFILE)
      GOTO 9999

* Grid file already open

  980 CONTINUE
      NOGO = .TRUE.

      CALL NBLANK (AGEO, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG1, N3)
      CALL NBLANK (AGEO, IFLAG2, N2)
      CALL NBLANK (CARD, IFLAG2, N4)

      WRITE (LUOUT, 985) AGEO(N1:N2)
  985 FORMAT (/, 5X, ' *********** ERROR ***********', /,
     +        5X, ' The grid file', /,
     +        6X, '''', A, '''', /,
     +        5X, ' has already been opened!  This file',
     +                                      ' will not be reopened.', /,
     +        5X, ' *****************************', /)

      GOTO 9999
      END

      SUBROUTINE PRINT (LU, NCONV, NAME, VRSION, IDLA, IMLA, SLA,
     +                  IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE)

* This subroutine prints out the actual results using
* a pretty format - not the same as the input file format (if there
* is one).  This subroutine is used by type-1 format input and
* interactive input.

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

      REAL VRSION
      REAL SLA, SLO, GHT
      INTEGER LU, NCONV, IPAGE
      INTEGER IDLA, IMLA, IDLO, IMLO
      CHARACTER*80 NAME
      CHARACTER*15 RESP
      LOGICAL PAGE

      IF (NCONV .EQ. 1) THEN

********************
* FIRST PAGE HEADING
********************

        WRITE (LU,10) IPAGE
        WRITE (LU,5)
        WRITE (LU,26)
        WRITE (LU,7) VRSION
        WRITE (LU,8)
      ENDIF

      IF (PAGE) THEN
        IF (IPAGE .GT. 1) THEN
          WRITE (LU,2)
          WRITE (LU,10) IPAGE
        ENDIF
      ENDIF

   10 FORMAT (70(' '), 'Page ', I4, /)
    5 FORMAT (20X, '      Geoid Height Computation')
   26 FORMAT (20X, '        GEOID93 Interpolation'   )
    7 FORMAT (20X, '         GEOID Program ', /, 20X,
     +             '          Version ', F4.2 )
    8 FORMAT (20X, '                          ', /, 1X, 79('=') )

*   2 FORMAT ('1')
    2 FORMAT ('')

      WRITE (LU,921) NCONV, RESP
      IF (NAME .NE. '    ') WRITE (LU,922) NAME
      WRITE (LU,900)
      WRITE (LU,927)
      WRITE (LU,923) IDLA, IMLA, SLA, IDLO, IMLO, SLO, GHT

  921 FORMAT (//, 27X, 'Computation #: ', I4, 8X, 'Region: ', A15, /)
  922 FORMAT (2X, 'Station Name:  ', A80, /)
  900 FORMAT (22X, 'Latitude', 17X, 'Longitude', 7X, 'Geoid Height')
  927 FORMAT (65X, '(meters)', /)
  923 FORMAT (19X, I2, 1X, I2.2, F9.5, 10X, I3, 1X, I2.2, F9.5,
     +        5X, F9.3, /)

      RETURN
      END

      SUBROUTINE PRINT2 (LU, NCONV, VRSION, GHT)

* This subroutine prints out the actual results using
* a free format - the same as the input file format.  This is used
* for type 2 format.

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

      REAL VRSION, GHT
      INTEGER LU, NCONV

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

* Write header record to identify source of geoid height and value

      IF (NCONV .EQ. 1) THEN
        WRITE (LU, 10) VRSION
   10   FORMAT ('GEOID      Version', F5.2,
     +          ' - NGS GEOID93 Model')
      ENDIF

* In this format, the variable CARD contains the image of the input
* card.  This is written to the output file instead of using the
* latitude, longitude, and name variables.  The geoid height
* overwrites whatever is in columns 33-40

      WRITE (LU,100) CARD(1:32), GHT, CARD(41:80)
  100 FORMAT (A32, F8.3, A40)

      RETURN
      END

      SUBROUTINE PRINT3 (ITYPE, LU, GHT, RESP)

* This subroutine prints out the actual results into
* a *84* or *86* record.  This is used for type 3 format.

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

      CHARACTER*15 B15
      PARAMETER (B15 = '               ')

      REAL GHT
      INTEGER LU, IGHT, ITYPE
      CHARACTER*125 CRD125
      CHARACTER*80 CRD80,card2
      CHARACTER*15 RESP
      CHARACTER*10 ASTD2
      CHARACTER*4 ASTD

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

*** Standard deviation is assumed to be 0.5 meters

      DATA ASTD /'   5'/, ASTD2 /'    0.5000'/

* In this format, the variable CARD contains the image of the input
* card.  This is written to the output file instead of using the
* latitude, longitude, and name variables.  The geoid height
* overwrites whatever is in columns 33-40

      IF (ITYPE .EQ. 3) THEN

* Old Blue Book

        CRD80(1:7) = CARD(1:7)
        CRD80(8:9) = '84'
        CRD80(10:14) = CARD(10:14)
        CRD80(15:20) = 'NGS  '
        CRD80(21:71) = 'PROGRAM GEOID      '//RESP//B15//'  '
        IGHT = NINT(GHT*10.)
        WRITE (CRD80(72:76),'(I5)') IGHT
        CRD80(77:80) = ASTD
        WRITE (LU,'(A80)') CRD80

      ELSEIF (ITYPE .EQ. 4) THEN

* New Blue Book  *84*

        CRD80(1:7) = CARD(1:7)
        CRD80(8:9) = '84'
        CRD80(10:14) = CARD(10:14)
        CRD80(15:20) = 'NGS  '
        CRD80(21:69) = 'PROGRAM GEOID      '//RESP//B15
        IGHT = NINT(GHT*100.)
        WRITE (CRD80(70:75),'(I6)') IGHT
        CRD80(76:80) = ASTD//'0'
        WRITE (LU,'(A80)') CRD80

*******************************************************************

      ELSEIF (ITYPE .EQ. 5) THEN

* New Blue Book  *86*

*** look ahead for an exisiting 86 record
*** kludgy because no clean way in standard fortran to check if at eof

        read(nin,'(a80)',end=66565) card2
        if(card2(8:9).eq.'86'.and.card2(11:14).eq.card(11:14)) then
          crd80=card2
          go to 66566
        else
          crd80='                                                      '
          backspace(nin)
          go to 66566
        endif
        stop 34287
*************** eof error trap ***************************
66565   crd80='                                                        '

66566   CRD80(1:7) = CARD(1:7)
        CRD80(8:9) = '86'
        CRD80(10:14) = CARD(10:14)
        IGHT = NINT(GHT*1000.0)
        WRITE (CRD80(36:42),'(I7)') IGHT
        CRD80(43:43) = 'D'
        WRITE (LU,'(A80)') CRD80

*******************************************************************

      ELSEIF (ITYPE .EQ. 6) THEN

* Extended Blue Book

        CRD125(1:7) = CARD(1:7)
        CRD125(8:9) = '84'
        CRD125(10:14) = CARD(10:14)
        CRD125(15:20) = 'NGS  '
        CRD125(21:71) = 'PROGRAM GEOID      '//RESP//B15//'  '
        IGHT = NINT(GHT*10.)
        WRITE (CRD125(72:76),'(I5)') IGHT
        CRD125(77:80) = ASTD
        CRD125(81:101) = B15//'     '
        WRITE (CRD125(102:115),'(F14.4)') GHT
        CRD125(116:125) = ASTD2
        WRITE (LU,'(A125)') CRD125
      ENDIF

      RETURN
      END

      SUBROUTINE REPORT (LU, SMGHT, BGGHT, AGHT, VGHT, SDGHT,
     +                   IPAGE, PAGE)

* This subroutine prints out the statistics

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

      DOUBLE PRECISION AGHT, VGHT, SDGHT
      REAL SMGHT, BGGHT
      INTEGER LU, IPAGE
      LOGICAL PAGE

************
* NEW PAGE
************

      WRITE (LU,2)
*   2 FORMAT ('1')
    2 FORMAT ('')

      IF (PAGE) THEN

******************
* NUMBER THIS PAGE
******************

        WRITE (LU,15)  IPAGE
   15   FORMAT (70(' '), 'Page ', I4, /)
      ENDIF

        WRITE (LU,90)
   90   FORMAT (28X, 'Geoid Height Computation')
        WRITE (LU,100)
  100   FORMAT (//, 30X, 'Statistics for Region', /, 1X, 79('='), /)
        WRITE (LU,910)
  910   FORMAT (/,'  Geoid Height           ',
     +          8X, 'MIN', 6X, 'MAX' )
        WRITE (LU,110) SMGHT, BGGHT
  110   FORMAT ('  Range of Height (meters)  ', 2F9.3)

        WRITE (LU,130) AGHT
  130   FORMAT (/, 10X, 'Mean Geoid Height        ', F9.3)
        WRITE (LU,131) VGHT
  131   FORMAT (   10X, 'Variance of Geoid Height ', F9.3)
        WRITE (LU,132) SDGHT
  132   FORMAT (   10X, 'Std. Dev. of Geoid Height', F9.3)

      RETURN
      END

      SUBROUTINE TYPE1 (NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +                  XPT, YPT, EOF, NOPT)

* Read a record from a file of type 1. In this type there is a station
* name (or blanks) in columns 1-40, and free-format latitude and
* longitude values in columns 41-80.  By free format we mean that the
* numbers making up the degrees, minutes and seconds of latitude,
* degrees, minutes, seconds of longitude must appear in that order in
* columns 41 through 80 but are not restricted to any specific columns.
* The latitude and longitude may be either (1) integer degrees, integer
* minutes, decimal seconds, or (2) integer degrees, decimal minutes, or
* (3) decimal degrees.

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*1 DOT, BLK
      PARAMETER ( DOT = '.', BLK = ' ' )
      CHARACTER*80  B80
      CHARACTER*20 B20
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD
      REAL SLA, SLO, RCARD, RMLA, RMLO
      INTEGER IFLAG1, IFLAG2, N1, N2
      INTEGER IDOT, IBLK, LENG, IERR
      INTEGER IDLA, IMLA, IDLO, IMLO
      CHARACTER*80 NAME
      CHARACTER*40 DUMLA, DUMLO
      LOGICAL EOF, NOPT

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA IFLAG1 /1/, IFLAG2 /2/

***********************************
* FOR INPUT FILE OF ITYPE = 1
***********************************

    1 READ (NIN,'(A80)',END=9999) CARD
      READ (CARD(1:40), '(A40)') NAME

* Check for blank line

      IF (CARD .EQ. B80) GOTO 1

* Find position of the first decimal point (to indicate the last
* number in the latitude)

      IDOT = INDEX(CARD(41:80), DOT)

* Error - no decimal point

      IF (IDOT .EQ. 0) GOTO 9980

* find position of the first blank after the first decimal point (to
* indicate the blank after the last number in the latitude)

      IDOT = IDOT + 40
      IBLK = INDEX(CARD(IDOT+1:80), BLK)
      IBLK = IBLK + IDOT

      DUMLA = CARD(41:IBLK)
      LENG = IBLK - 41
      RDLA = DCARD( DUMLA, LENG, IERR )
      IF (IERR .NE. 0) GOTO 9950
      IF (LENG .GT. 0) THEN

        RMLA = RCARD( DUMLA, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9950

        IF (LENG .GT. 0) THEN
          SLA  = RCARD( DUMLA, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9950
        ELSE
          SLA = 0.E0
        ENDIF

      ELSE
        RMLA = 0.E0
        SLA = 0.E0
      ENDIF

* Check for illogical values

      IF (RDLA .LT.   0.D0) GOTO 9940
      IF (RDLA .GT.  90.D0) GOTO 9950
      IF (RMLA .LT. 0.E0  .OR.  RMLA .GT. 60.E0) GOTO 9950
      IF ( SLA .LT. 0.E0  .OR.   SLA .GT. 60.E0) GOTO 9950

***********
* LONGITUDE
***********

      DUMLO = CARD(IBLK+1:80)
      CALL NBLANK (DUMLO, IFLAG2, N2)
      LENG = N2
      RDLO = DCARD( DUMLO, LENG, IERR )
      IF (IERR .NE. 0) GOTO 9960
      IF (LENG .GT. 0) THEN

        RMLO = RCARD( DUMLO, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9960

        IF (LENG .GT. 0) THEN
          SLO  = RCARD( DUMLO, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9960
        ELSE
          SLO = 0.E0
        ENDIF

      ELSE
        RMLO = 0.E0
        SLO = 0.E0
      ENDIF

* Check for illogical values

      IF (RDLO .LT.   0.D0) GOTO 9940
      IF (RDLO .GT. 360.D0) GOTO 9960
      IF (RMLO .LT. 0.E0  .OR.  RMLO .GT. 60.E0) GOTO 9960
      IF ( SLO .LT. 0.E0  .OR.   SLO .GT. 60.E0) GOTO 9960

* Calculate decimal degrees

      YPT = RDLA + DBLE(RMLA)/60.D0 + DBLE(SLA)/3600.D0
      XPT = RDLO + DBLE(RMLO)/60.D0 + DBLE(SLO)/3600.D0

* Get degrees, minutes, seconds

      CALL HMS (YPT, IDLA, IMLA, SLA)
      CALL HMS (XPT, IDLO, IMLO, SLO)

 9000 RETURN

* Error messages

 9940 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9945) CARD(N1:N2)
 9945 FORMAT (' ERROR - in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude and Longitudes must be positive!', /,
     +        '         Longitude is positive west.', /)
      NOPT = .TRUE.
      GOTO 9000

 9950 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9955) CARD(N1:N2)
 9955 FORMAT (' ERROR - Illogical values for latitude',
     +        ' in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude must be between 0 and 90 degrees.', /,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9960 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9965) CARD(N1:N2)
 9965 FORMAT (' ERROR - Illogical values for longitude',
     +        ' in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Longitude must be between 0 and 360 degrees.',/,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9980 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9985) CARD(N1:N2)
 9985 FORMAT (' ERROR - The following record does not have a decimal',
     +        ' point in the latitude.', /,
     +        9X, '''', A, '''', /,
     +        '         In the free format a decimal point is used',
     +        ' to determine what is', /,
     +        '         the last number in the latitude.  Please',
     +        ' correct this record', /,
     +        '         and check all of the data in this file to',
     +        ' ensure that it follows', /,
     +        '         the correct format.', /)
      NOPT = .TRUE.
      GOTO 9000

 9999 CONTINUE
      EOF = .TRUE.
      GOTO 9000
      END

      SUBROUTINE TYPE2 (IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +                  XPT, YPT, EOF, NOPT)

* Read a record from a file of type 2. In this type there is free-format
* latitude and longitude values in columns 1-32, and a station name
* (or blanks) in columns 41-80.  Columns 33-40 will contain the geoid
* height in the output record and are ignored in the input record.
* By free format we mean that the numbers making up the degrees,
* minutes and seconds of latitude, degrees, minutes, seconds of
* longitude must appear in that order in columns 1 through 32 but are
* not restricted to any specific columns.  The latitude and longitude
* may be either (1) integer degrees, integer minutes, decimal seconds,
* or (2) integer degrees, decimal minutes, or (3) decimal degrees.

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)
      CHARACTER*1 DOT, BLK
      PARAMETER ( DOT = '.', BLK = ' ' )
      CHARACTER*20 B20
      CHARACTER*80 B80
      PARAMETER (B20 = '                   ', B80 = B20//B20//B20//B20)

      DOUBLE PRECISION XPT, YPT, RDLA, RDLO, DCARD
      REAL SLA, SLO, RCARD, RMLA, RMLO
      INTEGER IFLAG1, IFLAG2, N1, N2
      INTEGER IDOT, IBLK, LENG, IERR
      INTEGER IDLA, IMLA, IDLO, IMLO
      CHARACTER*32 DUMLA, DUMLO
      LOGICAL EOF, NOPT

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA IFLAG1 /1/, IFLAG2 /2/

***********************************
* FOR INPUT FILE OF ITYPE = 2
***********************************

    1 READ (NIN,'(A80)',END=9999) CARD

* Check for blank line

      IF (CARD .EQ. B80) GOTO 1

* Find position of the first decimal point (to indicate the last
* number in the latitude)

      IDOT = INDEX(CARD(1:32), DOT)

* Error - no decimal point

      IF (IDOT .EQ. 0) GOTO 9980

* find position of the first blank after the first decimal point (to
* indicate the blank after the last number in the latitude)

      IBLK = INDEX(CARD(IDOT+1:32), BLK)
      IBLK = IBLK + IDOT

      DUMLA = CARD(1:IBLK)
      LENG = IBLK - 1
      RDLA = DCARD( DUMLA, LENG, IERR )
      IF (IERR .NE. 0) GOTO 9950
      IF (LENG .GT. 0) THEN

        RMLA = RCARD( DUMLA, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9950

        IF (LENG .GT. 0) THEN
          SLA  = RCARD( DUMLA, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9950
        ELSE
          SLA = 0.E0
        ENDIF

      ELSE
        RMLA = 0.E0
        SLA = 0.E0
      ENDIF

* Check for illogical values

      IF (RDLA .LT.   0.D0) GOTO 9940
      IF (RDLA .GT.  90.D0) GOTO 9950
      IF (RMLA .LT. 0.E0  .OR.  RMLA .GT. 60.E0) GOTO 9950
      IF ( SLA .LT. 0.E0  .OR.   SLA .GT. 60.E0) GOTO 9950

***********
* LONGITUDE
***********

      DUMLO = CARD(IBLK+1:32)
      CALL NBLANK (DUMLO, IFLAG2, N2)
      LENG = N2
      RDLO = DCARD( DUMLO, LENG, IERR )
      IF (IERR .NE. 0) GOTO 9960
      IF (LENG .GT. 0) THEN

        RMLO = RCARD( DUMLO, LENG, IERR )
        IF (IERR .NE. 0) GOTO 9960

        IF (LENG .GT. 0) THEN
          SLO  = RCARD( DUMLO, LENG, IERR )
          IF (IERR .NE. 0) GOTO 9960
        ELSE
          SLO = 0.E0
        ENDIF

      ELSE
        RMLO = 0.E0
        SLO = 0.E0
      ENDIF

* Check for illogical values

      IF (RDLO .LT.   0.D0) GOTO 9940
      IF (RDLO .GT. 360.D0) GOTO 9960
      IF (RMLO .LT. 0.E0  .OR.  RMLO .GT. 60.E0) GOTO 9960
      IF ( SLO .LT. 0.E0  .OR.   SLO .GT. 60.E0) GOTO 9960

* Calculate decimal degrees

      YPT = RDLA + DBLE(RMLA)/60.D0 + DBLE(SLA)/3600.D0
      XPT = RDLO + DBLE(RMLO)/60.D0 + DBLE(SLO)/3600.D0

* Get degrees, minutes, seconds

      CALL HMS (YPT, IDLA, IMLA, SLA)
      CALL HMS (XPT, IDLO, IMLO, SLO)

 9000 RETURN

* Error messages

 9940 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9945) CARD(N1:N2)
 9945 FORMAT (' ERROR - in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude and Longitudes must be positive!', /,
     +        '         Longitude is positive west.', /)
      NOPT = .TRUE.
      GOTO 9000

 9950 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9955) CARD(N1:N2)
 9955 FORMAT (' ERROR - Illogical values for latitude',
     +        ' in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude must be between 0 and 90 degrees.', /,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9960 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9965) CARD(N1:N2)
 9965 FORMAT (' ERROR - Illogical values for longitude',
     +        ' in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Longitude must be between 0 and 360 degrees.',/,
     +        '         Minutes and seconds must be between 0',
     +                                                    ' and 60.', /)
      NOPT = .TRUE.
      GOTO 9000

 9980 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9985) CARD(N1:N2)
 9985 FORMAT (' ERROR - The following record does not have a decimal',
     +        ' point in the latitude.', /,
     +        9X, '''', A, '''', /,
     +        '         In the free format a decimal point is used',
     +        ' to determine what is', /,
     +        '         the last number in the latitude.  Please',
     +        ' correct this record', /,
     +        '         and check all of the data in this file to',
     +        ' ensure that it follows', /,
     +        '         the correct format.', /)
      NOPT = .TRUE.
      GOTO 9000

 9999 CONTINUE
      EOF = .TRUE.
      GOTO 9000
      END

      SUBROUTINE TYPE3 (ITYPE, NAME, IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +                  XPT, YPT, EOF, NOPT)

* Read a record from a file of type 3 ('old' Blue Book) or
* of type 4 or 5 ('new' Blue Book) or type 6 (extended Blue Book).
* These formats are defined in 'Input Formats and Specifications of
* the' National Geodetic Survey Data Base', Volume 1. Horizontal Control
* Data.  For type 3 the publishing date was 1980.  For type 4 the
* publishing date was January, 1989.
*** Type 5 uses the new *86* record
*  Type 6 is for internal NGS use.
* This publication is available from the National Geodetic Survey for a
* fee by calling (301) 443-8631

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      DOUBLE PRECISION XPT, YPT
      REAL SLA, SLO
      INTEGER ITYPE
      INTEGER IDLA, IMLA, IDLO, IMLO
      INTEGER IFLAG1, IFLAG2, N1, N2
      CHARACTER*80 NAME
      CHARACTER*1 DIRLAT, DIRLON
      LOGICAL EOF, NOPT

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

      DATA IFLAG1 /1/, IFLAG2 /2/

***********************************
* FOR INPUT FILE OF ITYPE = 3
***********************************

      READ (NIN,6000,END=9999) CARD
 6000 FORMAT (A80)

*** strip any preexisting 84 and 86  (retained in extended format)

      IF((CARD(8:9).NE.'84'.AND.CARD(8:9).NE.'86').AND.ITYPE.NE.6)THEN

* For ITYPE = 3 and ITYPE 4 copy input non-*80*, non-*84* record
* to output

        WRITE (NOUT,6000) CARD
      ENDIF

      IF (CARD(8:9) .EQ. '80') THEN

* The station names and locations are in the *80* records

        IF (ITYPE .EQ. 3) THEN
          READ (CARD,6003) NAME, IDLA, IMLA, SLA, DIRLAT,
     +                     IDLO, IMLO, SLO, DIRLON
 6003     FORMAT (BZ, 10X, 3X, 1X, A30, I2, I2, F7.5, A1,
     +            I3, I2, F7.5, A1)
        ELSEIF (ITYPE .EQ. 4 ) THEN
          READ (CARD,6004) NAME, IDLA, IMLA, SLA, DIRLAT,
     +                     IDLO, IMLO, SLO, DIRLON
 6004     FORMAT (BZ, 10X, 4X, A30, I2, I2, F7.5, A1,
     +            I3, I2, F7.5, A1)

        ELSEIF (itype.eq.5) THEN
          READ (CARD,6004) NAME, IDLA, IMLA, SLA, DIRLAT,
     +                     IDLO, IMLO, SLO, DIRLON

        ELSEIF (ITYPE .EQ. 6) THEN
          READ (CARD,6005) NAME, IDLA, IMLA, SLA, DIRLAT,
     +                     IDLO, IMLO, SLO, DIRLON
 6005     FORMAT (BZ, 9X, 5X, A30, I2, I2, F7.5, A1,
     +            I3, I2, F7.5, A1)
        ENDIF

* Check for illogical values

        IF (DIRLAT .NE. 'N'  .AND.  DIRLAT .NE. 'n') GOTO 9940
        IF (IDLA .LT.   0) GOTO 9940
        IF (IDLA .GT.  90) GOTO 9950
        IF (IMLA .LT. 0     .OR.  IMLA .GT. 60   ) GOTO 9950
        IF ( SLA .LT. 0.E0  .OR.   SLA .GT. 60.E0) GOTO 9950
        IF (DIRLON .NE. 'W'  .AND.  DIRLON .NE. 'w') GOTO 9940
        IF (IDLO .LT.   0) GOTO 9940
        IF (IDLO .GT. 360) GOTO 9950
        IF (IMLO .LT. 0     .OR.  IMLO .GT. 60   ) GOTO 9950
        IF ( SLO .LT. 0.E0  .OR.   SLO .GT. 60.E0) GOTO 9950

      ELSE
        NOPT = .TRUE.
        GOTO 9000
      ENDIF

      YPT = DBLE(IDLA) + DBLE(IMLA)/60.D0 + DBLE(SLA)/3600.D0
      XPT = DBLE(IDLO) + DBLE(IMLO)/60.D0 + DBLE(SLO)/3600.D0

 9000 RETURN

* Error messages

 9940 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9945) CARD(N1:N2)
 9945 FORMAT (' ERROR - in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         Latitude and Longitudes must be positive!', /,
     +        '         Longitude is positive west.', /)
      NOPT = .TRUE.
      GOTO 9000

 9950 CONTINUE
      CALL NBLANK (CARD, IFLAG1, N1)
      CALL NBLANK (CARD, IFLAG2, N2)
      WRITE (LUOUT,9955) CARD(N1:N2)
 9955 FORMAT (' ERROR - Illogical values for latitude or longitude',
     +        ' in the following record:', /,
     +        9X, '''', A, '''', /,
     +        '         This record will be skipped.', /)
      NOPT = .TRUE.
      GOTO 9000

 9999 CONTINUE
      EOF = .TRUE.
      GOTO 9000
      END

      SUBROUTINE WRTPT (ITYPE, NCONV, VRSION, NAME,
     +                  IDLA, IMLA, SLA, IDLO, IMLO, SLO,
     +                  GHT, RESP, IPAGE, PAGE, SCREEN)

*** Write the computations to output file
***(and screen).

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

      INTEGER MXAREA
      PARAMETER (MXAREA =10)

      REAL VRSION, GHT
      REAL SLA, SLO
      INTEGER ITYPE, NCONV, IPAGE
      INTEGER IDLA, IMLA, IDLO, IMLO
      CHARACTER*80 NAME
      CHARACTER*15 RESP
      LOGICAL PAGE, SCREEN

      CHARACTER*80 CARD
      COMMON /CURNT/ CARD

      INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
      COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(MXAREA)

*********************
* PAGE NUMBER COUNTER
*********************

* this is where you change how many on a page

        IF ( MOD(NCONV,5) .EQ. 0  .OR.  NCONV .EQ. 1) THEN
          PAGE = .TRUE.
          IPAGE = IPAGE + 1
        ENDIF

***********************
** WRITE TO OUTPUT FILE
***********************

        IF (ITYPE .EQ. 0) THEN

**************************************
* ONLY INTERACTIVE USE - NO INPUT FILE
**************************************

          CALL PRINT (NOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA,
     +                IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE)

        ELSEIF (ITYPE .EQ. 1) THEN

**************************
* FOR FREE FORMAT TYPE 1
**************************

          CALL PRINT (NOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA,
     +                IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE)

        ELSEIF (ITYPE .EQ. 2) THEN

**************************
* FOR FREE FORMAT TYPE 2
**************************

          CALL PRINT2 (NOUT, NCONV, VRSION, GHT)

        ELSEIF (ITYPE .EQ. 3  .OR.  ITYPE .EQ. 4  .OR.
     +          ITYPE .EQ. 5  .or.   ITYPE .EQ. 6) THEN

**********************************************************
* FOR ITYPE = 3 THE HORIZONTAL BLUE BOOK OLD SPECIFICATION
* FOR ITYPE = 4 THE HORIZONTAL BLUE BOOK NEW SPECIFICATION *84*
* FOR ITYPE = 5 THE HORIZONTAL BLUE BOOK NEW SPECIFICATION *86*
* FOR ITYPE = 6 - THE EXTENDED BLUE BOOK SPECIFICATION
**********************************************************

          CALL PRINT3 (ITYPE, NOUT, GHT, RESP)

        ENDIF

*******************
* FOR SCREEN OUTPUT
*******************

        IF (SCREEN) THEN
          CALL PRINT (LUOUT, NCONV, NAME, VRSION, IDLA, IMLA, SLA,
     +                IDLO, IMLO, SLO, GHT, RESP, IPAGE, PAGE)
        ENDIF

      RETURN
      END

      CHARACTER*(*) FUNCTION CCARD (CHLINE, LENG, IERR)

*** Read a character variable from a line of card image.
*** LENG is the length of the card
*** blanks are the delimiters of the character variable

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

      INTEGER LENG, IERR, I, J, ILENG
      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 (assume length of 2)

        IF ( I .GE. LENG) THEN

         CCARD = '  '
         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 character string, it can be any length
* up to the length of the line

      ILENG = J - I + 1

      IF (ILENG .GT. LENG) THEN
        STOP 'CCARD'
      ENDIF

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

      READ (CHLINE(I:J), 55, ERR=9999) CCARD
   55 FORMAT (A)

* 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

      DOUBLE PRECISION FUNCTION DCARD (CHLINE, LENG, IERR)

*** Read a double precision number from a line of card image.
*** LENG is the length of the card
*** blanks are the delimiters of the REAL*8 variable

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

      DOUBLE PRECISION VAR
      INTEGER LENG, IERR, I, J, ILENG
      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
          DCARD = 0.0D0
          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. 20) THEN
        STOP 'DCARD'
      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 (F20.0)
      DCARD = 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

      INTEGER FUNCTION ICARD (CHLINE, LENG, IERR)

*** Read an integer from a line of card image.
*** LENG is the length of the card
*** blanks are the delimiters of the integer

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

      INTEGER LENG, IERR, I, J, ILENG, IVAR
      CHARACTER*80 CHLINE

      IERR = 0

* Find first non-blank character

* DO WHILE line character is bland, 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
          ICARD = 0
          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 integer string, it cannot be greater
* than 13 characters

      ILENG = J - I + 1

      IF (ILENG .GT. 13) THEN
        STOP 'ICARD'
      ENDIF

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

      READ (CHLINE(I:J), 55, ERR=9999) IVAR
   55 FORMAT (I13)
      ICARD = IVAR

* Now reset the values for 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

      REAL FUNCTION QTERP2 (X, F0, F1, F2)

*** Linear Quadratic interpolation from equally spaced (real) values
*** Uses Newton-Gregory forward polynomial
*** X ranged from 0 through 2.  (thus S = X)

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

      REAL X, F0, F1, F2
      REAL DF0, DF1, D2F0

*** forward differences

      DF0  = F1  - F0
      DF1  = F2  - F1
      D2F0 = DF1 - DF0

      QTERP2 = F0 + X*DF0 + 0.5*X*(X-1.)*D2F0

      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
