	 PROGRAM UTMS
*
*     THIS PROGRAM CONVERTS GPS TO UNIVERSIAL TRANSVERSE MERACTOR COORDINATES
*     AND VICE VERSA FOR THE NAD27 AND NAD83 DATUM.
*     THIS PROGRAM WAS WRITTEN BY E. CARLSON
*     SUBROUTINES TMGRID, TCONST, TMGEOD, TCONPC,
*     WERE WRITTEN BY T. VINCENTY, NGS, IN JULY 1984
*     AND LAST UPDATED IN FEBUARY 1986.
*
*****************************************************************     *
*                  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 DOUBLE PRECISION(A-H,O-Z)
	CHARACTER*60 INFILE
	INTEGER*2 UNITS,DATNUM,INUM,ICODE
	COMMON/CONST/RAD,ER,RF,ESQ,PI
	COMMON/DATUM/DATNUM,UNITS
 
	PI=4.D0*DATAN(1.D0)
	RAD=180.D0/PI
c
c  1 geo to UTM, 
c  2 UTM to geo,
c
	WRITE(*,*) 'STARTING CORPSUTM'
	READ(5,15) INUM
   15 FORMAT(I1)
c 
c  NAD27 EQUALS NUMBER ONE (1) 
c  NAD83 EQUALS NUMBER TWO (2) 
c
	READ(5,15) DATNUM
c
c  UNITS OF MEASURE
c
	READ(5,15) UNITS
c
c  INPUT FILE
c
	READ(5,16) INFILE
   16 FORMAT(A60)
c
c  ZONE
c
	READ(5,17) ICODE
   17 FORMAT(I5)
c
c  OPEN FILES
c
	WRITE(*,*) 'CORPSUTM INPUT PARAMETERS READ SUCCESSFULLY'
	OPEN(91,FILE=INFILE,ACCESS='SEQUENTIAL',FORM='BINARY',
     &ERR=25,STATUS='OLD')
	WRITE(*,*) 'CORPSUTM INPUT FILE OPENED SUCCESSFULLY'
	OPEN(93,FILE='INT$OUT',ACCESS='SEQUENTIAL',FORM='BINARY',
     &ERR=25,STATUS='UNKNOWN')
	WRITE(*,*) 'CORPSUTM OUTPUT FILE OPENED SUCCESSFULLY'
**  
**     FIND THE RIGHT SIMI MAJOR AXIS AND FLATTING 
**
	 IF(DATNUM.EQ.1) THEN
**      FOR THE NAD 27 DATUM
**
		ER=6378206.4D0
		RF=294.978698D0
		F=1.D0/RF
		ESQ=(F+F-F*F)
 
	ELSEIF (DATNUM.EQ.2) THEN
**     FOR THE NAD83 DATUM 
**   
		ER=6378137.D0
		RF=298.257222101D0
		F=1.D0/RF
		ESQ=(F+F-F*F)
	 ENDIF
***
***   USE THE NUM TO DO THE CORRECT FUNCTION    

	IF(INUM.EQ.1) THEN
	 CALL BUTGP(ICODE)
	ELSEIF(INUM.EQ.2) THEN
	 CALL BUTPC(ICODE)
	ENDIF
 25   CLOSE(91)
	CLOSE(93)
	WRITE(*,*) 'CORPSUTM INPUT/OUTPUT FILES CLOSED'
	WRITE(*,*) 'CORPSUTM COMPLETE'
	STOP
	END
******************************************************************
	SUBROUTINE BUTPC(ICODE)
***************************************************************
	IMPLICIT REAL*8 (A-H,O-Z)
	REAL*8 NORTH,EAST
	REAL*8 CORPSD(29)
	INTEGER*2 FLAG,ICODE
	INTEGER*2 DATNUM,UNITS
	CHARACTER*25 NAME
	COMMON/XY/NORTH,EAST
	COMMON/CORPS/NAME,CORPSD,FLAG
	COMMON/DATUM/DATNUM,UNITS
	 
   40 READ(91,END=99) NAME,(CORPSD(II),II=1,29),FLAG
	IF(DATNUM.EQ.1) THEN
	  NORTH = CORPSD(14)
	  EAST = CORPSD(13)
	ELSE
	  NORTH = CORPSD(18)
	  EAST = CORPSD(17)
	ENDIF
	IF(UNITS.EQ.1) THEN
	  NORTH = NORTH * (1200.0D0/3937.0D0)
	  EAST =  EAST * (1200.0D0/3937.0D0)
	ENDIF
	IF(UNITS.EQ.2) THEN
	  NORTH =  NORTH * (381.0D0/1250.0D0)
	  EAST =  EAST * (381.0D0/1250.0D0)
	ENDIF
	CALL DRUTGP(ICODE)
	GO TO 40
   99 RETURN
	END
