C     $Id: convrt.f 3304 2011-01-17 15:25:59Z brideout $
C
      SUBROUTINE CONVRT(I,GDLAT,GDALT,GCLAT,RKM)
C
C     jmh - 11/79  ans fortran 66
C
C     Converts between geodetic and geocentric coordinates. the
C     reference geoid is that adopted by the iau in 1964. a=6378.16,
C     b=6356.7746, f=1/298.25. the equations for conversion from
C     geocentric to geodetic are from astron. j., vol 66, 1961, p. 15.
C
C       Input:
C             I - 1 geodetic to geocentric, 2 geocentric to geodetic
C       Input, Output:
C         GDLAT - geodetic latitude (degrees)
C         GDALT - altitude above geoid (km)
C         GCLAT - geocentric latitude (degrees)
C           RKM - geocentric radial distance (km)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION GCLAT,GDALT,GDLAT,RKM
      INTEGER I
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION A,A2,A4,A6,A8,AB2,C2CL,C4CL,CCL,CL2,COSBET,
     *                 COSLAT,DLTCL,DTR,EP2,GCL,GDL,RER,RGEOID,S2CL,
     *                 S4CL,S6CL,S8CL,SB2,SCL,SINBET,SINLAT,SL2,X,Y
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC DATAN2,DCOS,DMIN1,DSIN,DSQRT
C     ..
C     .. Data statements ..
C
      DATA A/6378.16D0/,AB2/1.0067397D0/,EP2/.0067397D0/
      DATA DTR/.0174532925199D0/
C     ..
C
      IF (I.EQ.2) THEN
C
C        .....geocentic to geodetic.....
         RER = RKM/A
         A2 = ((-1.4127348D-8/RER+.94339131D-8)/RER+.33523288D-2)/RER
         A4 = (((-1.2545063D-10/RER+.11760996D-9)/RER+.11238084D-4)/RER-
     *        .2814244D-5)/RER
         A6 = ((54.939685D-9/RER-28.301730D-9)/RER+3.5435979D-9)/RER
         A8 = (((320.D0/RER-252.D0)/RER+64.D0)/RER-5.D0)/RER*
     *        .98008304D-12
         GCL = DTR*GCLAT
         CCL = DCOS(GCL)
         SCL = DSIN(GCL)
         S2CL = 2.D0*SCL*CCL
         C2CL = 2.D0*CCL*CCL - 1.D0
         S4CL = 2.D0*S2CL*C2CL
         C4CL = 2.D0*C2CL*C2CL - 1.D0
         S8CL = 2.D0*S4CL*C4CL
         S6CL = S2CL*C4CL + C2CL*S4CL
         DLTCL = S2CL*A2 + S4CL*A4 + S6CL*A6 + S8CL*A8
         GDLAT = GCLAT + DLTCL/DTR
         GDALT = RKM - A/DSQRT(1.D0+EP2*SCL*SCL)
      ELSE
C
C        .....geodetic to geocentric.....
         GDL = DTR*GDLAT
         SINLAT = DSIN(GDL)
         COSLAT = DCOS(GDL)
         SL2 = SINLAT*SINLAT
         CL2 = AB2*COSLAT
         CL2 = CL2*CL2
         SINBET = SINLAT/DSQRT(SL2+CL2)
         SB2 = DMIN1(SINBET*SINBET,1.D0)
         COSBET = DSQRT(1.D0-SB2)
         RGEOID = A/DSQRT(1.D0+EP2*SB2)
         X = RGEOID*COSBET + GDALT*COSLAT
         Y = RGEOID*SINBET + GDALT*SINLAT
         RKM = DSQRT(X*X+Y*Y)
         GCLAT = DATAN2(Y,X)/DTR
      END IF
      RETURN
C
      END