SUBROUTINE MERARC(PHI,A,B,S)                                              
C                                                                               
C                                                                               
C       THIS ROUTINE COMPUTES THE MERIDIAN ARC LENGTH FROM THE EQUATOR          
C        TO LATITUDE PHI ON AN ELLIPSOID OF REVOLUTION DEFINED BY ITS           
C        SEMI-MAJOR AXIS A AND ITS SEMI-MINOR AXIS B. THE COMPUTED ARC          
C        LENGTH IS ACCURATE TO APPROXIMATELY 10 MICROMETRES OVER THE            
C        ENTIRE RANGE: EQUATOR TO POLE.                                         
C                                                                               
C        INPUT: PHI-ELLIPSOIDAL LATITUDE IN RADIANS.(MAY BE POSITIVE            
C                   (NORTH) OR NEGATIVE(SOUTH OF EQUATOR))                      
C                                                                               
C               A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE ELLIPSOID.            
C                   (THE COMPUTED ARC LENGTH WILL BE OUTPUT IN THE SAME         
C                   UNITS AS A AND B).                                          
C                                                                               
C        OUTPUT: S -MERIDIAN ARC LENGTH                                         
C                                                                               
C                                                                               
C                                              WRITTEN BY R.R.STEEVES.          
C                                                    MAY,1977                   
C                                                                               
      IMPLICIT REAL*8(A-Z)                                                      
      E2=(A*A-B*B)/(A*A)                                                        
      E4=E2*E2                                                                  
      E6=E4*E2                                                                  
      E8=E6*E2                                                                  
      A0=1.0D0-E2/4.D0-3.D0*E4/64.D0-5.D0*E6/256.D0-175.D0*E8/16384.D0          
      A2=3.D0/8.D0*(E2+E4/4.D0+15.D0*E6/128.D0-455.D0*E8/4096.D0)               
      A4=15.D0/256.D0*(E4+3.D0*E6/4.D0-77.D0*E8/128.D0)                         
      A6=35.D0/3072.D0*(E6-41.D0*E8/32.D0)                                      
      A8=-315.D0*E8/131072.D0                                                   
      S=A*(A0*PHI-A2*DSIN(2.D0*PHI)+A4*DSIN(4.D0*PHI)-A6*DSIN(6.D0*PHI)         
     1  +A8*DSIN(8.D0*PHI))                                                     
      RETURN                                                                    
      END