	  PROGRAM NADCON

***********************************************************************
*                                                                     *
* PROGRAM :   NADCON                                                  *
*                                                                     *
* PURPOSE:    COMPUTATION PROGRAM TO CONVERT (OR TRANSFORM)           *
*             POSITIONAL DATA (E.G., LATITUDES AND LONGITUDES) FROM   *
*             THE NORTH AMERICAN DATUM OF 1927 (NAD 27) TO THE        *
*             NORTH AMERICAN DATUM OF 1983 (NAD 83).  THIS PROGRAM    *
*             CAN COMPUTE FROM FROM EITHER DATUM TO THE OTHER.        *
*                                                                     *
*             THE ACTUAL COMPUTATION IS PERFORMED AS AN INTERPOLATION *
*             FROM A REGULARLY-SPACED GRID OF POINTS OBTAINED FROM THE*
*             FITTING OF A MINIMUM-CURVATURE SURFACE TO THE ACTUAL    *
*             SHIFT DATA RESULTING FROM THE NAD 83 ADJUSTMENT.        *
*                                                                     *
*             THE INTERPOLATION IS ACCOMPLISHED BY LOCALLY FITTING    *
*             A CURVED POLYNOMIAL SURFACE TO THE FOUR DATA POINTS     *
*             DEFINING THE SQUARE WHICH SURROUND THE (X,Y) PAIR       *
*             WHERE THE INTERPOLATION IS TO TAKE PLACE.               *
*                                                                     *
*             THE POLYNOMIAL SURFACE IS DEFINED BY:                   *
*                                                                     *
*                         A+BX+CY+DXY=Z                               *
*                                                                     *
*             THE PROGRAM REQUIRES THAT THE USER SPECIFY:             *
*                                                                     *
*             1)  THE NAME OF AN OUTPUT FILE                          *
*                                                                     *
*             2)  THE NAME OF AN INPUT FILE (IF AVAILABLE).           *
*                                                                     *
*                                                                     *
*                                                                     *
*             ESTIMATES OF DATUM SHIFTS IN TERMS OF METERS ARE        *
*             COMPUTED FROM THE SHIFT ESTIMATES USING ELLIPSOIDAL     *
*             SCALING.                                                *
*                                                                     *
*             THIS PROGRAM ALLOWS FOR EITHER NGS STANDARD HORIZONTAL  *
*             DATA FORMATS AS SPECIFIED IN THE FGCC PUBLICATION,      *
*             COMMONLY KNOWN AS THE 'HORIZONTAL BLUE BOOK' (SEE       *
*             SUBROUTINE TYPE3), 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:  2.00                                                 *
*                                                                     *
* VERSION DATE:  JULY 1, 1991                                         *
*                                                                     *
*        AUTHOR:   WARREN T. DEWHURST, PH.D.                          *
*                    LIEUTENANT COMMANDER, NOAA                       *
*                  ALICE R. DREW                                      *
*                    SENIOR GEODESIST, HORIZONTAL NETWORK BRANCH      *
*                  NATIONAL GEODETIC SURVEY, NOS, NOAA                *
*                  ROCKVILLE, MD   20852                              *
***********************************************************************

