SUBROUTINE TMXYPL(X,Y,A,B,SF,XO,CMRAD,PHI,OLAM)                           
C                                                                               
C                                                                               
C        SUBROUTINE TMXYPL COMPUTES THE GEOGRAPHIC COORDINATES- LATITUDE        
C      AND LONGITUDE - GIVEN THE X,Y COORDINATES OF THE TRANSVERSE              
C      MERCATOR PROJECTION. THE EQUATIONS USED TO COMPUTE THE LONGITUDE         
C      AND LATITUDE ARE FROM THOMAS (1952). SUBROUTINE FPLAT IS USED            
C      TO COMPUTE THE FOOT-POINT LATITUDE.                                      
C                                                                               
C                                                                               
C        INPUT;                                                                 
C                                                                               
C              X   -EASTING COORDINATE OF THE TRANSVERSE MERCATOR               
C                   PROJECTION.                                                 
C                                                                               
C              Y   -NORTHING COORDINATE OF THE TRANSVERSE MERCATOR              
C                   PROJECTION.                                                 
C                                                                               
C              A   -SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID.                 
C                                                                               
C              B   -SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID.                 
C                                                                               
C              SF - SCALE OF THE CENTRAL MERIDIAN.                              
C                                                                               
C              XO - FALSE EASTING OF THE CENTRAL MERIDIAN.                      
C                                                                               
C              CMRAD - THE CENTRAL MERIDIAN,IN RADIANS.                         
C                                                                               
C        OUTPUT:                                                                
C                                                                               
C              PHI -LATITUDE OF THE POINT IN RADIANS                            
C                                                                               
C              OLAM-LONGITUDE OF THE POINT IN RADIANS                           
C                                                                               
C                                              WRITTEN BY R.R.STEEVES.          
C                                                    MAY,1977                   
C                                                                               
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      X=(X-XO)/SF                                                               
      Y=Y/SF                                                                    
      E=DSQRT((A**2-B**2)/A**2)                                                 
      CALL FPLAT(A,B,Y,PHI1)                                                    
      T=DTAN(PHI1)                                                              
      SP=DSIN(PHI1)                                                             
      CP=DCOS(PHI1)                                                             
      ETA=DSQRT((A**2-B**2)/B**2*CP**2)                                         
      DN=A/DSQRT(1.0D0-E**2*SP**2)                                              
      DM=A*(1.0D0-E**2)/DSQRT((1.0D0-E**2*SP**2)**3)                            
      PHI=PHI1-T*X**2/2.0D0/DM/DN+T*X**4/24.0D0/DM/DN**3*(5.0D0+3.0D0*          
     1   T**2+ETA**2-4.0D0*ETA**4-9.0D0*ETA**2*T**2)-T*X**6/720.0D0/DM/         
     2   DN**5*(61.0D0+90.0D0*T**2+46.0D0*ETA**2+45.0D0*T**4-252.0D0*T**        
     3   2*ETA**2-3.0D0*ETA**4+100.0D0*ETA**6-66.0D0*T**2*ETA**4-90.0D0         
     4   *T**4*ETA**2+88.0D0*ETA**8+225.0D0*T**4*ETA**4+84.0D0*T**2*            
     5   ETA**6-192.0D0*T**2*ETA**8)                                            
      PHI=PHI+T*X**8/40320.0D0/DM/DN**7*(1385.0D0+3633.0D0*T**2+4095.0D0        
     1    *T**4+1575.0D0*T**6)                                                  
      DLAM=(X/DN-(X/DN)**3/6.0D0*(1.0D0+2.0D0*T**2+ETA**2)+(X/DN)**5/           
     1   120.0D0*(5.0D0+6.0D0*ETA**2+28.0D0*T**2-3.0D0*ETA**4+8.0D0*T**2        
     2   *ETA**2+24.0D0*T**4-4.0D0*ETA**6+4.0D0*T**2*ETA**4+24.0D0*T**2*        
     3   ETA**6)-(X/DN)**7/5040.0D0*(61.0D0+662.0D0*T**2+1320.0D0*T**4+         
     4   720.0D0*T**6))/CP                                                      
      OLAM=CMRAD+DLAM                                                           
      X=X*SF+XO                                                                 
      Y=Y*SF                                                                    
      RETURN                                                                    
      END