program deflec

***********************************************************************
*                                                                     *
* program :   deflec                                                  *
*                                                                     *
* purpose:    computation program to interpolate deflection           *
*             values from a grid of regularly-spaced estimates.       *
*             the input requires latitude and longitude values        *
*             referrenced to the north american datum of 1983 (nad83),*
*             for the conterminous U.S.                               *
*             input data referrenced to the north american datum of   *
*             1927 (nad27) should be converted to nad 83 via nadcon   *
*             prior to using this program, for the USA.  NAD 27 is    *
*             acceptable for MEXICO and the CARIBBEAN.                *
*                                                                     *
*             the actual computation is performed as an interpolation *
*             from a regularly-spaced grid of points obtained from    *
*             the deflec96 model computed by dr. dennis g. milbert    *
*             and dr. dru a. smith.                                   *
*                                                                     *
*         VERSION 3.10:                                               *
*             MODELS DMEX97 (computed by Dr. Dru A. Smith) and        *    
*             DCAR97 (computed by Dru A. Smith and                    *
*             Howard J. Small (NIMA)) have been added                 *
*             NOTE THAT THE MODELS DMEX97 AND DCAR97 ARE              *
*             ALL RELATIVE TO THE GRS-80 ELLIPSOID, CENTERED          *
*             AT THE ORIGIN OF THE ITRF94(1996.0) REFERENCE FRAME.    *
*             YOU CAN *STILL* INPUT LAT/LON IN NAD83 AND EXPECT TO    *
*             GET THE SAME VALUE AS IF YOU HAD USED GRS-80/ITRF94     *
*             LAT/LON VALUES, HOWEVER IT *MUST* BE POINTED OUT THAT   *
*             THE GEOID HEIGHTS OF THOSE MODELS ARE *NOT* GEOID       *
*             HEIGHTS ABOVE NAD83!  (THE REASON THIS WORKS IS BECAUSE *
*             THE HORIZONTAL SHIFT FROM NAD83 TO ITRF94 IS ABOUT 1 OR *
*             2 METERS, WHICH WILL HAVE *NO* IMPACT ON THE GEOID      *
*             HEIGHT VALUE; HOWEVER THE VERTICAL SHIFT FROM ITRF94    *
*             TO NAD83 IS ALSO 1-2 METERS, SO IT IS IMPORTANT TO      *
*             REMEMBER WHICH ELLIPSOID THE GEOIDS REFER TO.           *
*                                                                     *
*                                                                     *
*             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:  3.00   (20-oct-96)                                   *
* version code:  3.10   ( 8-apr-97)                                   *
*                                                                     *
* version date:  may 7, 1997                                          *
*                july 16, 1997 -- retained version # as 3.10          *
*                                                                     *
*   added support for inexact header in GRIDS()                       *
*   increaded MXAREA to 18                                            *
*                                                                     *
*       authors: dru a. smith, ph.d.                                  *
*                dennis g. milbert, ph.d.                             *
*                systems development division                         *
*                national geodetic survey, nos, noaa                  *
*                silver spring, md   20910                            *
*                                                                     *
* original authors: (software)                                        *
*       author:  warren t. dewhurst, ph.d.    (software)              *
*                  lieutenant commander, noaa                         *
*                alice r. drew                (software)              *
*                  senior geodesist, horizontal network branch        *
***********************************************************************

***********************************************************************
*                                                                     *
*                  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 = 3.10e0, mxarea =18)

      double precision xsmall, xbig, ysmall, ybig
      double precision axii, vxii
      double precision sdxii, sdxii2
      real xii, smxii, bgxii
      double precision aeta, veta
      double precision sdeta, sdeta2
      real eta, smeta, bgeta
      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(2*mxarea)
      
* v 3.10: differentiate between deflec96, dmex97 and dcar97, das 
      logical ldmex,ldcar
      common/scienc/ldmex,ldcar     

**********************
* initialize variables
**********************

      call initl (screen, page, name, ipage, itype,
     +            smxii, bgxii, axii, vxii, sdxii, sdxii2,
     +            smeta, bgeta, aeta, veta, sdeta, sdeta2,
     +            xsmall, xbig, ysmall, ybig)

**************************
* print header information
**************************

      call headr (vrsion)

***********************
* open deflec data files
***********************
      
      ldmex=.false.
      ldcar=.false.
      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,
     +            xii, sdxii, sdxii2, smxii, bgxii,
     +            eta, sdeta, sdeta2, smeta, bgeta,
     +            xsmall, xbig, ysmall, ybig,
     +            page, screen, name)

