C*********************************************************************
C  PROGRAM TO CONVERT FROM LATITUDE AND LONGITUDE IN DEG-MIN-SEC
C  TO UTM COORDINATES (IN METERS).  SUBROUTINES WERE PROVIDED BY
C  GARTH D. NEWTON AND ARE DESCRIBED IN U.S. GEOLOGICAL SURVEY
C  BULLETIN 1642, 1985.
C
C  NOTE: THE REQUIRED USE OF DOUBLE PRECISION
C        WRITTEN BY WRDANSKIN@DCASCR
C        LAST REVISION     1/6/87
C*********************************************************************
      PROGRAM MAIN
      IMPLICIT DOUBLE PRECISION (A-Z)
      INTEGER LATRED,LONRED,J
      DIMENSION LAT(3),LON(3),LATRED(3),LONRED(3)
C  INITIALIZE DATA DEPENDENT UPON THE CHOICE OF ORIGIN, E.G. SETS ZONE
      LAT(1)=37.
      LAT(2)=0.
      LAT(3)=0.
      LON(1)=118.
      LON(2)=0.
      LON(3)=0.
C  CONVERT TO DECIMAL DEGREES
      CALL DEGREE(LAT,ALAT)
      CALL DEGREE(LON,ALON)
C  SET GLOBAL PARAMETERS
      CALL DATA(ALAT,ALON)
C  INPUT LATITUDE AND LONGITUDE IN DEGREES-MINUTES-SECONDS
C  AND COMPUTE GRID COORDINATES.  THEN COMPUTE NEW LATITUDE
C  AND LONGITUDE FROM GRID COORDINATES.  THEN COMPUTE NEW
C  GRID COORDINATES AND COMPARE ALL RESULTS.  PRINT THE
C  FINAL ONE
C  USE A DO-LOOP OF NCOUNT LAT/LONG PAIRS
      READ(5,*)NCOUNT
      WRITE(6,900)
      DO 100 I=1,NCOUNT
      READ(5,*)LAT,LON
      IF(LAT(1).EQ.0.) GOTO 200
C  SAVE ORIGINAL VALUES FOR FINAL PRINTOUT
      DO 10 J=1,3
      LATRED(J)=INT(LAT(J))
      LONRED(J)=INT(LON(J))
   10 CONTINUE
C  CONVERT FROM DEG,MIN,SEC TO DEGREES
      CALL DEGREE(LAT,ALAT)
      CALL DEGREE(LON,ALON)
C  COMPUTE UTM GRID COORDINATES
      CALL XYTRAN(ALAT,ALON,XUTM,YUTM)
C  PRINT RESULTS
C      WRITE(6,901)LATRED,LONRED,XUTM,YUTM
      WRITE(6,902)LATRED,LONRED,XUTM,YUTM
  100 CONTINUE
  200 CONTINUE
C  FORMAT STATEMENTS
  900 FORMAT('   LAT      LONG       XUTM      YUTM')
  901 FORMAT(1X,3I3,1X,3I3,2F10.0)
  902 FORMAT(1X,3I3,1X,3I3,2F15.4)
      STOP
      END
C********************************************************************
C      CONVERT DEG-MIN-SEC TO DEGREES
C
      SUBROUTINE DEGREE(A,B)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      DIMENSION A(3)
      B=A(1)+A(2)/60.+A(3)/3600.
      RETURN
      END
C*********************************************************************
C     CONVERT DEGREES TO DEG-MIN-SEC
C
      SUBROUTINE DECIMAL(A,B)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      DIMENSION A(3)
      A(1)=AINT(B+.5/3600.)
      A(2)=AINT((B-A(1))*60.+.5/60.)
      A(3)=AINT((B-A(1)-A(2)/60.)*3600.+.5)
      RETURN
      END
C***********************************************************************
C      COMPUTE THE LATITUDE AND LONGITUDE GIVEN THE
C      RECTANGULAR COORDINATES X AND Y
C
C1     CHANGED AUTMI TO INVERSE 01/09/85
      SUBROUTINE INVERSE(LAT,LON,X,Y)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0
C
      X=(X-500000.)/M0
      Y=Y/M0
C
C----- ITERATE TO COMPUTE MERIDIONAL DISTANCE
C     MAKE INITIAL GUESS
      PHI=SINV(Y,PHI0)
      OLD=Y
C BEGIN ITERATIONS
      DO 10 I=1,10
      PHI=SINV(Y,PHI)
C     IF OLD VALUE EQUALS NEW VALUE STOP
      IF(OLD.EQ.Y0) GOTO 11
      OLD=Y0
10    CONTINUE
11    S=DSIN(PHI)
      C=DCOS(PHI)
      T=S/C
      T2=T*T
      N2=EPS*C*C
      N=A/SQRT(1.D0-E2*S*S)
      AA0=X/N
      AA1=AA0**3*(1.D0+2.D0*T2+N2)/6.D0
      AA2=AA0**5*(5.D0+28.D0*T2+24.D0*T2*T2+6.D0*N2+8.D0*T2*N2)
     #/120.D0
C
C      COMPUTE LONGITUDE (RADIANS)
CC
      LON=LAM0-(1.D0/C)*(AA0-AA1+AA2)