***********************************************************
	SUBROUTINE BUTGP(ICODE)
**********************************************************
	IMPLICIT REAL*8 (A-H,O-Z)
	REAL*8 CORPSD(29),DDLAT,DDLON
	INTEGER*2 FLAG
	INTEGER*2 DATNUM,UNITS,II,ICODE
	CHARACTER*25 NAME
	COMMON/LATLON/DDLAT,DDLON
	COMMON/CORPS/NAME,CORPSD,FLAG
	COMMON/DATUM/DATNUM,UNITS 

   40 READ(91,END=99) NAME,(CORPSD(II),II=1,29),FLAG
	IF(DATNUM.EQ.1) THEN
	     DDLAT = CORPSD(9)
	     DDLON = CORPSD(10)
	ELSE
	     DDLAT = CORPSD(11)
	     DDLON = CORPSD(12)
	ENDIF
	CALL DRGPUT(ICODE)
	GO TO 40
   99 RETURN
	END
************************************************************************
	 SUBROUTINE DRGPUT(ICODE)
**********************************************************************
*
*
*      THIS IS THE DRIVER TO COMPUTE UTM NORTHINGS AND EASTINGS
*      FOR EACH PRIMARY ZONE AND THE AJACENT ZONE IF THE LONGITUDE
*      IS WITH 5 MINUTES OF THE ZONE BOUNDARIES
*
*      THE OUTPUT IS FOR THE DATA SHEET PROGRAM
*
*      VARIABLES
*      CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS
*      ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS)
*      RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD
*      ESQ= E SQUARED
*      RAD = RADIAN CONVERSION FACTOR
*      CM = CENTRAL MERIDIAN ( COMPUTED USEING THE LONGITUDE)
*      SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM)
*      OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM)
*      R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING
*                            MERIDIONAL DISTANCE FROM LATITUDE
*      SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR )
*           FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE
*           ( ALWAYS ZERO FOR UTM)
*

	 IMPLICIT REAL*8 (A-H,O-Z)

	 CHARACTER*4 ZONE
	 INTEGER*4 FOUND
	 REAL*8 DDLAT,DDLON,FI,LAM,LCM,UCM
	 REAL*8 KP
	 REAL*8 GRIDAZ
	 REAL*8 NORTH,EAST
	 INTEGER*2 ICODE,FLAG

	 COMMON/CONST/RAD,ER,RF,ESQ,PI
	 COMMON/LATLON/DDLAT,DDLON
*
*      CONVERT THE LATITUDE AND LONGITUDE TO PI AND LAM
*
	 FI=(DDLAT)/RAD
	 LAM=(DDLON)/RAD

*
*      FIND THE ZONE FOR LONGITUDE LESS THAN 180 DEGREES
*
	IF(DDLON.LT.180) THEN
	   IZ=DDLON/6
	   IZ= 30 -IZ
	   IF(ICODE.NE.IZ) THEN
		FLAG = 8
	   ELSE
		FLAG = 0
	   ENDIF
	   IZ = ICODE
	   ICM=(183-(6*IZ))
	   CM=DBLE(FLOAT(ICM))/RAD
	   UCM=(ICM+3)/RAD
	   LCM=(ICM-3)/RAD
	ENDIF
*
*       FIND THE ZONE FOR LONGITUDE GREATER THAN 180 DEGREES
*
*        IF(DDLON.GE.180) THEN
*           IZ=(LOD)/6
*           IZ= 90 - IZ
*           ICM=(543 - (6*IZ))
*           CM= DBLE(FLOAT(ICM))/RAD
*           UCM=(ICM+3)/RAD
*           LCM=(ICM-3)/RAD
*        ENDIF

	 TOL=(5.0D0/60.0D0)/RAD

	 FN=0.D0
	 FE=500000.0D0
	 SF=0.9996D0
	 OR=0.0D0

	 FOUND=0

	 CALL TCONST(ER,RF,SF,OR,ESQ,EPS,R,A,B,C,U,V,W,SO,
     &            CM,FE,FN)