***********************************************
* finished with all computations - write report
***********************************************

      call endrep (ipage, nconv, itype,
     +             xii, axii, vxii, sdxii, sdxii2, smxii, bgxii,
     +             eta, aeta, veta, sdeta, sdeta2, smeta, bgeta,
     +             xsmall, xbig, ysmall, ybig)

*****************
* close all files
*****************

      do 1010 ifile = 1, 2*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 =18)
      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(2*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.', /)
        WRITE (LUOUT,1111)
 1111   FORMAT (' *************** Which latitude/longitude to use? ',
     +          ' ***************',/,
     +    ' You may enter latitudes and longitudes in either the',
     +    ' NAD 83 or the',/,
     +    ' ITRF94/GRS-80 system, with no change in the interpolated',
     +    ' deflection value.',/,
     +    ' Do not use NAD 27 in the United States.',/,
     +    ' Transform to NAD 83 first,',
     +    ' using NADCON software.',/)
      endif

        write (luout,*) ' What is its geodetic 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 geodetic 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 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 deflection 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 (3(/), t4, '      ', /,
     +       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, '       ', 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,
     +                   xii, axii, vxii, sdxii, sdxii2, smxii, bgxii,
     +                   eta, aeta, veta, sdeta, sdeta2, smeta, bgeta,
     +                   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 =18)

      double precision xsmall, xbig, ysmall, ybig
      double precision axii, vxii
      double precision sdxii, sdxii2
      real xii, smxii, bgxii
      double precision aeta, veta
      double precision sdeta, sdeta2
      real eta, smeta, bgeta
      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(2*mxarea)

      page = .true.
      ipage = ipage + 1

*******************
* do the statistics
*******************

      if (nconv .ge. 2) then

* calculate mean, variance, standard deviation for deflection
* estimates (arc sec).

*******************
* xi deflection (arc sec)
*******************

        axii = sdxii/dble(nconv)
        vxii = sdxii2/dble(nconv-1) - sdxii**2/dble( nconv*(nconv-1 ))

        if (vxii .gt. 1.0d-6) then
          sdxii = dsqrt(vxii)
        else
          vxii = 0.0d0
          sdxii = 0.0d0
        endif

*******************
* eta deflection (arc sec)
*******************

        aeta = sdeta/dble(nconv)
        veta = sdeta2/dble(nconv-1) - sdeta**2/dble( nconv*(nconv-1 ))

        if (veta .gt. 1.0d-6) then
          sdeta = dsqrt(veta)
        else
          veta = 0.0d0
          sdeta = 0.0d0
        endif
      endif

      if (nconv .lt. 2) then
        axii = xii
        vxii = 0.0d0
        sdxii = 0.0d0
        aeta = eta
        veta = 0.0d0
        sdeta = 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, smxii, bgxii, axii, vxii, sdxii,
     +                   smeta, bgeta, aeta, veta, sdeta, 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, smxii, bgxii, axii, vxii, sdxii,
     +                   smeta, bgeta, aeta, veta, sdeta, 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 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, smxii, bgxii, axii, vxii, sdxii,
     +                   smeta, bgeta, aeta, veta, sdeta, 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) then
      
* itype = 3, itype = 4, and itype = 5 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 deflec.

