SUBROUTINE TMPLXY(PHI,DLAM,A,B,SF,XO,CMRAD,X,Y)                           
C                                                                               
C                                                                               
C                                                                               
C        SUBROUTINE TMPLXY COMPUTES THE X,Y COORDINATES FOR THE                 
C      TRANSVERSE MERCATOR PROJECTION GIVEN THE GEOGRAPHIC COORDINATES          
C      -LATITUDE AND LONGITUDE. THE EQUATIONS USED TO COMPUTE X AND Y           
C      ARE FROM THOMAS (1952). SUBROUTINE MERARC IS USED TO COMPUTE             
C      THE MERIDIAN ARC LENGTH.                                                 
C                                                                               
C                                                                               
C        INPUT:                                                                 
C                                                                               
C              PHI -LATITUDE IN RADIANS                                         
C                                                                               
C              DLAM-LONGITUDE OF POINT MINUS LONGITUDE OF CENTRAL               
C                   MERIDIAN (IN RADIANS) FOR LONGITUDE POSITIVE EAST.          
C                                                                               
C              A   -SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID.                 
C                                                                               
C              B   -SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID.                 
C                                                                               
C              XO - FALSE EASTING OF THE CENTRAL MERIDIAN.                      
C                                                                               
C              SF - SCALE OF THE CENTRAL MERIDIAN.                              
C                                                                               
C              CMRAD - THE CENTRAL MERIDIAN,IN RADIANS.                         
C                                                                               
C        OUTPUT:                                                                
C                                                                               
C              X   -EASTING COORDINATE OF THE TRANSVERSE MERCATOR               
C                   PROJECTION.                                                 
C                                                                               
C              Y   -NORTHING COORDINATE OF THE TRANSVERSE MERCATOR              
C                   PROJECTION.                                                 
C                                                                               
C                                  WRITTEN BY R.R.STEEVES AND G.BOWIE           
C                                             MAY,1977                          
C                                                                               
C                                                                               
      IMPLICIT REAL*8(A-H,O-Z)                                                  
      SP=DSIN(PHI)                                                              
      CP=DCOS(PHI)                                                              
      T=DTAN(PHI)                                                               
      E=DSQRT((A**2-B**2)/A**2)                                                 
      ETA=DSQRT((A**2-B**2)/B**2*CP**2)                                         
      CALL MERARC(PHI,A,B,SPHI)                                                 
      DN=A/DSQRT(1.0D0-E**2*SP**2)                                              
      X=DN*(DLAM*CP+DLAM**3*CP**3/6.0D0*(1.0D0-T**2+ETA**2)+DLAM**5*CP          
     1   **5/120.0D0*(5.0D0-18.0D0*T**2+T**4+14.0D0*ETA**2-58.0D0*T**2*         
     2   ETA**2+13.D0*ETA**4+4.D0*ETA**6-64.D0*ETA**4*T**2-24.D0*ETA**6         
     3   *T**2) +DLAM**7/5040.0D0*CP**7*(61.0D0-479.0D0*T**2+179.0D0*           
     4   T**4-T**6))                                                            
      Y=SPHI+DN*   (DLAM**2/2.0D0*SP*CP+DLAM**4/24.0D0*SP*CP**3*(5.0D0-         
     1   T**2+9.0D0*ETA**2+4.0D0*ETA**4)+DLAM**6/720.0D0*SP*CP**5*              
     2   (61.0D0-58.0D0*T**2+T**4+270.0D0*ETA**2-330.0D0*T**2*ETA**2            
     3   +445.D0*ETA**4+324.D0*ETA**6-680.D0*ETA**4*T**2+88.D0*ETA**8           
     4   -600.D0*ETA**6*T**2-192.D0*ETA**8*T**2)+DLAM**8/40320.D0*SP*CP*        
     5   *7*(1385.D0-3111.D0*T**2+543.D0*T**4-T**6))                            
      X=SF*X+XO                                                                 
      Y=SF*Y                                                                    
      RETURN                                                                    
      END