*
*      COMPUTE THE NORTH AND EASTINGS

*

 200   CALL TMGRID(FI,LAM,NORTH,EAST,CONV,KP,ER,ESQ,EPS,CM,
     &            FE,FN,SF,SO,R,A,B,C,U,V,W)


*
*      WRITE THE ZONE NUMBER
*
	 IF (IZ.GT.9) THEN
	 WRITE(ZONE,600) IZ
 600     FORMAT(1X,I2)
	 ELSE
	 WRITE(ZONE,605) IZ
 605     FORMAT(1X,I2.2)
	 ENDIF

*
*      WRITE THE OUTPUT TO THE PLANE FILE FOR THE DATA SHEET
*      PROGRAM
*
	 CALL DATAUT(NORTH,EAST,CONV,KP,IZ,FLAG)

*
*      DO THE TEST TO SEE IF THE LONGITUDE IS WITHIN 5 MINUTES
*      OF THE BOUNDARIES FOR THE ZONE AND IF SO COMPUTE THE
*      NORHT AND EASTING FOR THE ADJACENT ZONE
*
C       IF(FOUND.NE.0) THEN
C          RETURN
C       ENDIF
C
C       IF(DABS(UCM-LAM).LE.TOL) THEN
C          CM=DBLE(FLOAT(ICM+6))/RAD
C          IZ=IZ-1
C          IF(IZ.EQ.0) IZ=60
C          FOUND=FOUND+1
C          GO TO 200
C       ENDIF
C
C       IF(DABS(LCM-LAM).LE.TOL) THEN
C          CM=DBLE(FLOAT(ICM-6))/RAD
C          IZ=IZ+1
C          IF(IZ.EQ.61) IZ=1
C          FOUND=FOUND+1
C          GO TO 200
C       ENDIF
C

	 RETURN
	 END

*******************************************************************
	 SUBROUTINE DRUTGP(ICODE)
*
*
*
*      THIS IS THE DRIVER TO COMPUTE LATITUDES AND LONGITUDES FROM
*      THE UTMs FOR EACH ZONE
*
*
*      VARIABLES
*      CARDR = A MODIFIED 80 RECORD CARD WITH A LENGTH OF 211 COLS
*      ER = EQUATORIAL RADIUS OF THE ELLIPSOID (SEMI-MAJOR AXIS)
*      RF = RECIPROCAL OF FLATTING OF THE ELLIPSOD
*      ESQ= E SQUARED
*      RAD = RADIAN CONVERSION FACTOR
*      CM = CENTRAL MERIDIAN ( COMPUTED USEING THE LONGITUDE)
*      SF = SCALE FACTOR OF CENTRAL MERIDIAN ( ALWAYS .9996 FOR UTM)
*      OR = SOUTHERNMOST PARALLEL OF LATITUDE ( ALWAYS ZERO FOR UTM)
*      R, A, B, C, U, V, W = ELLIPSOID CONSTANTS USED FOR COMPUTING
*                            MERIDIONAL DISTANCE FROM LATITUDE
*      SO = MERIDIONAL DISTANCE (MULTIPLIED BY SCALE FACTOR )
*           FROM THE EQUATOR TO THE SOUTHERNMOST PARALLEL OF LATITUDE
*           ( ALWAYS ZERO FOR UTM)
*

	 IMPLICIT REAL*8 (A-H,O-Z)
	 CHARACTER*4 ZONE
	 INTEGER*4 FOUND
	 INTEGER*2 ICODE

	 REAL*8 LAT,LON,LCM,UCM
	 REAL*8 KP
	 REAL*8 NORTH,EAST
 
	 COMMON/CONST/RAD,ER,RF,ESQ,PI
	 COMMON/XY/NORTH,EAST
 

*
*      FIND THE CENTRAL MERIDAIN IF THE ZONE NUMBER IS LESS THAN 30
*
	IF(ICODE.LT.30) THEN
	   IZ=ICODE
	   ICM=(183-(6*IZ))
	   CM=DBLE(FLOAT(ICM))/RAD
	   UCM=(ICM+3)/RAD
	   LCM=(ICM-3)/RAD
	ENDIF