*     implicit real (a-h, o-z)
*     implicit integer (i-n)
*     implicit undefined (a-z)

      integer mxarea, itype, ios
      parameter (mxarea =18)

      character*1 ans

      integer luin, luout, nout, nin, napar, luarea
      common /inout/ luin, luout, nout, nin, napar, luarea(2*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) format'
      write (luout,*) '  5) 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. 5  .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, both the input latitude and the',
     +                                 ' longitude are expressed in', /,
     +          ' degrees, minutes, and seconds regardless of the',
     +                                           ' method of input.', /)
        write (luout, 140)
  140   format ('                           Computation #:    1',
     +                              '      Region: deflec96SCU(NGS)', //,
     +          '  Station Name:  AAA', //,
     +          '     Latitude           ',
     +                               'Longitude     Xi Deflec.',
     +          '   Eta Deflec.  Hor. Laplace', /,
     +          '                                    ',
     +          '   (arc sec)    (arc sec)    (arc sec)', //,
     +          ' 34 26 39.99984     98 53',
     +          ' 19.99968       -10.90        -5.50         3.77', /)

      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 24 characters',
     +                              ' (columns 33-56) will hold the', /,
     +          ' calculated deflections (Xi, Eta) and the horizontal',
     +                                 ' Laplace correction.  The ', /,
     +          ' rest of the input record (columns 57-80) may contain',
     +                                 ' the station name or be ', /,
     +          ' blank.  The output will be in the same format as the',
     +                                 ' input but will contain ', /,
     +          ' the deflections and Laplace corrections.', /)

        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:', //,
     +             ' DEFLEC Version 3.00 - NGS deflec96 Model', /,
     +          ' <-------- Columns 1-32 -------->',
     +          '<---- Columns 33-56---->',
     +                     '<---- Columns 57-80--->', /,
     +             ' 45 45 45.55555 111 11 11.11111      6.08   -1.40',
     +             '    1.44one', /,
     +             ' 25 55.5555555   76 56.6666666       2.98   22.22',
     +             '  -10.80two', /,
     +             ' 34.444444444    98.888888888      -10.90   -5.50',
     +             '    3.77three', /)

      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 *85* (Defle', /,
     +          'ction) records.  DEFLEC uses information from the',
     +                              ' *80* and *85* 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', /,
     +          ' *85* record is modified by DEFLEC - the *80*',
     +                          ' records and all other records are')
        write (luout, 311)
  311   format (' passed through without change to the output file.',
     +                              '  On the *85* records, the SSN', /,
     +          ' is used to determine which is the corresponding',
     +                                      ' the *80* record.  The', /,
     +          ' information to the right of the SSN (columns',
     +                               ' 41-80) in the *85* record is', /,
     +          ' overwritten by DEFLEC.  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,',
     +                             'Silver Spring, MD 20910.', /)

        write (luout,931)
        read (luin,'(a1)') ans
        write (luout,2)

        write (luout, 330)
  330   format (' The following input example is *80* and *85* records',
     +                                  ' from an ''old'' Blue Book', /,
     +          ' file:', //,
     +          '004560*80*096 KNOXVILLE CT HSE',
     +          '              411906578  N0930548534  W 277  MIA33', /,
     +          '004560*85*096 ABC   bogus deflection              ',
     +          '               10S999  280E999')

        write (luout, 340)
  340   format (/, ' The following example is of the output *85*',
     +                          ' record with the new, interpolated', /,
     +          ' deflections.', //,
     +          '004560*85*096 NGS   program DEFLEC deflc96NCL(NGS)',
     +          '               15N100  325E100',/)

      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 *85* (Defle', /,
     +          'ction) records.  DEFLEC uses information from the',
     +                              ' *80* and *85* 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', /,
     +          ' *85* record is modified by DEFLEC - the *80*',
     +                          ' records and all other records are')
        write (luout, 411)
  411   format (' passed through without change to the output file.',
     +                              '  on the *85* records, the SSN', /,
     +          ' is used to determine which is the corresponding',
     +                                      ' the *80* record.  The', /,
     +          ' information to the right of the SSN (columns',
     +                               ' 41-80) in the *85* record is', /,
     +          ' overwritten by DEFLEC.  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,',
     +                              ' Silver Spring, MD 20910.', /)

        write (luout,931)
        read (luin,'(a1)') ans
        write (luout,2)

        write (luout, 430)
  430   format (' The following input example is *80* and *85* records',
     +                                   ' from a ''new'' Blue Book', /,
     +          ' file:', //,
     +          '004560*80*0096KNOXVILLE CT HSE              ',
     +                        '411906578  N0930548534  W 277  MIA33', /,
     +          '004560*85*0096ABC   bogus deflection              ',
     +          '               10S999  280E999')

        write (luout, 440)
  440   format (/, ' The following example is of the output *85*',
     +                          ' record with the new, interpolated', /,
     +          ' deflection.', //,
     +          '004560*85*0096NGS   program DEFLEC deflc96NCL(NGS)',
     +          '               15N100  325E100',/)

      elseif (itype .eq. 5) then

        write (luout,2)
        write (luout,510)
  510   format (' NGS Horizontal ''Extended'' Blue Book format',
     +                         ' (BIGADJUST format) - *80* (Control', /,
     +          ' Point) and *85* (Deflection) records.  Only the',
     +                                 ' *80* records are used from', /,
     +          ' the input file.  The output consists only of *85*'
     +                               ' records.  The SSN''s for the', /,
     +          ' *85* records are taken from the *80* records - that',
     +                                       ' is, there is an *85*', /,
     +          ' record for each *80* record, regardless of any',
     +                           ' *85* 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'' or ''2'' or ''3'' -',
     +                ' 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 =18)

      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(2*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 deflection',
     +                    ' 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 ) then

* itype=3 is ngs old horizontal blue book
* itype=3 is ngs new horizontal blue book
* itype=5 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 deflec grids 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 deflection 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 '.xii' and '.eta' extensions are 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 =18)
      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(2*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

*** patch for inexact headers (due to 2' spacing)   8-oct-96  dgm

          idx1=idnint(dx1*3600.d0)
          dx(NAREA) = dble(idx1)/3600.d0

          idy1=idnint(dy1*3600.d0)
          dy(NAREA) = dble(idy1)/3600.d0

*****     dx(narea) = dx1     ***** old code
*****     dy(narea) = dy1     ***** old code

          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 =18)

      real vrsion
      character*1 ans

      integer luin, luout, nout, nin, napar, luarea
      common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea)

      write (luout,920)
  920 format (10x, '                   Welcome', /,
     +        10x, '                   to  the', /,
     +        10x, '           National Geodetic Survey', /,
     +        10x, '       Deflection Interpolation Program', //,
     +        10x, '             NGS DEFLEC96 Model', //,
     +        10x, '         (can also interpolate DMEX97,',/,
     +        10x, '            and DCAR97 models)',/)
                                
      write (luout,930) vrsion

  930 format ( /,10x, '               (Version', f5.2, ')', /,
     +       10x,     '                July 16, 1997', /,
     +       10x,     '               Dru Smith, Ph.D.',/,
     +       10x,     '             Dennis Milbert, Ph.D.',/,
     +       10x,     '         Systems Development Divison',/)

      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,
     +                  smxii, bgxii, axii, vxii, sdxii, sdxii2,
     +                  smeta, bgeta, aeta, veta, sdeta, sdeta2,
     +                  xsmall, xbig, ysmall, ybig)