***********************************************************************
*                                                                     *
*                  DISCLAIMER                                         *
*                                                                     *
*   THIS PROGRAM AND SUPPORTING INFORMATION IS FURNISHED BY THE       *
* GOVERNMENT OF THE UNITED STATES OF AMERICA, AND IS ACCEPTED AND     *
* USED BY THE RECIPIENT WITH THE UNDERSTANDING THAT THE UNITED STATES *
* GOVERNMENT MAKES NO WARRANTIES, EXPRESS OR IMPLIED, CONCERNING THE  *
* ACCURACY, COMPLETENESS, RELIABILITY, OR SUITABILITY OF THIS         *
* PROGRAM, OF ITS CONSTITUENT PARTS, OR OF ANY SUPPORTING DATA.       *
*                                                                     *
*   THE GOVERNMENT OF THE UNITED STATES OF AMERICA SHALL BE UNDER NO  *
* LIABILITY WHATSOEVER RESULTING FROM ANY USE OF THIS PROGRAM.  THIS  *
* PROGRAM SHOULD NOT BE RELIED UPON AS THE SOLE BASIS FOR SOLVING A   *
* PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO PERSON   *
* OR PROPERTY.                                                        *
*                                                                     *
*   THIS PROGRAM IS PROPERTY OF THE GOVERNMENT OF THE UNITED STATES   *
* OF AMERICA.  THEREFORE, THE RECIPIENT FURTHER AGREES NOT TO ASSERT  *
* PROPRIETARY RIGHTS THEREIN AND NOT TO REPRESENT THIS PROGRAM TO     *
* ANYONE AS BEING OTHER THAN A GOVERNMENT PROGRAM.                    *
*                                                                     *
***********************************************************************
C
C THE SUBROUTINES HAVE BEEN ALTERED AND REORDERED FOR
C THE INTERGRATION WITH THE CORPSCON 3.0 PROGRAM.....
C
C LAST MODIFICATION:
C       2 DECEMBER 1991
C       Parameter  Change 2 June 92
C       Parameter Change 24 February 94 - changed 
C            corpsd(22)->corpsd(24)
C MODIFIED BY:
C       JEFF WALKER
C       U.S. ARMY TOPOGRAPHIC ENGINEERING CENTER
C       FORT BELVOIR, VA 22060-5546
C       (703) 355-2766
C
C

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

	  DOUBLE PRECISION VRSION
	  INTEGER MXAREA
	  PARAMETER (VRSION = 2.00D0, MXAREA = 7)

	  DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS
	  INTEGER IFILE
	  LOGICAL NODATA, NOGO

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

	  INTEGER LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA
	  COMMON /INOUT/ LUIN, LUOUT, NOUT, NIN, NAPAR, LUAREA(2*MXAREA)
	  DOUBLE PRECISION CORPSD(29)
	  INTEGER*2 FLAG, KEY
	  CHARACTER*25 NAME
	  CHARACTER*60 INFILE

	  WRITE(*,*) 'STARTING CORPSNAD'

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

	  CALL INITL

*******************************************************
* OPEN NADCON DATA FILES (LATITUDE AND LONGITUDE GRIDS)
*******************************************************

	  CALL NGRIDS (NODATA)
	  IF (NODATA) GOTO 9999

C********************************
C READ DATA FROM CORPSCON PROGRAM
C********************************

	  READ(5,10) INFILE
   10 FORMAT(A60)
	  READ(5,15) KEY
   15 FORMAT(I1)
	WRITE(*,*) 'CORPSNAD INPUT PARAMETERS READ SUCCESSFULLY'

C********************************
C OPEN BINARY INPUT/OUTPUT FILES
C FROM THE CORPSCON PROGRAM
C********************************

	  OPEN(91,FILE=INFILE,ACCESS='SEQUENTIAL',FORM='BINARY',
     &ERR=25,STATUS='OLD')
	WRITE(*,*) 'CORPSNAD INPUT FILE OPENED SUCCESSFULLY'
	  OPEN(93,FILE='INT$OUT',ACCESS='SEQUENTIAL',FORM='BINARY',
     &ERR=25,STATUS='UNKNOWN')
	WRITE(*,*) 'CORPSNAD OUTPUT FILE OPENED SUCCESSFULLY'

*********************************
* LOOP (ONCE FOR EACH CONVERSION)
C UNTIL END OF BINARY FILE (91)
*********************************

  20  READ (91,END=25) NAME,(CORPSD(II),II=1,29),FLAG
	  IF(KEY.EQ.1) THEN
	YPT = CORPSD(9)
	XPT = -CORPSD(10)
	  ELSE
	YPT = CORPSD(11)
	XPT = -CORPSD(12)
	  ENDIF

************************
* DO THE TRANSFORMATION
************************

	  NOGO = .FALSE.
	  CALL TRANSF (NOGO, XPT, YPT, XPT2, YPT2,
     +               DLAM, DLOM, DLAS, DLOS, KEY)

	  IF (NOGO) THEN
	FLAG = IOR(FLAG,4)
	IF(KEY.EQ.1) THEN
		CORPSD(11) = CORPSD(9)
		CORPSD(12) = CORPSD(10)
	ELSE
		CORPSD(9) = CORPSD(11)
		CORPSD(10) = CORPSD(12)
	ENDIF
	CORPSD(21) = 0
	CORPSD(22) = 0
	  ELSE