*
*       FIND THE CENTRAL MERIDAN IF THE ZONE NUMBER IS LARGER THAN 30
*
	IF(ICODE.GE.30) THEN
	   IZ=ICODE
	   ICM=(543 - (6*IZ))
	   CM= DBLE(FLOAT(ICM))/RAD
	   UCM=(ICM+3)/RAD
	   LCM=(ICM-3)/RAD
	ENDIF

	 TOL=(5.0D0/60.0D0)/RAD

	 FN=0.D0
	 FE=500000.0D0
	 SF=0.9996D0
	 OR=0.0D0

	 FOUND=0

	 CALL TCONPC (SF,OR,EPS,R,SO,V0,V2,V4,V6,ER,ESQ,RF)

*
*      COMPUTE THE LATITUDES AND LONGITUDES
*

 200   CALL TMGEOD (NORTH,EAST,LAT,LON,EPS,CM,FE,SF,SO,R,V0,V2,
     &              V4,V6,FN,ER,ESQ,CONV,KP)

*
*      WRITE THE ZONE NUMBER
*
	 IF (IZ.GT.9) THEN
	 WRITE(ZONE,600) IZ
 600     FORMAT(1X,I2)
	 ELSE
	 WRITE(ZONE,605) IZ
 605     FORMAT(1X,I2.2)
	 ENDIF

*
*      WRITE THE OUTPUT TO THE PLANE FILE FOR THE DATA SHEET
*      PROGRAM
*
	 CALL DATAGP(LAT,LON,IZ,CONV,KP)

*
*      DO THE TEST TO SEE IF THE LONGITUDE IS WITHIN 5 MINUTES
*      OF THE BOUNDARIES FOR THE ZONE AND IF SO COMPUTE THE
*      NORHT AND EASTING FOR THE ADJACENT ZONE
*
C       IF(FOUND.NE.0) THEN
C          RETURN
C       ENDIF
C
C       IF(DABS(UCM-LAM).LE.TOL) THEN
C          CM=DBLE(FLOAT(ICM+6))/RAD
C          IZ=IZ-1
C          IF(IZ.EQ.0) IZ=60
C          FOUND=FOUND+1
C          GO TO 200
C       ENDIF
C
C       IF(DABS(LCM-LAM).LE.TOL) THEN
C          CM=DBLE(FLOAT(ICM-6))/RAD
C          IZ=IZ+1
C          IF(IZ.EQ.61) IZ=1
C          FOUND=FOUND+1
C          GO TO 200
C       ENDIF
C
	 RETURN
	 END

**********************************************************************
	SUBROUTINE DATAUT(NORTH,EAST,CONV,KP,ZONE,FLAG2)
*********************************************************************
*
	IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER*25 NAME
	REAL*8 CORPSD(29)
	INTEGER*2 FLAG,FLAG2
	INTEGER*4 IDEG,IMIN
	REAL*8 NORTH,EAST,CONV,KP,CSEC
	INTEGER*4 ZONE
	INTEGER*2 DATNUM,UNITS
	COMMON/DATUM/DATNUM,UNITS
	COMMON/CORPS/NAME,CORPSD,FLAG

	CALL TODMS(CONV,IDEG,IMIN,CSEC)
	IF(UNITS.EQ.1) THEN
	  NORTH = NORTH * (3937.0D0/1200.0D0)
	  EAST = EAST * (3937.0D0/1200.0D0)
	ENDIF
	IF(UNITS.EQ.2) THEN
	  NORTH = NORTH * (1250.0D0/381.0D0)
	  EAST = EAST * (1250.0D0/381.0D0)
	ENDIF
	FLAG = IOR(FLAG,FLAG2)
	IF(DATNUM.EQ.1) THEN
	  CORPSD(13) = EAST
	  CORPSD(14) = NORTH
	  IF(CONV.LT.0) THEN
		CORPSD(15) = DBLE(IDEG - IMIN/60.0D0 - CSEC/3600.0D0)
	  ELSE
		CORPSD(15) = DBLE(IDEG + IMIN/60.0D0 + CSEC/3600.0D0)
	  ENDIF
	  CORPSD(16) = KP
	ELSE
	  CORPSD(17) = EAST
	  CORPSD(18) = NORTH
	  IF(CONV.LT.0) THEN
		CORPSD(19) = DBLE(IDEG - IMIN/60.0D0 - CSEC/3600.0D0)
	  ELSE
		CORPSD(19) = DBLE(IDEG + IMIN/60.0D0 + CSEC/3600.0D0)
	  ENDIF
	  CORPSD(20) = KP
	ENDIF
	WRITE (93) NAME,CORPSD,FLAG
	RETURN
	END
