SUBROUTINE VININ (AE,F,XLAT1,XLONG1,XLAT2,XLONG2,DIST,AZ1)                
      IMPLICIT REAL * 8 (A - H, O - Z)                                          
C  STATEMENT FUNCTIONS:                                                         
      SIN(A) = DSIN(A)                                                          
      COS(A) = DCOS(A)                                                          
      ACOS(A)=DARCOS(A)                                                         
      TAN(A) = DSIN(A)/DCOS(A)                                                  
      ATAN(A) = DATAN(A)                                                        
      ATAN2(A,B)=DATAN2(A,B)                                                    
      SQRT(A) = DSQRT(A)                                                        
      ABS(A) = DABS(A)                                                          
      PI=ACOS(-1.D0)                                                            
C$                                                                              
C NAME        VININ                                                             
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     INVERSE GEODETIC PROBLEM                                          
C                                                                               
C AUTHOR/DATE P.DELORME 780830 USING ALGORITHM OF T.VINCENTY                    
C                                                                               
C EXTERNALS   DSIN,DCOS,DARCOS,DATAN,DATN2,DSQRT,DABS                           
C                                                                               
C CALLING     CALL VININ(AE,F,XLAT1,XLONG1,XLAT2,XLONG2,DIST,AZ1)               
C                                                                               
C PARAMETERS  AE = ELLIPSOID SEMI-MAJOR AXIS (M)                                
C             F = INVERSE OF ELLIPSOID FLATTENING                               
C             XLAT1 = INPUT LAT OF POINT 1 (DEG)                                
C             XLONG1 = INPUT LONG OF POINT 1 (DEG)                              
C             XLAT2 = INPUT LAT OF POINT 2 (DEG)                                
C             XLONG2 = INPUT LONG OF POINT 2 (DEG)                              
C             DIST = OUTPUT GEODESIC DISTANCE (M)                               
C             AZ1 = OUTPUT FORWARD AZIMUTH (DEG)                                
C             AZ2 = OUTPUT BACK AZIMUTH (DEG) (IF REQUIRED ADD AZ2              
C                   TO END OF ARGUMENT LIST AND ACTIVATE APPENDIX)              
C                                                                               
C LANGUAGE    FORTRAN                                                           
C$                                                                              
C                                                                               
      FL = 1D0 / F                                                              
      FUZ = 20.0D0 * DEPS(DARG)                                                 
      BE = AE *(1.0 - FL)                                                       
C                                                                               
C   CONVERT LATITUDE AND LONGITUDE OF POINTS 1 AND 2 TO RADIANS.                
C                                                                               
      ALAT1 = (XLAT1 * PI) / 180D0                                              
      ALON1 = (XLONG1 * PI) / 180D0                                             
      ALAT2 = (XLAT2 * PI) / 180D0                                              
      ALON2 = (XLONG2 * PI) / 180D0                                             
C                                                                               
C   CALCULATE THE REDUCED LATITUDE AND ITS TRIG FUNCTIONS:                      
C                                                                               
      TU1 = (1D0 - FL) * TAN(ALAT1)                                             
      TU2 = (1D0 - FL) * TAN(ALAT2)                                             
      U1 = ATAN(TU1)                                                            
      U2 = ATAN(TU2)                                                            
      SU1 = SIN(U1)                                                             
      SU2 = SIN(U2)                                                             
      CU1 = COS(U1)                                                             
      CU2 = COS(U2)                                                             
C                                                                               
C   1ST APPROX:  DIFF IN LONG ON ELLIPSOID = DIFF IN LONG ON SPHERE:            
C                                                                               
      DL = ALON2 - ALON1                                                        
      XDL = DL                                                                  
C                                                                               
C                                                                               
C   *********************  BEGINNING OF ITERATION  *****************            
C                                                                               
C                                                                               
100   CONTINUE                                                                  
      SDL = SIN(DL)                                                             
      CDL = COS(DL)                                                             
      CS = (SU1 * SU2) + (CU1 * CU2 * CDL)                                      
      SS = SQRT(1D0 - (CS**2))                                                  
      SIG = ATAN2(SS , CS)                                                      
      IF(ABS(SS) .LT. FUZ)  SS = FUZ                                            
      SA = (CU1 * CU2 * SDL)/SS                                                 
      CA = SQRT(1D0 - (SA**2))                                                  
      C2SM = CS - (2D0 * SU1 * SU2) / (CA**2)                                   
C                                                                               
C   CALCULATE THE DIFFERENCE IN LONGITUDE ON THE AUXILIARY SPHERE:              
C                                                                               
      C = (FL / 16D0) * (CA**2) * (4D0 + (FL * (4D0 - 3D0 * (CA**2))))          
      DL1 = XDL + (1D0 - C) * FL * SA * (SIG + C * SS * (C2SM + C               
     $                      * CS * (-1D0 + 2D0 * (C2SM**2))))                   
      IF(ABS(DL1 - DL) .LE. 10D-10)  GOTO 45                                    
      DL = DL1                                                                  
      GOTO 100                                                                  
45    CONTINUE                                                                  
C                                                                               
C                                                                               
C   *********************  END OF ITERATION  *********************              
C                                                                               
C                                                                               
C                                                                               
C   CALCULATE U SQUARE:                                                         
C                                                                               
      U = (CA**2) * ((AE**2) - (BE**2)) / (BE**2)                               
C                                                                               
C   CALCULATE A:                                                                
C                                                                               
      A = 1D0 + (U / 256D0) * (64D0 + U * (-12D0 + 5D0 *U))                     
C                                                                               
C   CALCULATE B:                                                                
C                                                                               
      B = (U / 512D0) * (128D0 + U * (-64D0 + (37D0 * U)))                      
C                                                                               
C   CALCULATE DELTA SIGMA:                                                      
C                                                                               
      DSIG = B * SS * (C2SM + 0.25D0 * B * CS * (-1D0 +2D0*(C2SM**2)))          
C                                                                               
C   CALCULATE THE GEODESIC:                                                     
C                                                                               
      DIST = BE * A * (SIG - DSIG)                                              
C                                                                               
C   CALCULATE THE FORWARD AZIMUTH:                                              
C                                                                               
      SDL1 = SIN(DL1)                                                           
      CDL1 = COS(DL1)                                                           
      AZ1 = ATAN2((CU2 * SDL1),(CU1 * SU2 - SU1 * CU2 * CDL1))                  
     $           * 180D0 / PI                                                   
      IF (AZ1 .GT. 0D0)  GOTO 70                                                
      AZ1 = AZ1 + 360D0                                                         
70    CONTINUE                                                                  
C                                                                               
C                                                                               
C  APPENDIX   **************  CALCULATE THE BACKWARD AZIMUTH:  *****            
C                                                                               
C     AZ2 = ATAN2((-1D0 * CU1 * SDL1),(SU1 * CU2 - CU1 * SU2 * CDL1             
C    $           * 180D0 / PI                                                   
C     IF(AZ2 .GT. 0D0)  GOTO 75                                                 
C     AZ2 = AZ2 + 360D0                                                         
C75   CONTINUE                                                                  
C                                                                               
C                                                                               
C  *******************************************************************          
C                                                                               
C                                                                               
C                                                                               
      RETURN                                                                    
      END