c       DISTM = SQRT(DLAM**2 + DLOM**2)
	IF(KEY.EQ.1) THEN
		CORPSD(11) = YPT2
		CORPSD(12) = -XPT2
	ELSE
		CORPSD(9) = YPT
		CORPSD(10) = -XPT
	ENDIF
	CORPSD(21) = DLAM
	CORPSD(22) = DLOM
	ENDIF
	WRITE (93) NAME,(CORPSD(II),II=1,29),FLAG
	GOTO 20

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

  25  DO 1010 IFILE = 1, 2*NAREA
	CLOSE (LUAREA(IFILE), STATUS='KEEP')
 1010 CONTINUE
C      WRITE(LUOUT,*) 'ERROR OPENING FILES'
	  CLOSE(91)
	  CLOSE(93)
	WRITE(*,*) 'CORPSNAD INPUT/OUTPUT FILES CLOSED'
	WRITE(*,*) 'CORPSNAD COMPLETE'
 9999 STOP
	  END

	  SUBROUTINE COEFF (TEE1, TEE2, TEE3, TEE4, AY, BEE, CEE, DEE)

**********************************************************************
** SUBROUTINE COEFF: GENERATES COEFFICIENTS FOR SURFACE FUNCTION     *
**********************************************************************

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

	  DOUBLE PRECISION AY, BEE, CEE, DEE
	  DOUBLE PRECISION TEE1, TEE2, TEE3, TEE4

	  AY = TEE1
	  BEE = TEE3 - TEE1
	  CEE = TEE2 - TEE1
	  DEE = TEE4 - TEE3 - TEE2 + TEE1

	  RETURN
	  END

**********************************************************************
	  SUBROUTINE DGRIDS
**********************************************************************
* This subroutine opens the NADCON grids using the default grid
* names and locations.  The default names of the grid areas are
* given in DAREAS and the default base file locations are in DFILES

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

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

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

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

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

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

	  DATA DUM / B80 /

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

	  DATA DFILES /'conus', 'hawaii', 'prvi',
     +             'stlrnc', 'stgeorge', 'stpaul', 'alaska'/

	  DATA DAREAS /'Conus', 'Hawaii', 'P.R. and V.I.',
     +             'St. Lawrence I.', 'St. George I.', 'St. Paul I.'
     +             ,'Alaska'/

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

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

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

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

	IF (.NOT. NOGO) THEN

* Set of files opened OK and variables read

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

	  WRITE (LUOUT,120) NAREA,AREAS(NAREA),YMIN1,XMIN1,YMAX1,XMAX1
  120     FORMAT (2X, I2, 2X, A15, 4(2X, F10.5))

	ENDIF
  140 CONTINUE

  999 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 INTRP.     *
**********************************************************************
* 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 DOUBLE PRECISION (A-H, O-Z)
	  IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

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

	  NOGO = .FALSE.

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

	  IF (XPT .GE. XMAX  .OR.  XPT .LE. XMIN   .OR.
     +    YPT .GE. YMAX  .OR.  YPT .LE. YMIN ) THEN
		NOGO = .TRUE.
c   WRITE (LUOUT,*) '***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 = ( XPT - XMIN )/DX + 1.D0
	  YGRID = ( YPT - YMIN )/DY + 1.D0

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

	  IROW = IDINT(YGRID)
	  JCOL = IDINT(XGRID)

  200 RETURN
	  END


**********************************************************************
	  SUBROUTINE INITL
**********************************************************************
*** This subroutine initializes all the variables needed in NADCON

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

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

	  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

	  RETURN
	  END

**********************************************************************
	  SUBROUTINE INTRP (IAREA, IROW, NC, JCOL, XGRID, YGRID,
     +                  XPT, YPT, XPT2, YPT2, DLOS, DLAS, DLAM, DLOM)