*** this subroutine initializes all the variables needed in deflec

*     implicit real (a-h, o-z)
*     implicit integer (i-n)
*     implicit undefined (a-z)

      integer mxarea
      parameter (mxarea =18)
      character*20 b20
      character*80 b80
      parameter (b20 = '                   ', b80 = b20//b20//b20//b20)

      double precision xsmall, xbig, ysmall, ybig
      double precision axii, vxii
      double precision sdxii, sdxii2
      real smxii, bgxii
      double precision aeta, veta
      double precision sdeta, sdeta2
      real smeta, bgeta
      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(2*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 2*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

      smxii =  1.0e10
      bgxii = -1.0e10

      axii = 0.0d0
      vxii = 0.0d0
      sdxii = 0.0d0
      sdxii2 = 0.0d0

      smeta =  1.0e10
      bgeta = -1.0e10

      aeta = 0.0d0
      veta = 0.0d0
      sdeta = 0.0d0
      sdeta2 = 0.0d0

      xsmall =  1.0d10
      xbig   = -1.0d10
      ysmall =  1.0d10
      ybig   = -1.0d10

      return
      end

      subroutine interp (nogo, resp, xpt, ypt, xii, eta, 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 =18)
      character*20 b20
      character*80 b80
      parameter (b20 = '                   ', b80 = b20//b20//b20//b20)

      double precision xpt, ypt
      real buf, xii, eta
      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, ifile1, ifile2, 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(2*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)

      ifile1= luarea(2*iarea-1)
      ifile2= luarea(2*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 deflection files !!!

****************
* xii deflection
****************

*** lower boundary

                read (ifile1,rec=jy+1) idum, (buf(j), j=1,nc0)
                tee1 = buf(ix)
                tee2 = buf(ix+1)
                tee3 = buf(ix+2)

*** middle row

                read (ifile1,rec=jy1+1) idum, (buf(j), j=1,nc0)
                tee4 = buf(ix)
                tee5 = buf(ix+1)
                tee6 = buf(ix+2)

*** upper boundary

                read (ifile1,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)

        xii = qterp2(y, fx0, fx1, fx2)

****************
* eta deflection
****************

*** lower boundary

                read (ifile2,rec=jy+1) idum, (buf(j), j=1,nc0)
                tee1 = buf(ix)
                tee2 = buf(ix+1)
                tee3 = buf(ix+2)

*** middle row

                read (ifile2,rec=jy1+1) idum, (buf(j), j=1,nc0)
                tee4 = buf(ix)
                tee5 = buf(ix+1)
                tee6 = buf(ix+2)

*** upper boundary

                read (ifile2,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)

        eta = 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 deflec.

*     implicit real (a-h, o-z)
*     implicit integer (i-n)
*     implicit undefined (a-z)

      integer mxarea
      parameter (mxarea =18)
      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(2*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 ''deflec.out''.'
      read (luin,'(a80)') ofile
      if (ofile .eq. b80) ofile = 'deflec.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'
        write (luout,*) '  5) 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. 5  .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,
     +                  xii, sdxii, sdxii2, smxii, bgxii,
     +                  eta, sdeta, sdeta2, smeta, bgeta,
     +                  xsmall, xbig, ysmall, ybig,
     +                  page, screen, name)

**********************************************************************
* this subroutine loops through the input data (either an input data *
* file or interactively), calculates the deflection 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 =18)

      double precision xsmall, xbig, ysmall, ybig, xpt, ypt
      real vrsion, sla, slo
      double precision sdxii, sdxii2
      real  xii, smxii, bgxii, xii2
      double precision sdeta, sdeta2
      real  eta, smeta, bgeta, eta2
      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(2*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, xii, eta, 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
* arc sec....

        xii2   = ( dble(xii) )**2
        sdxii2 = sdxii2 + xii2
        sdxii  = sdxii  + dble(xii)

        eta2   = ( dble(eta) )**2
        sdeta2 = sdeta2 + eta2
        sdeta  = sdeta  + dble(eta)

* then the ranges

        xsmall = dmin1( xsmall, xpt)
        xbig   = dmax1( xbig  , xpt)
        ysmall = dmin1( ysmall, ypt)
        ybig   = dmax1( ybig  , ypt)

        smxii = min( smxii, xii)
        bgxii = max( bgxii, xii)

        smeta = min( smeta, eta)
        bgeta = max( bgeta, eta)

**********************************
** write to output file and screen
**********************************

        call wrtpt (itype, nconv, vrsion, name,
     +              idla, imla, sla, idlo, imlo, slo,
     +              xii, eta, 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 deflec grids.
* two file are 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 =18)

      character*1 ans
      logical nodata

      integer luin, luout, nout, nin, napar, luarea
      common /inout/ luin, luout, nout, nin, napar, luarea(2*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
      
      
*** dmex97 and dcar97 flags
      logical ldmex,ldcar
      common/scienc/ldmex,ldcar
      
* initialize

      nodata = .false.
      narea = 0
      write (luout,100)
  100 format (' DEFLEC is now opening the files containing the',
     +        ' gridded data.') 

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

* NOTE: ALL INDICATIONS ARE THAT DGRIDS IS DISABLED,
* thus, extracted it, april-8-1997, das
      
*** added support for models DMEX97, DCAR97
      if(ldmex) write(luout,*) ' **************************************'
      if(ldmex) write(luout,*) ' DMEX97 models encountered ************'
      if(ldmex) write(luout,*) ' **************************************'
      if(ldmex) write(luout,*) ' '  
     
      if(ldcar) write(luout,*) ' **************************************'
      if(ldcar) write(luout,*) ' DCAR97 models encountered ************'
      if(ldcar) write(luout,*) ' **************************************'
      if(ldcar) write(luout,*) ' '  
 
      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 two data files

      implicit real (a-h, o-z)
      implicit integer (i-n)
*     implicit undefined (a-z)

      integer mxarea
      parameter (mxarea =18)
      character*23 b23
      parameter (b23 = '                       ')
      character*69 b69
      parameter (b69 = b23//b23//b23)

      real dx1, dy1, dx2, dy2
      real xmax1, xmin1, ymax1, ymin1
      real x01, y01, angle1, x02, y02, angle2
      integer iflag1, iflag2, n1, n2, n3, n4
      integer itemp, lrecl, ixii, ieta, ifile, ios
      integer nc1, nr1, nz1, nc2, nr2, nz2
      character*80 card
      character*69 axii, aeta
      character*65 afile
      character*56 rident
      character*8 pgm
      logical gflag, nogo, oflag, eflag1, eflag2

      integer luin, luout, nout, nin, napar, luarea
      common /inout/ luin, luout, nout, nin, napar, luarea(2*mxarea)
      
      logical ldmex,ldcar
      common/scienc/ldmex,ldcar
      
      data iflag1 /1/, iflag2 /2/
      data oflag /.false./, eflag1 /.false./, eflag2 /.false./

* initialize

      nogo = .false.

* form complete names of grid files

      call nblank (afile, iflag2, n2)
      if (n2 .eq. 0) stop 'Logical coding error in OPENF'
      axii = b69
      axii(1:n2) = afile
      axii(n2+1:n2+4) = '.xii'
      aeta = b69
      aeta(1:n2) = afile
      aeta(n2+1:n2+4) = '.eta'

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

* xi deflection (seconds) grid file

      lrecl = 256
      ixii = 2*itemp - 1
      ifile = ixii + 10
      luarea(ixii) = ifile
      inquire (file=axii, exist=eflag1, opened=oflag)
      if (.not. eflag1) goto 100
      if (oflag) goto 980
      open (ifile,file=axii,form='unformatted',status='old',
     +       access='direct',recl=lrecl,err=940,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=axii,form='unformatted',status='old',
     +       access='direct',recl=lrecl,err=940,iostat=ios)

* eta deflection (seconds) grid file

  100 lrecl = 256
      ieta = 2*itemp
      ifile = ieta + 10
      luarea(ieta) = ifile
      inquire (file=aeta, exist=eflag2, opened=oflag)
      if (.not. eflag1) goto 910
      if (.not. eflag2) goto 920
      if (oflag) goto 980
      open (ifile,file=aeta,form='unformatted',status='old',
     +       access='direct',recl=lrecl,err=940,iostat=ios)
      read (ifile,rec=1) rident, pgm, nc2, nr2, nz2, x02, dx2,
     +                   y02, dy2, angle2
      close (ifile)

      lrecl = 4*(nc2+1)
      open (ifile,file=aeta,form='unformatted',status='old',
     +       access='direct',recl=lrecl,err=940,iostat=ios)

* check to see if the two files have the same variables

      if ( (nc2 .ne. nc1)  .or.  (nr2 .ne. nr1)  .or.
     +     (nz2 .ne. nz1)  .or.
     +     (x02 .ne. x01)  .or.  (dx2 .ne. dx1)  .or.
     +     (y02 .ne. y01)  .or.  (dy2 .ne. dy1)  .or.
     +     (angle2 .ne. angle1) ) goto 960
      
      
c set flag if a dmex97 or dcar97 model is encountered, das
c only do this once the 2 grid files have been checked as matching
c this means we don't need to check "aeta", just "axii"
      
      if(axii(1:4).eq.'dmex' .or. axii(1:4).eq.'DMEX')ldmex=.true.
      if(axii(1:4).eq.'dcar' .or. axii(1:4).eq.'DCAR')ldcar=.true.
      
* 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 files do not exist

  910 continue
      nogo = .true.
      if (gflag) then
        call nblank (axii, iflag1, n1)
        call nblank (axii, iflag2, n2)
        call nblank (card, iflag1, n3)
        call nblank (card, iflag2, n4)
        if (eflag2) then

* .eta exists, .xii does not exist

          write (luout, 925) axii(n1:n2), card(n3:n4), aeta(n1:n2)
        else

* neither .eta nor .xii exist

          write (luout, 915) axii(n1:n2), aeta(n1:n2), card(n3:n4)
  915 FORMAT(' Files "',a,',',a,'" in record "',a,'" do not exist')
        endif
      endif
      close ( luarea(ixii) )
      close ( luarea(ieta) )
      goto 9999

  920 continue
      nogo = .true.
      if (gflag) then

* .xii exists, .eta does not exist

        call nblank (axii, iflag1, n1)
        call nblank (axii, iflag2, n2)
        call nblank (card, iflag1, n3)
        call nblank (card, iflag2, n4)
        write (luout, 925) aeta(n1:n2), card(n3:n4), axii(n1:n2)
  925   format (/, 5x, ' *********** WARNING ***********', /,
     +          5x, ' The grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' from record:', /,
     +          6x, '''', a, '''', /,
     +          5x, ' does not exist!  However, the grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' does exist!', /,
     +          5x, ' *******************************', /)
      endif
      close ( luarea(ixii) )
      close ( luarea(ieta) )
      goto 9999

* grid files already open

  940 continue
      nogo = .true.
      if (gflag) then
        call nblank (axii, iflag1, n1)
        call nblank (axii, iflag2, n2)
        call nblank (card, iflag1, n3)
        call nblank (card, iflag2, n4)
        write (luout, 945) axii(n1:n2), aeta(n1:n2), card(n3:n4), ios
  945   format (/, 5x, ' *********** WARNING ***********', /,
     +          5x, ' The grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' and/or grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' from record:', /,
     +          6x, '''', a, '''', /,
     +          5x, ' cannot be opened!',
     +                      '  These files will be ignored.', 5x, i5, /,
     +          5x, ' *******************************', /)
      endif
      close ( luarea(ixii) )
      close ( luarea(ieta) )
      goto 9999

* grid files do not agree

  960 continue
      nogo = .true.
      call nblank (axii, iflag1, n1)
      call nblank (axii, iflag2, n2)
      write (luout, 965) axii(n1:n2), aeta(n1:n2)
  965 format (/, 5x, ' *********** ERROR ***********', /,
     +        5x, ' The header information in grid files', /,
     +        6x, '''', a, '''', /,
     +        5x, ' and', /,
     +        6x, '''', a, '''', /,
     +        5x, ' do not agree!  One or both of these files must',
     +                                             ' be corrupted.', /,
     +        5x, ' These files will be ignored.', /,
     +        5x, ' *****************************', /)
      close ( luarea(ixii) )
      close ( luarea(ieta) )
      goto 9999

* grid files already open

  980 continue
      nogo = .true.
      if (gflag) then
        call nblank (axii, iflag1, n1)
        call nblank (axii, iflag2, n2)
        write (luout, 985) axii(n1:n2), aeta(n1:n2)
  985   format (/, 5x, ' *********** ERROR ***********', /,
     +          5x, ' The grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' and the grid file', /,
     +          6x, '''', a, '''', /,
     +          5x, ' have already been opened!  These files',
     +                                      ' will not be reopened.', /,
     +          5x, ' *****************************', /)
      endif
      goto 9999
      end

      subroutine print (lu, nconv, name, vrsion, idla, imla, sla,
     +                  idlo, imlo, slo, xii, eta, resp, ipage, page)

* this subroutine prints out the 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 rad
      parameter(rad=57.29577951)

      real vrsion
      real sla, slo, xii, eta
      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,6)
        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, '         Deflection Computation')
    6 format (20x, ' ')
   26 format (20x, '  DEFLEC96/DMEX97/DCAR97 Interpolation'   )
    7 format (20x, '             DEFLEC 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)
      hlap=-tan( (idla+imla/60.0+sla/3600.0)/rad )*eta
      write (lu,923) idla, imla, sla, idlo, imlo, slo, xii, eta, hlap

  921 format (//, 27x, 'Computation #: ', i4, 8x, 'Region: ', a15, /)
  922 format (2x, 'Station Name:  ', a80, /)
  900 format ( 5x, 'Latitude', 11x, 'Longitude',
     +         4x, ' Xi Deflec.', 3x, 'Eta Deflec.',
     +         2x, 'Hor. Laplace')
  927 format (39x, '(arc sec)    (arc sec)    (arc sec)', /)
  923 format ( 1x, i2, 1x, i2.2, f9.5, 4x, i3, 1x, i2.2, f9.5,
     +        3(4x, f9.2), /)

      return
      end

      subroutine print2 (lu, nconv, vrsion, idla, imla, sla, xii, eta)

* this subroutine prints out the actual transformation 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 rad
      parameter(rad=57.29577951)

      real vrsion, xii, eta
      real sla
      integer idla, imla
      integer lu, nconv

      character*80 card
      common /curnt/ card

* write header record to identify source of deflection and value

      if (nconv .eq. 1) then
        write (lu, 10) vrsion
   10   format ('DEFLEC Version', f5.2,
     +          ' - NGS DEFLEC96/DMEX97/DCAR97 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 xi deflection
* overwrites whatever is in columns 33-40

      hlap=-tan( (idla+imla/60.0+sla/3600.0)/rad )*eta
      write (lu,100) card(1:32), xii, eta, hlap, card(41:64)
  100 format (a, 3f8.2, a)

      return
      end

      subroutine print3 (itype, lu, xii, eta, resp)

* this subroutine prints out the actual transformation results into
* a *85* 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 xii,xii2
      integer ixii
      real eta,eta2
      integer ieta
      integer lu, itype
      character*125 crd125
      character*80 crd80
      character*15 resp
      character*3 astd
      character*1 ans, aew

      character*80 card
      common /curnt/ card

c - dmex97/dcar97 flags
      logical ldmex,ldcar
      common/scienc/ldmex,ldcar

*** standard deviation is assumed to be 1.0 arc sec

      data astd /'100'/

* 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 deflection
* overwrites whatever is in columns 33-40

      if (itype .eq. 3) then

* old blue book

        crd80=' '
        crd80(1:7) = card(1:7)
        crd80(8:9) = '85'
        crd80(10:14) = card(10:14)
        crd80(15:20) = 'NGS  '
        crd80(21:62) = 'program DEFLEC '//resp
        call asgnns(xii,xii2,ans)
        ixii = nint(xii2*100.0)
        write (crd80(63:67),'(i5)') ixii
        crd80(68:68) = ans
        crd80(69:71) = astd
        call asgnew(eta,eta2,aew)
        ieta = nint(eta2*100.0)
        write (crd80(72:76),'(i5)') ieta
        crd80(77:77) = aew
        crd80(78:80) = astd
        write (lu,'(a80)') crd80

      elseif (itype .eq. 4) then

* new blue book

c - NOTE, as of april 8, 1997, NGS blue book record *85*
c - has NO place for the deflection model used in the
c - compuation of the deflection.  Thus, until that
c - policy is changed, the ldmex and ldcar variables
c - need not be used

c - when changed, this line will assign the code for DEFLEC96
c - when changed, this line will assign the code for DMEX97
c - when changed, this line will assign the code for DCAR97
c - SEE program geoid.for (version 3.1 or later) for details

        crd80=' '
        crd80(1:7) = card(1:7)
        crd80(8:9) = '85'
        crd80(10:14) = card(10:14)
        crd80(15:20) = 'NGS  '
        crd80(21:62) = 'program DEFLEC '//resp
        call asgnns(xii,xii2,ans)
        ixii = nint(xii2*100.0)
        write (crd80(63:67),'(i5)') ixii
        crd80(68:68) = ans
        crd80(69:71) = astd
        call asgnew(eta,eta2,aew)
        ieta = nint(eta2*100.0)
        write (crd80(72:76),'(i5)') ieta
        crd80(77:77) = aew
        crd80(78:80) = astd
        write (lu,'(a80)') crd80

      elseif (itype .eq. 5) then

* extended blue book

        crd125=' '
        crd125(1:7) = card(1:7)
        crd125(8:9) = '85'
        crd125(10:14) = card(10:14)
        crd125(15:20) = 'NGS  '
        crd125(21:62) = 'program DEFLEC '//resp
        call asgnns(xii,xii2,ans)
        ixii = nint(xii2*100.0)
        write (crd125(63:67),'(i5)') ixii
        crd125(68:68) = ans
        crd125(69:71) = astd
        call asgnew(eta,eta2,aew)
        ieta = nint(eta2*100.0)
        write (crd125(72:76),'(i5)') ieta
        crd125(77:77) = aew
        crd125(78:80) = astd
        crd125(81:125) = ' '
        write (lu,'(a125)') crd125
      endif

      return
      end

      subroutine asgnns(val,val2,ans)

*** strip minus and assign positive north, negative south

      real val,val2
      character*1 ans

      if(val.ge.0.0) then
        val2= val
        ans='N'
      else
        val2=-val
        ans='S'
      endif

      return
      end

      subroutine asgnew(val,val2,aew)

*** strip minus and assign positive east, negative west

      real val,val2
      character*1 aew

*** 'E' means that astro is "east" of the geodetic
*** gridded values are stored positive east

      if(val.ge.0.0) then
        val2= val
        aew='E'

*** 'W' means that astro is "west" of the geodetic

      else
        val2=-val
        aew='W'
      endif

      return
      end

      subroutine report (lu, smxii, bgxii, axii, vxii, sdxii,
     +                       smeta, bgeta, aeta, veta, sdeta,
     +                   ipage, page)

* this subroutine prints out the statistics for the transformations

*     implicit real (a-h, o-z)
*     implicit integer (i-n)
*     implicit undefined (a-z)

      double precision axii, vxii, sdxii
      real smxii, bgxii
      double precision aeta, veta, sdeta
      real smeta, bgeta
      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 (29x, 'Deflection Computation')
        write (lu,100)
  100   format (/, 30x, 'Statistics for Region', /, 1x, 79('='), /)
        write (lu,900)
  900   format (38x, 'Xi                    Eta')
        write (lu,910)
  910   format ('  Deflections          ',
     +         10x, 'Min', 6x, 'Max', 10x, 'Min', 6x, 'Max' )

        write (lu,110) smxii, bgxii,smeta, bgeta
  110   format ('  Range of values (arc sec) ', 2f9.2, 4x,2f9.2)
        write (lu,130) axii,aeta
  130   format (/, 4x, 'Mean deflections       ',5x,f9.2,14x,f9.2)
        write (lu,131) vxii,veta
  131   format (   4x, 'Variance of deflections',5x,f9.2,14x,f9.2)
        write (lu,132) sdxii,sdeta
  132   format (   4x, 'Std. dev. of deflections',4x,f9.2,14x,f9.2)

      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 =18)
      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(2*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 deflection
* 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 =18)
      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(2*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 ('new' blue book) or type 5 (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 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 =18)

      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(2*mxarea)

      data iflag1 /1/, iflag2 /2/

***********************************
* for input file of itype = 3
***********************************

      read (nin,6000,end=9999) card
 6000 format (a80)

      if (card(8:9) .ne. '85'  .and.  itype .ne. 5) then

* for itype = 3 and itype 4 copy input non-*80*, non-*85* 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,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,
     +                  xii, eta, 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 =18)

      real vrsion, xii, eta
      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(2*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, xii, eta, 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, xii, eta, resp, ipage, page)

        elseif (itype .eq. 2) then

**************************
* for free format type 2
**************************

          call print2 (nout, nconv, vrsion, idla, imla, sla, xii, eta)

        elseif (itype .eq. 3  .or.  itype .eq. 4  .or.
     +          itype .eq. 5) then

**********************************************************
* for itype = 3 the horizontal blue book old specification
* for itype = 4 the horizontal blue book new specification
* for itype = 4 - the extended blue book specification
**********************************************************

          call print3 (itype, nout, xii, eta, resp)

        endif

*******************
* for screen output
*******************

        if (screen) then
          call print (luout, nconv, name, vrsion, idla, imla, sla,
     +                idlo, imlo, slo, xii, eta, 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