      SUBROUTINE GCONV2 (A,F,X,Y,Z,PHI,LAMBDA,H)
*+
*  - - - - - - - - - - -
*   G C O N V 2
*  - - - - - - - - - - -
*
*  This routine is part of the International Earth Rotation and
*  Reference Systems Service (IERS) Conventions software collection.
*
*  This subroutine transforms from geocentric rectangular to 
*  geodetic coordinates.
*
*  In general, Class 1, 2, and 3 models represent physical effects that
*  act on geodetic parameters while canonical models provide lower-level
*  representations or basic computations that are used by Class 1, 2, or
*  3 models.
* 
*  Status:  Class 3 model
*
*     Class 1 models are those recommended to be used a priori in the
*     reduction of raw space geodetic data in order to determine
*     geodetic parameter estimates.
*     Class 2 models are those that eliminate an observational
*     singularity and are purely conventional in nature.
*     Class 3 models are those that are not required as either Class
*     1 or 2.
*     Canonical models are accepted as is and cannot be classified as
*     a Class 1, 2, or 3 model.
*
*  Given:
*     a           d         Equatorial Radius of the Earth (Note 1)
*     f           d         Flattening form factor (Note 2)
*     x           d         Rectangular X coordinate (Note 3)
*     y           d         Rectangular Y coordinate (Note 3)
*     z           d         Rectangular Z coordinate (Note 3)
*
*  Returned:
*     phi         d         Latitude coordinate on the ellipsoid (Note 4)
*     lambda      d         Longitude coordinate on the ellipsoid (Note 4)
*     h           d         Height coordinate on the ellipsoid (Note 3)
*
*  Notes:
*
*  1) The parameter given is from the 1980 Geodetic Reference System, which
*     was adopted at the XVII General Assembly of the International Union
*     of Geodesy and Geophysics (IUGG).  It is expressed in meters.
*
*  2) The parameter given is from the 1980 Geodetic Reference System, which
*     was adopted at the XVII General Assembly of the International Union
*     of Geodesy and Geophysics (IUGG).  It is a dimensionless quantity.
*
*  3) Coordinates are expressed in meters.
*
*  4) Coordinates are expressed in radians.
*
*  5) This routine is closely based on the GCONV2H subroutine by
*     Toshio Fukushima (See reference 1).
*
*  6) This version of the routine uses the GRS 1980 ellipsoid parameters
*     as the given default. The user may choose to use other ellipsoids as input parameters.
*
*  Test case:
*     given input: x = 4075579.496D0 meters  Wettzell (TIGO) station
*                  y =  931853.192D0 meters
*                  z = 4801569.002D0 meters
*
*     expected output: phi    =   0.857728298603D0 radians
*                      lambda =   0.224779294628D0 radians
*                      h      = 665.9207D0 meters
*
*  References:
*
*     Fukushima, T., "Transformation from Cartesian to geodetic
*     coordinates accelerated by Halley's method", J. Geodesy (2006),
*     79(12): 689-693
*
*     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
*     IERS Technical Note No. 36, BKG (2010)
*
*  Revisions:
*     2006              T.   Fukushima             Original code
*     2010 March 19     B.E. Stetzler              Added header and copyright
*     2010 March 19     B.E. Stetzler              Initial standardization
*                                                  of code, capitalized
*                                                  variables for FORTRAN
*                                                  77 backwards compatibility
*     2010 March 22      B.E. Stetzler             Provided test case
*     2010 September 2   B.E. Stetzler             Corrected F to match Table
*                                                  1.2 of IERS Conventions
*                                                  (2010) and updated test case
*-----------------------------------------------------------------------

      IMPLICIT NONE
      DOUBLE PRECISION A,F,X,Y,Z,PHI,LAMBDA,H,AEPS2,E2,E4T,EP,EP2,AEP,
     .                 P2,ABSZ,P,S0,PN,ZP,C0,C02,C03,S02,S03,A02,
     .                 A0,A03,D0,F0,B0,S1,CP,S12,CP2

*  Pi
      DOUBLE PRECISION DPI
      PARAMETER ( DPI = 3.141592653589793238462643D0 )

*  -------------
*  Preliminaries
*  -------------

*	Default Values: GRS1980 System (recommended in IERS Conventions)

         A = 6378137D0
         F = 0.00335281068118D0

*  Functions of ellipsoid parameters

         AEPS2 = A*A*1D-32
         E2    = (2.0D0-F)*F
         E4T   = E2*E2*1.5D0
         EP2   = 1.0D0-E2
         EP    = SQRT(EP2)
         AEP   = A*EP