**********************************************************************
** DETERMINE SURFACE FUNCTION FOR THIS GRID SQUARE                   *
** AND INTERPOLATE A VALUE, ZEE, FOR XPT, YPT                        *
**********************************************************************

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

	  INTEGER MXAREA, MAXCOL
	  PARAMETER (MAXCOL = 600, MXAREA = 7)

	  DOUBLE PRECISION XPT, YPT, XPT2, YPT2, XGRID, YGRID
	  DOUBLE PRECISION DLOS, DLAS, DLAM, DLOM
	  DOUBLE PRECISION TEE1, TEE2, TEE3, TEE4, ZEE
	  INTEGER IROW, JCOL, NC, IAREA, IFILE, IDUM, J
	  REAL BUF(MAXCOL)

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

	  DOUBLE PRECISION AY1, BEE1, CEE1, DEE1, AY2, BEE2, CEE2, DEE2
	  SAVE AY1, BEE1, CEE1, DEE1, AY2, BEE2, CEE2, DEE2
	  INTEGER IROWL, JCOLL, IAREAL
	  SAVE IROWL, JCOLL, IAREAL

	  DATA IROWL / 0 /, JCOLL / 0 /, IAREAL / 0 /

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

	  IF ( IROW .NE. IROWL  .OR.  JCOL .NE. JCOLL  .OR.
     +    IAREA .NE. IAREAL ) THEN

* Lower boundary

	IFILE = LUAREA( 2*IAREA - 1 )
	READ (IFILE,REC=IROW+1) IDUM, (BUF(J), J=1,NC)
	TEE1 = DBLE( BUF(JCOL) )
*       TEE4 = DBLE( BUF(JCOL+1) )
	TEE3 = DBLE( BUF(JCOL+1) )

* Upper boundary

	READ (IFILE,REC=IROW+2) IDUM, (BUF(J), J=1,NC)
	TEE2 = DBLE( BUF(JCOL) )
*       TEE3 = DBLE( BUF(JCOL+1) )
	TEE4 = DBLE( BUF(JCOL+1) )

	CALL COEFF (TEE1, TEE2, TEE3, TEE4, AY1, BEE1, CEE1, DEE1)

	  ENDIF

	  CALL SURF (XGRID, YGRID, ZEE, AY1, BEE1, CEE1, DEE1, IROW, JCOL)
	  DLAS = ZEE

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

	  IF ( IROW .NE. IROWL  .OR.  JCOL .NE. JCOLL  .OR.
     +    IAREA .NE. IAREAL ) THEN


* Lower boundary

	IFILE = LUAREA( 2*IAREA )
	READ (IFILE,REC=IROW+1) IDUM, (BUF(J), J=1,NC)
	TEE1 = DBLE( BUF(JCOL) )
*       TEE4 = DBLE( BUF(JCOL+1) )
	TEE3 = DBLE( BUF(JCOL+1) )

* Upper boundary

	READ (IFILE,REC=IROW+2) IDUM, (BUF(J), J=1,NC)
	TEE2 = DBLE( BUF(JCOL) )
*       TEE3 = DBLE( BUF(JCOL+1) )
	TEE4 = DBLE( BUF(JCOL+1) )

	CALL COEFF (TEE1, TEE2, TEE3, TEE4, AY2, BEE2, CEE2, DEE2)

	  ENDIF

	  CALL SURF (XGRID, YGRID, ZEE, AY2, BEE2, CEE2, DEE2, IROW, JCOL)
	  DLOS = ZEE

**************************
* COMPUTE THE NAD 83 VALUES
**************************

	  YPT2 = YPT + DLAS/3600.D0

* Longitude is positive west in this subroutine

	  XPT2 = XPT - DLOS/3600.D0

*********************************************************************
* USE THE NEW ELLIPSOIDAL VARIABLES TO COMPUTE THE SHIFTS IN METERS
*********************************************************************

	  CALL METERS (YPT, XPT, YPT2, XPT2, DLAM, DLOM)

* Update the last-value variables

	  IROWL = IROW
	  JCOLL = JCOL
	  IAREAL = IAREA

	  RETURN
	  END

**********************************************************************
	  SUBROUTINE METERS (LAT1, LONG1, LAT2, LONG2, LATMTR, LONMTR)
