SUBROUTINE XYZPLH(X,Y,Z,XO,YO,ZO,A,B,PHI,RLAM,H)                          
C                                                                               
C        THIS ROUTINE COMPUTES THE ELLIPSOIDAL COORDINATES PHI,RLAM,H           
C      GIVEN THE CARTESION COORDINATES X,Y,Z.                                   
C                                                                               
C        INPUT:                                                                 
C              X,Y,Z-CARTESION COORDINATES OF THE POINT IN METRES.              
C              XO,YO,ZO-TRANSLATION COMPONENTS FROM THE ORIGIN OF THE           
C                       CARTESIAN COORDINATE SYTEM (X,Y,Z)TO THE CENTER         
C                       OF THE REFERENCE ELLIPSOID. (IN METRES.)                
C              A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE              
C                  ELLIPSOID IN METRES.                                         
C                                                                               
C        OUTPUT:                                                                
C              PHI-ELLIPSOIDAL LATITUDE IN RADIANS.                             
C              RLAM-ELLIPSOIDAL LONGITUDE IN RADIANS.                           
C                   (POSITIVE EAST OF GREENWICH)                                
C              H-ELLIPSOIDAL HEIGHT IN METRES.                                  
C                                                                               
C                                           WRITTEN BY R.R.STEEVES              
C                                                JUNE,1977                      
C                                                                               
C                                                                               
      IMPLICIT REAL*8(A-Z)                                                      
      E2=(A*A-B*B)/(A*A)                                                        
      XP=X-XO                                                                   
      YP=Y-YO                                                                   
      ZP=Z-ZO                                                                   
      S=DSQRT(XP**2+YP**2)                                                      
      RLAM=DATAN2(YP,XP)                                                        
      ZPS=ZP/S                                                                  
      H=DSQRT(XP**2+YP**2+ZP**2)-A                                              
      PHI=DATAN(ZPS/(1.D0-E2*A/(A+H)))                                          
    1 N=A/DSQRT(1.D0-E2*DSIN(PHI)**2)                                           
      HP=H                                                                      
      PHIP=PHI                                                                  
      H=S/ DCOS(PHI)-N                                                          
      PHI=DATAN(ZPS/(1.D0-E2*N/(N+H)))                                          
      IF(DABS(PHIP-PHI).GT.1.D-11.OR.DABS(HP-H).GT.1.D-5)GO TO 1                
      RETURN                                                                    
      END