*********************************************************************
	SUBROUTINE DATAGP(LAT,LON,ZONE,CONV,KP)
**********************************************************************
***
	IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER*25 NAME
	REAL*8 CORPSD(29)
	INTEGER*2 FLAG
	REAL*8 LAT,LON,CONV,KP,CSEC,SLAT,SLON
	INTEGER*4 ZONE
	INTEGER*2 DATNUM,UNITS
	INTEGER*4 IDEG,IMIN,LD,LM,LOD,LOM
	COMMON/DATUM/DATNUM,UNITS
	COMMON/CORPS/NAME,CORPSD,FLAG

	CALL TODMS(LAT,LD,LM,SLAT)
	CALL TODMS(LON,LOD,LOM,SLON)
	CALL TODMS(CONV,IDEG,IMIN,CSEC)
	IF(DATNUM.EQ.1) THEN
		 CORPSD(9) = DBLE(LD+LM/60.0D0+SLAT/3600.0D0)
		 CORPSD(10) = DBLE(LOD+LOM/60.0D0+SLON/3600.0D0)
		 IF(CONV.LT.0) THEN
			CORPSD(15) = DBLE(IDEG - IMIN/60.0D0 - CSEC/3600.0D0)
		 ELSE
			CORPSD(15) = DBLE(IDEG + IMIN/60.0D0 + CSEC/3600.0D0)
		 ENDIF
		 CORPSD(16) = KP
	 ELSE
		 CORPSD(11) = DBLE(LD+LM/60.0D0+SLAT/3600.0D0)
		 CORPSD(12) = DBLE(LOD+LOM/60.0D0+SLON/3600.0D0)
		 IF(CONV.LT.0) THEN
			CORPSD(19) = DBLE(IDEG - IMIN/60.0D0 - CSEC/3600.0D0)
		 ELSE
			CORPSD(19) = DBLE(IDEG + IMIN/60.0D0 + CSEC/3600.0D0)
		 ENDIF
		 CORPSD(20) = KP
	ENDIF
	 WRITE (93) NAME,CORPSD,FLAG
	RETURN
	END
*********************************************************************
	SUBROUTINE TODMS(RAD,IDG,MIN,SEC)
*********************************************************************
C     RADIANS TO DEGREES,MINUTES AND SECONDS
C
	REAL*8 RAD,SEC,RHOSEC
	DATA RHOSEC/2.062648062471D05/
	SEC=RAD*RHOSEC
	IDG=SEC/3600.D0
	SEC=SEC-DBLE(IDG*3600)
	MIN=SEC/60.D0
	SEC=SEC-DBLE(MIN*60)
	IF((60.D0-DABS(SEC)).GT.5.D-6) GO TO 100
	SEC=SEC-DSIGN(60.D0,SEC)
	MIN=MIN+ISIGN(1,MIN)
  100 IF(IABS(MIN).LT.60) GO TO 101
	MIN=MIN-ISIGN(60,MIN)
	IDG=IDG+ISIGN(1,IDG)
  101 MIN=IABS(MIN)
	SEC=DABS(SEC)
	IF(RAD.GE.0.D0) GO TO 102
	IF(IDG.EQ.0) MIN=-MIN
	IF(IDG.EQ.0.AND.MIN.EQ.0)SEC=-SEC
  102 RETURN
	END
********************************************************************
	SUBROUTINE VERDMS(VAL,ID,IM,S,ISIGN)
********************************************************************
***
*** CONVERT RADIANS TO DEG, MIN, SEC
 
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	COMMON/CONST/RAD,ER,RF,ESQ,PI
 
    1 IF(VAL.GT.PI) THEN
	VAL=VAL-PI-PI
	GO TO 1
	ENDIF
 
    2 IF(VAL.LT.-PI) THEN
	VAL=VAL+PI+PI
	GO TO 2
	ENDIF
 
	IF(VAL.LT.0.D0) THEN
	ISIGN=-1
	ELSE
	ISIGN=+1
	ENDIF
 
	S=DABS(VAL*RAD)
	ID=IDINT(S)
	S=(S-ID)*60.D0
	IM=IDINT(S)
	S=(S-IM)*60.D0
 
	RETURN
	END
	SUBROUTINE TMGRID(FI,LAM,NORTH,EAST,CONV,KP,ER,ESQ,EPS,CM,
     &                  FE,FN,SF,SO,R,A,B,C,U,V,W)
 
	IMPLICIT DOUBLE PRECISION(A-H,K-Z)