**********************************************************************
* This subroutine computes the difference in two positions in meters.
*
* This method utilizes ellipsoidal rather than spherical
* parameters.  I believe that the original approach and code
* for this came from Ed McKay.
* The reference used by Ed McKay for this was:
*       'A Course in Higher Geodesy' by P.W. Zakatov, Israel Program
*       for Scientific Translations, Jerusalem, 1962
*
*       Warren T. Dewhurst
*       11/1/89
* Note that this subroutine is set up for +west longitude

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

* I think that these are GRS80 parameters

	  DOUBLE PRECISION AXIS, E2, RHOSEC
	  PARAMETER (AXIS = 6378137.0D0)
	  PARAMETER (E2 = 0.0066943800229D0)
	  PARAMETER (RHOSEC = 206264.806247D0)

	  DOUBLE PRECISION W, LM, LP, AVLAT
	  DOUBLE PRECISION LAT1S, LAT2S, LONG1S, LONG2S, LAT1, LAT2
	  DOUBLE PRECISION LONG1, LONG2, DLAT, DLONG
	  DOUBLE PRECISION LATMTR, LONMTR


*     LAT1  = (LATSEC + 60.D0*( LATMIN + 60.D0*LATDEG) )/RHOSEC
*     LONG1 = (LONSEC + 60.D0*( LONMIN + 60.D0*LONDEG) )/RHOSEC
*     LAT2  = (LATSEC + 60.D0*( LATMIN + 60.D0*LATDEG) )/RHOSEC
*     LONG2 = (LONSEC + 60.D0*( LONMIN + 60.D0*LONDEG) )/RHOSEC

* Change into sec.ddd and convert to +west longitude

	  LAT1S =    LAT1*60.D0*60.D0/RHOSEC
	  LONG1S = -LONG1*60.D0*60.D0/RHOSEC
	  LAT2S =    LAT2*60.D0*60.D0/RHOSEC
	  LONG2S = -LONG2*60.D0*60.D0/RHOSEC

	  DLAT  = ( LAT2S -  LAT1S)*RHOSEC
	  DLONG = (LONG2S - LONG1S)*RHOSEC

	  AVLAT = (LAT1S + LAT2S)/2.0D0

	  W  = DSQRT(1.0D0 - E2*DSIN(AVLAT)**2)
	  LM = AXIS*(1.0D0 - E2)/(W**3*RHOSEC)
	  LP = AXIS*DCOS(AVLAT)/(W*RHOSEC)

	  LATMTR = LM*DLAT
	  LONMTR = LP*DLONG

	  RETURN
	  END

**********************************************************************
	  SUBROUTINE NGRIDS (NODATA)
**********************************************************************
* This subroutine opens the NADCON grids which contain datum shifts.
* A total of two files are necessary for each area; 1 for each latitude
* and longitude shift table (gridded data set) expressed in arc seconds.

* 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 DOUBLE PRECISION (A-H, O-Z)
	  IMPLICIT INTEGER (I-N)
*     IMPLICIT UNDEFINED (A-Z)

	  INTEGER MXAREA
	  PARAMETER (MXAREA = 7)

	  LOGICAL NODATA

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

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

* Initialize

	  NODATA = .FALSE.
	  NAREA = 0
	  CALL DGRIDS

	  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, GFLAG, NOGO, DX, DY,
     +                   XMAX1, XMIN1, YMAX1, YMIN1, NC1, CARD)
**********************************************************************
*** Given base name of gridded data files, open the two data files

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

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

	  DOUBLE PRECISION XMAX1, XMIN1, YMAX1, YMIN1, DX, DY
	  REAL DX1, DY1, DX2, DY2
	  REAL X01, Y01, ANGLE1, X02, Y02, ANGLE2
	  INTEGER IFLAG1, IFLAG2, N1, N2, N3, N4
	  INTEGER ITEMP, LRECL, ILA, ILO, IFILE, IOS
	  INTEGER NC1, NR1, NZ1, NC2, NR2, NZ2
	  CHARACTER*80 CARD
	  CHARACTER*69 ALAS, ALOS
	  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)

	  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'
	  ALAS = B69
	  ALAS(1:N2) = AFILE
	  ALAS(N2+1:N2+4) = '.las'
	  ALOS = B69
	  ALOS(1:N2) = AFILE
	  ALOS(N2+1:N2+4) = '.los'

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