C
      T3=T*T2
      T4=T*T3
      T5=T*T4
      T6=T*T5
      N4=N2*N2
      N6=N4*N2
      N8=N6*N2
      AA1=T*(1.D0+N2)*AA0*AA0
      AA2=T*(1.D0+N2)*AA0**4*(5.D0+3.D0*T2+N2-4.D0*N4-9.D0*N2*T2)
      AA3=T*(1.D0+N2)*AA0**6*(61.D0+90.D0*T2+46.D0*N2+45.D0*T4
     # -252.D0*T2*N2
     #-3.D0*N4+100*N6-66*T2*N4-90.D0*T4*N2+88.D0*N8
     #+225.D0*T4*N4+84.D0*T2*N6-192.D0*T2*N8)
      AA4=T*(1.D0+N2)*AA0**8*(1385.D0+3633.D0*T2+4095.D0*T4+1575.D0*T6)
C
C    COMPUTE LATITUDE (RADIANS)
C
      LAT=PHI-AA1/2.D0+AA2/24.D0-AA3/720.D0+AA4/40320.D0
C
C     CONVERT TO DEGREES
      LAT=LAT /D2R
      LON=LON /D2R
C
      RETURN
      END
C*********************************************************************
C      INITIALIZE CONSTANTS RELATED TO THE LATITUDE AND LONGITUDE
C      OF THE AREA BEING MAPPED
C
C1     CHANGE FROM UTM TO DATA 01/09/85
      SUBROUTINE DATA(BASLAT,LAM)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0
      A=6378206.4
      B=6356584.8
      M0=0.9996
      D2R=.0174533
      E2=(A*A-B*B)/(A*A)
      EPS=(A*A-B*B)/(B*B)
      E4=E2*E2
      E6=E4*E2
      E8=E6*E2
      AE=A*(1.0-E2)
      A0=AE*(1.+(3./4.)*E2+(45./64.)*E4+(175./256.)*E6
     1 +(11025./16384.)*E8)
      A1=-AE/2.*((3./4.)*E2+(15./16)*E4+(525./512.)*E6
     1 +(2205./2048.)*E8)
      A2=AE/4.*((15./64.)*E4+(105./256.)*E6+(2205./4096.)*E8)
      ERR=-AE/6.*((35./512.)*E6+(315./2048.)*E8)
C      PRINT *,'A-CONSTANTS ',A0,A1*2,A2*2,ERR*2
      PHI0=BASLAT*D2R
      IZONE=AINT((180.-LAM)/6.+1.)
      LAM0=(183.-(6.*FLOAT(IZONE)))*D2R
C      PRINT *,'ZONE NUMBER IS ',IZONE
C      PRINT *,'CENTRAL MERIDIAN IS ',LAM0/D2R
      RETURN
      END
C*********************************************************************
C     CALCULATE THE DISTANCE FROM THE EQUATOR ALONG THE MERIDION
      FUNCTION SS(PHI)
      IMPLICIT DOUBLE PRECISION (A-Z)
C--------------------------------------------------------------------
      COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0
      LAT=PHI*D2R
      SS=(A0*LAT + A1*DSIN(2.D0*LAT) + A2*DSIN(4.D0*LAT) +
     * ERR*DSIN(6.D0*LAT))
      RETURN
      END
C***********************************************************************
C      COMPUTE THE RECTANGULAR COORDINATES X AND Y FROM THE
C      LATITUDE AND LONGITUDE
C
C1     CHANGED FROM UTMF TO XYTRAN 01/09/85
      SUBROUTINE XYTRAN(LAT,LON,X,Y)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0
      D=LAM0-LON*D2R
      D2=D*D
      D3=D2*D
      D4=D3*D
      D5=D4*D
      D6=D5*D
      D7=D6*D
      D8=D7*D
C
      S=DSIN(LAT*D2R)
      C=DCOS(LAT*D2R)
      T=S/C
      T2=T*T
      T4=T2*T2
      T6=T4*T2
      C2=C*C
      C3=C2*C
      C4=C3*C
      C5=C4*C
      C6=C5*C
      C7=C6*C
C
      N2=EPS*C*C
      N4=N2*N2
      N=A/SQRT(1.-E2*S*S)
      AA0=N*C
      AA1=N*C3*(1.-T2+N2)
      AA2=N*C5*(5.-18.*T2+T4+14.*N2
     1-58.*T2*N2)
      AA3=N*S*C7*(61.-479.*T2+179.*T4-T6)
      B0=SS(LAT)
      B1=N*S*C
      B2=N*S*C3*(5.-T2+9.*N2+4.*N4)
      B3=N*S*C5*(61.-58.*T2+T4+270.*N2
     1-330.*T2*N2)
      B4=N*S*C7*(1385.-3111.*T2+543.*T4-T6)
C
      X=M0*(AA0*D+AA1*D3/6.+AA2*D5/120.+AA3*D7/5040.)+500000.
      Y=M0*(B0+B1*D2/2.+B2*D4/24.+B3*D6/720.+B4*D8/40320)
      RETURN
      END
C*********************************************************************
C      FUNCTION TO COMPUTE THE INVERSE OF FUNCTION SS.  IT FINDS THE
C      LATITUDE OF A POINT A DISTANCE Y FROM THE EQUATOR.  PHI IS A
C      INITIAL GUESS AT THE LATITUDE OR A POINT NEAR Y
C
      FUNCTION SINV(Y,PHI)
      IMPLICIT DOUBLE PRECISION (A-Z)
C
C---------------------------------------------------------------------
      COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0
      Y0=SS(PHI/D2R)
      DY=Y-Y0
      C=A*(1.-E2)
      D=(1.-E2*DSIN(PHI)*DSIN(PHI))
      C1=-3.*E2/(C*C)
      P1=D**(1.5)/C
      P2=C1*D*D*DCOS(PHI)
      SINV=(PHI+DY*P1+DY*DY*P2)
      RETURN
      END

