SUBROUTINE PLHXYZ(PHI,RLAM,H,XO,YO,ZO,A,B,X,Y,Z)                          
C                                                                               
C        THIS ROUTINE COMPUTES THE CARTESIAN COORDINATES X,Y,Z GIVEN THE        
C      ELLIPSOIDAL COORDINATES PHI,RLAM,H.                                      
C                                                                               
C        INPUT:                                                                 
C              PHI-ELLIPSOIDAL LATITUDE IN RADIANS.                             
C              RLAM-ELLIPSOIDAL LONGITUDE IN RADIANS.                           
C                   (POSITIVE EAST OF GREENWICH)                                
C              H-ELLIPSOIDAL HEIGHT 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              X,Y,Z-CARTESIAN COORDINATES OF THE POINT 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)                                                        
      SP=DSIN(PHI)                                                              
      CP=DCOS(PHI)                                                              
      N=A/DSQRT(1.D0-E2*SP**2)                                                  
      X=XO+(N+H)*CP*DCOS(RLAM)                                                  
      Y=YO+(N+H)*CP*DSIN(RLAM)                                                  
      Z=ZO+(N*(1.D0-E2)+H)*SP                                                   
      RETURN                                                                    
      END