C
C*****  TRANSVERSE MERCATOR PROJECTION
C       CONVERSION OF GEODETIC COORDINATES TO GRID COORDINATES
C*****  Programmed by T. Vincenty, NGS, in July 1984.
C*****************  SYMBOLS AND DEFINITIONS *************************
C   Latitude positive north, longitude positive west.  All angles are
C     in radian measure.
C   N, E are northing and easting coordinates respectively.
C   LAT, LON are latitude and longitude respectively.
C   CONV is convergence.
C   KP is point scale factor.
C   ER is equatorial radius of the ellipsoid (= major semiaxis).
C   ESQ is the square of first eccentricity of the ellipsoid.
C   EPS is the square of second eccentricity of the ellipsoid.
C   CM is the central meridian of the projection zone.
C   FE is false easting value at the central meridian.
C   FN is "false northing" at the southernmost latitude, usually zero.
C   SF is scale factor at the central meridian.
C   SO is meridional distance (multiplied by the scale factor) from
C     the equator to the southernmost parallel of latitude for the zone.
C   R is the radius of the rectifying sphere (used for computing
C     meridional distance from latitude and vice versa).
C   A, B, C, U, V, W are other precomputed constants for determination
C     of meridional distance from latitude and vice versa.
C
C   The formula used in this subroutine gives geodetic accuracy within
C   zones of 7 degrees in east-west extent.  Within State transverse
C   Mercator projection zones, several minor terms of the equations
C   may be omitted (see a separate NGS publication).  If programmed
C   in full, the subroutine can be used for computations in surveys
C   extending over two zones.
C
C*********************************************************************
	OM=FI + A*SIN(2.*FI) + B*SIN(4.*FI) + C*SIN(6.*FI)
	S=R*OM*SF
	SINFI=SIN(FI)
	COSFI=COS(FI)
	TN=SINFI/COSFI
	TS=TN**2
	ETS=EPS*COSFI**2
	L=(LAM-CM)*COSFI
	LS=L*L
	RN=SF*ER/SQRT(1.-ESQ*SINFI**2)
C
	A2=RN*TN/2.
	A4=(5.-TS+ETS*(9.+4.*ETS))/12.
	A6=(61.+TS*(TS-58.)+ETS*(270.-330.*TS))/360.
	A1=-RN
	A3=(1.-TS+ETS)/6.
	A5=(5.+TS*(TS-18.)+ETS*(14.-58.*TS))/120.
	A7=(61.-479.*TS+179.*TS**2-TS**3)/5040.
	NORTH=S-SO + A2*LS*(1.+LS*(A4+A6*LS)) +FN
	EAST=FE + A1*L*(1.+ LS*(A3+LS*(A5+A7*LS)))
C
C*** CONVERGENCE
	C1=-TN
	C3=(1.+3.*ETS+2.*ETS**2)/3.
	C5=(2.-TS)/15.
	CONV=C1*L*(1.+LS*(C3+C5*LS))
C
C*** POINT SCALE FACTOR
	F2=(1.+ETS)/2.
	F4=(5.-4.*TS+ETS*( 9.-24.*TS))/12.
	KP=SF*(1.+F2*LS*(1.+F4*LS))
C
	RETURN
	END
	SUBROUTINE TCONPC(SF,OR,EPS,R,SO,V0,V2,V4,V6,ER,ESQ,RF)
 