* ---------------------------------------------------------
*	Compute Coefficients of (Modified) Quartic Equation
*
*	Remark: Coefficients are rescaled by dividing by 'a'
* ---------------------------------------------------------

*      Compute distance from polar axis squared

       P2 = X*X + Y*Y

*      Compute longitude lambda 

       IF ( P2.NE.0D0 ) THEN
          LAMBDA = ATAN2(Y,X)
       ELSE
          LAMBDA = 0D0
       END IF

*      Ensure that Z-coordinate is unsigned
    
       ABSZ = ABS(Z)

*      Continue unless at the poles

       IF ( P2.GT.AEPS2 ) THEN

*      Compute distance from polar axis

         P = SQRT(P2)

*      Normalize
       
         S0 = ABSZ/A
         PN = P/A
         ZP = EP*S0

*     Prepare Newton correction factors.

         C0  = EP*PN
         C02 = C0*C0
         C03 = C02*C0
         S02 = S0*S0
         S03 = S02*S0
         A02 = C02+S02
         A0  = SQRT(A02)
         A03 = A02*A0
         D0  = ZP*A03 + E2*S03
         F0  = PN*A03 - E2*C03

*     Prepare Halley correction factor.

         B0 = E4T*S02*C02*PN*(A0-EP)
         S1 = D0*F0 - B0*S0
         CP = EP*(F0*F0-B0*C0)

*     Evaluate latitude and height.

         PHI = ATAN(S1/CP)
         S12 = S1*S1
         CP2 = CP*CP
         H   = (P*CP+ABSZ*S1-A*SQRT(EP2*S12+CP2))/SQRT(S12+CP2)
      ELSE

*     Special case: pole.

         PHI = DPI/2D0
         H   = ABSZ - AEP
      END IF

*  Restore sign of latitude.
      IF ( Z.LT.0D0 ) PHI = -PHI

      RETURN

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2008
*  IERS Conventions Center
*
*  ==================================
*  IERS Conventions Software License
*  ==================================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
*  WHICH APPLY TO ITS USE.
*
*  1. The Software is provided by the IERS Conventions Center ("the
*     Center").
*
*  2. Permission is granted to anyone to use the Software for any
*     purpose, including commercial applications, free of charge,
*     subject to the conditions and restrictions listed below.
*
*  3. You (the user) may adapt the Software and its algorithms for your
*     own purposes and you may distribute the resulting "derived work"
*     to others, provided that the derived work complies with the
*     following requirements:
*
*     a) Your work shall be clearly identified so that it cannot be
*        mistaken for IERS Conventions software and that it has been
*        neither distributed by nor endorsed by the Center.
*
*     b) Your work (including source code) must contain descriptions of
*        how the derived work is based upon and/or differs from the
*        original Software.
*
*     c) The name(s) of all modified routine(s) that you distribute
*        shall be changed.
* 
*     d) The origin of the IERS Conventions components of your derived
*        work must not be misrepresented; you must not claim that you
*        wrote the original Software.
*
*     e) The source code must be included for all routine(s) that you
*        distribute.  This notice must be reproduced intact in any
*        source distribution. 
*
*  4. In any published work produced by the user and which includes
*     results achieved by using the Software, you shall acknowledge
*     that the Software was used in obtaining those results.
*
*  5. The Software is provided to the user "as is" and the Center makes
*     no warranty as to its use or performance.   The Center does not
*     and cannot warrant the performance or results which the user may
*     obtain by using the Software.  The Center makes no warranties,
*     express or implied, as to non-infringement of third party rights,
*     merchantability, or fitness for any particular purpose.  In no
*     event will the Center be liable to the user for any consequential,
*     incidental, or special damages, including any lost profits or lost
*     savings, even if a Center representative has been advised of such
*     damages, or for any claim by any third party.
*
*  Correspondence concerning IERS Conventions software should be
*  addressed as follows:
*
*                     Gerard Petit
*     Internet email: gpetit[at]bipm.org
*     Postal address: IERS Conventions Center
*                     Time, frequency and gravimetry section, BIPM
*                     Pavillon de Breteuil
*                     92312 Sevres  FRANCE
*
*     or
*
*                     Brian Luzum
*     Internet email: brian.luzum[at]usno.navy.mil
*     Postal address: IERS Conventions Center
*                     Earth Orientation Department
*                     3450 Massachusetts Ave, NW
*                     Washington, DC 20392
*
*
*-----------------------------------------------------------------------
      END