* Seconds of latitude grid file

	  LRECL = 256
	  ILA = 2*ITEMP - 1
	  IFILE = ILA + 10
	  LUAREA(ILA) = IFILE
	  INQUIRE (FILE=ALAS, EXIST=EFLAG1, OPENED=OFLAG)
	  IF (.NOT. EFLAG1) GOTO 100
	  IF (OFLAG) GOTO 980
	  OPEN (IFILE,FILE=ALAS,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=ALAS,FORM='UNFORMATTED',STATUS='OLD',
     +       ACCESS='DIRECT',RECL=LRECL,ERR=940,IOSTAT=IOS)

* Seconds of longitude grid file

  100 LRECL = 256
	  ILO = 2*ITEMP
	  IFILE = ILO + 10
	  LUAREA(ILO) = IFILE
	  INQUIRE (FILE=ALOS, EXIST=EFLAG2, OPENED=OFLAG)
	  IF (.NOT. EFLAG1) GOTO 910
	  IF (.NOT. EFLAG2) GOTO 920
	  IF (OFLAG) GOTO 980
	  OPEN (IFILE,FILE=ALOS,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=ALOS,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

* Calculate values used in this program

	  XMIN1 = DBLE(X01)
	  YMIN1 = DBLE(Y01)
	  XMAX1 = DBLE(X01) + (NC1-1)*DBLE(DX1)
	  YMAX1 = DBLE(Y01) + (NR1-1)*DBLE(DY1)
	  DX = DBLE( ABS(DX1) )
	  DY = DBLE( 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 (ALAS, IFLAG1, N1)
	CALL NBLANK (ALAS, IFLAG2, N2)
	CALL NBLANK (CARD, IFLAG1, N3)
	CALL NBLANK (CARD, IFLAG2, N4)
	IF (EFLAG2) THEN

* .los exists, .las does not exist

	  WRITE (LUOUT, 925) ALAS(N1:N2), CARD(N3:N4), ALOS(N1:N2)
	ELSE

* neither .los nor .las exist

	  WRITE (LUOUT, 915) ALAS(N1:N2), ALOS(N1:N2), CARD(N3:N4)
  915     FORMAT (/, 5X, ' *********** WARNING ***********', /,
     +            5X, ' The grid files', /,
     +            6X, '''', A, '''', /,
     +            6X, '''', A, '''', /,
     +            5X, ' from record:', /, 
     +            6X, '''', A, '''', /,
     +            5X, ' do not exist!', /,
     +            5X, ' *******************************', /)
	ENDIF
	  ENDIF
	  CLOSE ( LUAREA(ILA) )
	  CLOSE ( LUAREA(ILO) )
	  GOTO 9999

  920 CONTINUE
	  NOGO = .TRUE.
	  IF (GFLAG) THEN

* .las exists, .los does not exist

	CALL NBLANK (ALAS, IFLAG1, N1)
	CALL NBLANK (ALAS, IFLAG2, N2)
	CALL NBLANK (CARD, IFLAG1, N3)
	CALL NBLANK (CARD, IFLAG2, N4)
	WRITE (LUOUT, 925) ALOS(N1:N2), CARD(N3:N4), ALAS(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(ILA) )
	  CLOSE ( LUAREA(ILO) )
	  GOTO 9999

* Grid file(s) already open

  940 CONTINUE
	  NOGO = .TRUE.
	  IF (GFLAG) THEN
	CALL NBLANK (ALAS, IFLAG1, N1)
	CALL NBLANK (ALAS, IFLAG2, N2)
	CALL NBLANK (CARD, IFLAG1, N3)
	CALL NBLANK (CARD, IFLAG2, N4)
	WRITE (LUOUT, 945) ALAS(N1:N2), ALOS(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(ILA) )
	  CLOSE ( LUAREA(ILO) )
	  GOTO 9999

* Grid files do not agree

  960 CONTINUE
	  NOGO = .TRUE.
	  CALL NBLANK (ALAS, IFLAG1, N1)
	  CALL NBLANK (ALAS, IFLAG2, N2)
	  WRITE (LUOUT, 965) ALAS(N1:N2), ALOS(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(ILA) )
	  CLOSE ( LUAREA(ILO) )
	  GOTO 9999

* Grid files already open

  980 CONTINUE
	  NOGO = .TRUE.
	  IF (GFLAG) THEN
	CALL NBLANK (ALAS, IFLAG1, N1)
	CALL NBLANK (ALAS, IFLAG2, N2)
	WRITE (LUOUT, 985) ALAS(N1:N2), ALOS(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 TO83 (NOGO, XPT, YPT, XPT2, YPT2,
     +                 DLAM, DLOM, DLAS, DLOS)
**********************************************************************
* This subroutine predicts the NAD 83 latitude and longitude values
* given the NAD 27 latitude and longitude values in degree decimal
* format.  In addition, the program returns the shift values between
* The datums in both arc secs and meters.

* All of the predictions are based upon a straight-forward interpolation
* of a gridded data set of datum shifts.  The datum shifts are assumed
* to be provided in the files opened in the NGRIDS subroutine.  The
* common AREAS contains the names of the valid areas while the common
* GDINFO contains the grid variables.  NAREA is the number of areas
* which had data files opened.  A total of two files are necessary for
* each area: one latitude and one longitude shift table (gridded data
* set) expressed in arc seconds.

* For this subroutine, it is important to remember that the
* input longitude is assumed to be positive east and the
* output longitude will be positive east.

*       Author:     Warren T. Dewhurst, PH. D.
*                   National Geodetic Survey
*                   November 1, 1989

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

	  INTEGER MXAREA
	  PARAMETER (MXAREA = 7)

	  DOUBLE PRECISION XPT, YPT, XPT2, YPT2
	  DOUBLE PRECISION XGRID, YGRID
	  DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS
	  DOUBLE PRECISION DX0, DY0, XMAX0, XMIN0, YMAX0, YMIN0
	  INTEGER IROW, JCOL, IAREA, I, NC0, ITYPE
	  INTEGER IFLAG1, IFLAG2, N1, N2
	  LOGICAL NOGO, FLAG

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

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

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

	  CHARACTER*80 CARD
	  COMMON /CURNT/ CARD

	  SAVE FLAG

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

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

	  NOGO  =  .FALSE.
*      ITYPE = 0

****************************************************
* 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)
	NC0 = NC(IAREA)

	CALL FGRID (XPT, YPT, DX0, DY0, XMAX0, XMIN0,
     +              YMAX0, YMIN0, XGRID, YGRID, IROW, JCOL, NOGO)
	IF (.NOT. NOGO) GOTO 200
	WRITE (LUOUT,28) YPT,XPT,IAREA
   28   FORMAT ('THE POINT ',F10.5,2X,F10.5,' IS NOT IN AREA ',I2)
  100 CONTINUE

* Not in any of the grid areas

	  NOGO = .TRUE.
	  GOTO 950

  200 CONTINUE

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

	CALL INTRP (IAREA, IROW, NC0, JCOL, XGRID, YGRID,
     +              XPT, YPT, XPT2, YPT2, DLOS, DLAS, DLAM, DLOM)

 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
	  GOTO 9999
	  END

**********************************************************************
	  SUBROUTINE TRANSF (NOGO, XPT, YPT, XPT2, YPT2,
     +                   DLAM, DLOM, DLAS, DLOS, KEY)
**********************************************************************
* This subroutine computes either the forward or inverse coordinate
* transformation depending upon the value of the integer variable 'key'

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

	  INTEGER MXAREA, ITMAX
	  DOUBLE PRECISION SMALL
	  PARAMETER (MXAREA = 7, ITMAX = 10, SMALL = 1.0D-9 )

	  DOUBLE PRECISION XPT, YPT, XPT2, YPT2
	  DOUBLE PRECISION XTEMP, YTEMP, XDIF, YDIF
	  DOUBLE PRECISION DLAM, DLOM, DLAS, DLOS
	  DOUBLE PRECISION DXLAST, DYLAST
	  INTEGER*2 KEY, NUM
	  LOGICAL NOGO

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

	  IF (KEY .EQ. 1) THEN

**********************
* FOR NAD 27 TO NAD 83
**********************

	CALL TO83 (NOGO, XPT, YPT, XPT2, YPT2,
     +             DLAM, DLOM, DLAS, DLOS)

	  ELSEIF (KEY .EQ. 2) THEN

***************************
* FOR NAD 83 TO NAD 27)
* THIS IS DONE BY ITERATION
***************************

	NUM = 0

**************************************************
* SET THE XPT,YPT TO TEMPORARY VALUES
* (REMEMBER, XPT AND YPT ARE REALLY NAD 83 VALUES)
**************************************************

	XTEMP = XPT
	YTEMP = YPT

**************************************************************
* PRETEND THAT THESE TEMPORARY VALUES ARE REALLY NAD 27 VALUES
* FOR A FIRST GUESS AND COMPUTE PSEUDO-NAD 83 COORDINATES
**************************************************************

  200   CONTINUE
	  NUM = NUM + 1

********************************
* CHECK THE NUMBER OF ITERATIONS
********************************

	  IF (NUM .GE. ITMAX) THEN
		WRITE (LUOUT,*) ' *** MAXIMUM ITERATIONS EXCEEDED!! ***'
		WRITE (LUOUT,*) ' *** CALL PROGRAMMER FOR HELP ***'
		WRITE (LUOUT,*) ' LATITUDE =', YTEMP, ' LONGITUDE =', XTEMP
		NOGO = .TRUE.
		GOTO 1000
	  ENDIF

	  CALL TO83 (NOGO, XTEMP, YTEMP, XPT2, YPT2,
     +               DLAM, DLOM, DLAS, DLOS)
	  DXLAST = DLOS
	  DYLAST = DLAS

**************************************
* COMPARE TO ACTUAL NAD 83 COORDINATES
**************************************

	  XDIF = XPT - XPT2
	  YDIF = YPT - YPT2

****************************************************************
* COMPUTE A NEW GUESS UNLESS THE DIFFERENCES ARE LESS THAN SMALL
* WHERE SMALL IS DEFINED (ABOVE) TO BE;  SMALL = 1.0D-9
****************************************************************

	  IF (NUM .EQ. 1) THEN
		IF (DABS(XDIF) .GT. SMALL) THEN
		  XTEMP = XPT - DLOS/3600.D0
		ENDIF
		IF (DABS(YDIF) .GT. SMALL) THEN
		  YTEMP = YPT - DLAS/3600.D0
		ENDIF
	  ELSE
		IF (DABS(XDIF) .GT. SMALL) THEN
		  XTEMP = XTEMP - (XPT2 - XPT)
		ENDIF
		IF (DABS(YDIF) .GT. SMALL) THEN
		  YTEMP = YTEMP - (YPT2 - YPT)
		ENDIF

	  ENDIF

	  IF (DABS(YDIF) .LE. SMALL  .AND.  DABS(XDIF) .LE. SMALL) THEN

******************************
* IF CONVERGED THEN LEAVE LOOP
******************************

		XPT = XTEMP
		YPT = YTEMP
		GOTO 1000
	  ENDIF

*******************************
* IF NOT CONVERGED THEN ITERATE
*******************************

	GOTO 200

	  ENDIF
 1000 RETURN
	  END

**********************************************************************
	SUBROUTINE SURF (XGRID, YGRID, ZEE, AY, BEE, CEE, DEE, IROW, JCOL)
**********************************************************************
** SUBROUTINE SURF: INTERPOLATES THE Z VALUE                         *
**********************************************************************

* Calculated the value of the grid at the point XPT, YPT.  The
* interpolation is done in the index coordinate system for convenience.

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

	  DOUBLE PRECISION XGRID, YGRID
	  DOUBLE PRECISION AY, BEE, CEE, DEE
	  DOUBLE PRECISION ZEE, ZEE1, ZEE2, ZEE3, ZEE4
	  INTEGER IROW, JCOL

	  ZEE1 = AY
	  ZEE2 = BEE*(XGRID - DBLE(JCOL) )
	  ZEE3 = CEE*(YGRID - DBLE(IROW) )
	  ZEE4 = DEE*(XGRID - DBLE(JCOL) )*(YGRID - DBLE(IROW) )
	  ZEE  = ZEE1 + ZEE2 + ZEE3 + ZEE4

	  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 DOUBLE PRECISION (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