***          TRANSVERSE MERCATOR PROJECTION               ***
*** CONVERSION OF GRID COORDS TO GEODETIC COORDS
*** REVISED SUBROUTINE OF T. VINCENTY  FEB. 25, 1985
************** SYMBOLS AND DEFINITIONS ***********************
*** ER IS THE SEMI-MAJOR AXIS FOR GRS-80
*** SF IS THE SCALE FACTOR AT THE CM
*** SO IS THE MERIDIANAL DISTANCE (TIMES THE SF) FROM THE
***       EQUATOR TO SOUTHERNMOST PARALLEL OF LAT. FOR THE ZONE
*** R IS THE RADIUS OF THE RECTIFYING SPHERE
*** U0,U2,U4,U6,V0,V2,V4,V6 ARE PRECOMPUTED CONSTANTS FOR
***   DETERMINATION OF MERIDIANAL DIST. FROM LATITUDE
*** OR IS THE SOUTHERNMOST PARALLEL OF LATITUDE FOR WHICH THE
***       NORTHING COORD IS ZERO AT THE CM
**************************************************************
 
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 
	F=1.D0/RF
	EPS=ESQ/(1.D0-ESQ)
	PR=(1.D0-F)*ER
	EN=(ER-PR)/(ER+PR)
	EN2=EN*EN
	EN3=EN*EN*EN
	EN4=EN2*EN2
 
	C2=-3.D0*EN/2.D0+9.D0*EN3/16.D0
	C4=15.D0*EN2/16.D0-15.D0*EN4/32.D0
	C6=-35.D0*EN3/48.D0
	C8=315.D0*EN4/512.D0
	U0=2.D0*(C2-2.D0*C4+3.D0*C6-4.D0*C8)
	U2=8.D0*(C4-4.D0*C6+10.D0*C8)
	U4=32.D0*(C6-6.D0*C8)
	U6=128.D0*C8
 
	C2=3.D0*EN/2.D0-27.D0*EN3/32.D0
	C4=21.D0*EN2/16.D0-55.D0*EN4/32.D0
	C6=151.D0*EN3/96.D0
	C8=1097.D0*EN4/512.D0
	V0=2.D0*(C2-2.D0*C4+3.D0*C6-4.D0*C8)
	V2=8.D0*(C4-4.D0*C6+10.D0*C8)
	V4=32.D0*(C6-6.D0*C8)
	V6=128.D0*C8
 
	R=ER*(1.D0-EN)*(1.D0-EN*EN)*(1.D0+2.25D0*EN*EN+
     &     (225.D0/64.D0)*EN4)
	COSOR=DCOS(OR)
	OMO=OR+DSIN(OR)*COSOR*(U0+U2*COSOR*COSOR+U4*COSOR**4+
     &    U6*COSOR**6)
	SO=SF*R*OMO
 
	RETURN
	END
*********************************************************************
	SUBROUTINE TCONST (ER,RF,SF,OR,ESQ,EPS,R,A,B,C,U,V,W,SO,
     &                   CM,FE,FN)
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C***** TRANSVERSE MERCATOR PROJECTION
C      PRECOMPUTATION OF CONSTANTS
C***** Programmed by T. Vincenty, NGS, in July 1984.
C******************** SYMBOLS AND DEFINITIONS  **********************
C   ER is equatorial radius of the ellipsoid (= major semiaxis).
C   RF is reciprocal of flattening of the ellipsoid.
C   SF is scale factor of the central meridian.
C   OR is southernmost parallel of latitude (in radians) for which
C     the northing coordinate is zero at the central meridian.
C   R, A, B, C, U, V, W are ellipsoid constants used for computing
C     meridional distance from latitude and vice versa.
C   SO is meridional distance (multiplied by the scale factor) from
C     the equator to the southernmost parallel of latitude.
C******************************************************************
C
	F=1./RF
	ESQ=(F+F-F**2)
	EPS=ESQ/(1.-ESQ)
	PR=(1.-F)*ER
	EN=(ER-PR)/(ER+PR)
	A=-1.5D0*EN + (9./16.)*EN**3
	B= 0.9375D0*EN**2 - (15./32.)*EN**4
	C=-(35./48.)*EN**3
	U=1.5D0*EN - (27./32.)*EN**3
	V=1.3125D0*EN**2 - (55./32.)*EN**4
	W=(151./96.)*EN**3
	R=ER*(1.-EN)*(1.-EN**2)*(1.+2.25D0*EN**2+(225./64.)*EN**4)
	OMO=OR + A*SIN(2.*OR) + B*SIN(4.*OR) + C*SIN(6.*OR)
	SO=SF*R*OMO
C
	RETURN
	END
	SUBROUTINE TMGEOD(N,E,LAT,LON,EPS,CM,FE,SF,SO,R,V0,V2,
     &                  V4,V6,FN,ER,ESQ,CONV,KP)
 
***          TRANSVERSE MERCATOR PROJECTION               ***
*** CONVERSION OF GRID COORDS TO GEODETIC COORDS
*** REVISED SUBROUTINE OF T. VINCENTY  FEB. 25, 1985
************** SYMBOLS AND DEFINITIONS ***********************
*** LATITUDE POSITIVE NORTH, LONGITUDE POSITIVE WEST.  ALL
***          ANGLES ARE IN RADIAN MEASURE.
*** LAT,LON ARE LAT. AND LONG. RESPECTIVELY
*** N,E ARE NORTHING AND EASTING COORDINATES RESPECTIVELY
*** K IS POINT SCALE FACTOR
*** ER IS THE SEMI-MAJOR AXIS FOR GRS-80
*** ESQ IS THE SQUARE OF THE 1ST ECCENTRICITY
*** E IS THE 1ST ECCENTRICITY
*** CM IS THE CENTRAL MERIDIAN OF THE PROJECTION ZONE
*** FE IS THE FALSE EASTING VALUE AT THE CM
*** CONV IS CONVERGENCE
*** EPS IS THE SQUARE OF THE 2ND ECCENTRICITY
*** SF IS THE SCALE FACTOR AT THE CM
*** SO IS THE MERIDIANAL DISTANCE (TIMES THE SF) FROM THE
***       EQUATOR TO SOUTHERNMOST PARALLEL OF LAT. FOR THE ZONE
*** R IS THE RADIUS OF THE RECTIFYING SPHERE
*** U0,U2,U4,U6,V0,V2,V4,V6 ARE PRECOMPUTED CONSTANTS FOR
***   DETERMINATION OF MERIDIANAL DIST. FROM LATITUDE
***
*** THE FORMULA USED IN THIS SUBROUTINE GIVES GEODETIC ACCURACY
*** WITHIN ZONES OF 7 DEGREES IN EAST-WEST EXTENT.  WITHIN STATE
*** TRANSVERSE MERCATOR PROJECTION ZONES, SEVERAL MINOR TERMS OF
*** THE EQUATIONS MAY BE OMMITTED (SEE A SEPARATE NGS PUBLICATION).
*** IF PROGRAMMED IN FULL, THE SUBROUTINE CAN BE USED FOR
*** COMPUTATIONS IN SURVEYS EXTENDING OVER TWO ZONES.
***********************************************************************
 
	IMPLICIT DOUBLE PRECISION(A-H,K-Z)
 
	OM=(N-FN+SO)/(R*SF)
	COSOM=DCOS(OM)
	FOOT=OM+DSIN(OM)*COSOM*(V0+V2*COSOM*COSOM+V4*COSOM**4+
     &     V6*COSOM**6)
	SINF=DSIN(FOOT)
	COSF=DCOS(FOOT)
	TN=SINF/COSF
	TS=TN*TN
	ETS=EPS*COSF*COSF
	RN=ER*SF/DSQRT(1.D0-ESQ*SINF*SINF)
	Q=(E-FE)/RN
	QS=Q*Q
	B2=-TN*(1.D0+ETS)/2.D0
	B4=-(5.D0+3.D0*TS+ETS*(1.D0-9.D0*TS)-4.D0*ETS*ETS)/12.D0
	B6=(61.D0+45.D0*TS*(2.D0+TS)+ETS*(46.D0-252.D0*TS-
     &    60.D0*TS*TS))/360.D0
	B1=1.D0
	B3=-(1.D0+TS+TS+ETS)/6.D0
	B5=(5.D0+TS*(28.D0+24.D0*TS)+ETS*(6.D0+8.D0*TS))/120.D0
	B7=-(61.D0+662.D0*TS+1320.D0*TS*TS+720.D0*TS**3)/5040.D0
	LAT=FOOT+B2*QS*(1.D0+QS*(B4+B6*QS))
	L=B1*Q*(1.D0+QS*(B3+QS*(B5+B7*QS)))
	LON=-L/COSF+CM

C*********************************************************************
C     COMPUTE CONVERENCE AND SCALE FACTOR 
	FI=LAT
	LAM = LON
	SINFI=SIN(FI)
	COSFI=COS(FI)
	L1=(LAM-CM)*COSFI
	LS=L1*L1
C
C*** CONVERGENCE
	C1=-TN
	C3=(1.+3.*ETS+2.*ETS**2)/3.
	C5=(2.-TS)/15.
	CONV=C1*L1*(1.+LS*(C3+C5*LS))
C
C*** POINT SCALE FACTOR
	F2=(1.+ETS)/2.
	F4=(5.-4.*TS+ETS*( 9.-24.*TS))/12.
	KP=SF*(1.+F2*LS*(1.+F4*LS))
 
	RETURN
	END
