SUBROUTINE W3FB10(DLAT1, DLON1, DLAT2, DLON2, BEARD, GCDKM)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:  W3FB10        LAT/LONG PAIR TO COMPASS BEARING, GCD
C   PRGMMR: CHASE            ORG: NMC421      DATE:88-10-26
C
C ABSTRACT: GIVEN A PAIR OF POINTS (1) AND (2) GIVEN BY LATITUDE AND
C   LONGITUDE, W3FB10 COMPUTES THE BEARING AND GREAT CIRCLE DISTANCE
C   FROM POINT (1) TO POINT (2) ASSUMING A SPHERICAL EARTH.  THE
C   NORTH AND SOUTH POLES ARE SPECIAL CASES.  IF LATITUDE OF POINT
C   (1) IS WITHIN 1E-10 DEGREES OF THE NORTH POLE, BEARING IS THE
C   NEGATIVE LONGITUDE OF POINT (2) BY CONVENTION.  IF LATITUDE OF
C   POINT (1) IS WITHIN 1E-10 DEGREES OF THE SOUTH POLE, BEARING IS
C   THE LONGITUDE OF POINT (2) BY CONVENTION.  IF POINT (2) IS WITHIN
C   1E-6 RADIANS OF THE ANTIPODE OF POINT (1), THE BEARING WILL BE
C   SET TO ZERO.  IF POINT (1) AND POINT (2) ARE WITHIN 1E-10 RADIANS
C   OF EACH OTHER, BOTH BEARING AND DISTANCE WILL BE SET TO ZERO.
C
C PROGRAM HISTORY LOG:
C   88-08-29  CHASE, P.
C   88-09-23  CHASE, P.   FIX DUMB SOUTH POLE ERROR
C   88-10-05  CHASE, P.   FIX BEARING AMBIGUITY
C   90-04-12  R.E.JONES   CONVERT TO CFT77 FORTRAN
C
C USAGE:    CALL W3FB10(DLAT1, DLON1, DLAT2, DLON2, BEARD, GCDKM)
C
C   INPUT ARGUMENT LIST:
C     DLAT1    - REAL  LATITUDE OF POINT (1) IN DEGREES NORTH.
C     DLON1    - REAL  LONGITUDE OF POINT (1) IN DEGREES EAST.
C     DLAT2    - REAL  LATITUDE OF POINT (2) IN DEGREES NORTH.
C     DLON2    - REAL  LONGITUDE OF POINT (2) IN DEGREES EAST.
C
C   OUTPUT ARGUMENT LIST:
C     BEARD    - REAL  BEARING OF POINT (2) FROM POINT (1) IN
C                COMPASS DEGREES WITH NORTH = 0.0, VALUES FROM
C                -180.0 TO +180.0 DEGREES.
C     GCDKM    - REAL  GREAT CIRCLE DISTANCE FROM POINT (1) TO
C                POINT (2) IN KILOMETERS.
C
C   SUBPROGRAMS CALLED: NONE
C
C REMARKS:  ACCORDING TO THE NMC HANDBOOK, THE EARTH'S RADIUS IS
C   6371.2 KILOMETERS.  THIS IS WHAT WE USE, EVEN THOUGH THE VALUE
C   RECOMMENDED BY THE SMITHSONIAN METEOROLOGICAL HANDBOOK IS
C   6371.221 KM.  (I WOULDN'T WANT YOU TO THINK THAT I DIDN'T KNOW
C   WHAT THE CORRECT VALUE WAS.)
C   METHOD:  THE POLES ARE SPECIAL CASES, AND HANDLED SEPARATELY.
C   OTHERWISE, FROM SPHERICAL TRIGONOMETRY, THE LAW OF COSINES IS USED
C   TO CALCULATE THE THIRD SIDE OF THE SPHERICAL TRIANGLE HAVING
C   SIDES FROM THE POLE TO POINTS (1) AND (2) (THE COLATITUDES).
C   THEN THE LAW OF SINES IS USED TO CALCULATE THE ANGLE AT POINT
C   (1).  A TEST IS APPLIED TO SEE WHETHER THE ARCSINE RESULT MAY BE
C   BE USED AS SUCH, GIVING AN ACUTE ANGLE AS THE BEARING, OR WHETHER
C   THE ARCSINE RESULT SHOULD BE SUBTRACTED FROM PI, GIVING AN OBTUSE
C   ANGLE AS THE BEARING.  THIS TEST IS DERIVED BY CONSTRUCTING A
C   RIGHT SPHERICAL TRIANGLE USING THE POLE, POINT (2), AND THE
C   MERIDIAN THROUGH POINT(1).  THE LATITUDE OF THE RIGHT-ANGLED
C   VERTEX THEN PROVIDES A TEST--IF LATITUDE (1) IS GREATER THAN THIS
C   LATITUDE, THE BEARING ANGLE MUST BE OBTUSE, OTHERWISE ACUTE.
C         IF THE TWO POINTS ARE WITHIN 1E-6 RADIANS OF EACH OTHER
C   A FLAT EARTH IS ASSUMED, AND THE FOUR-QUADRANT ARCTANGENT
C   FUNCTION IS USED TO FIND THE BEARING.  THE Y-DISPLACEMENT IS
C   THE DIFFERENCE IN LATITUDE AND THE X-DISPLACEMENT IS THE
C   DIFFERENCE IN LONGITUDE TIMES COSINE LATITUDE, BOTH IN RADIANS.
C   DISTANCE IS THEN THE DIAGONAL.
C         FUNDAMENTAL TRIGONOMETRIC IDENTITIES ARE USED FREELY, SUCH
C   AS THAT COS(X) = SIN(PI/2 - X), ETC.  SEE ALMOST ANY MATHEMATICAL
C   HANDBOOK, SUCH AS THE C.R.C. STANDARD MATH TABLES UNDER 'RELATIONS
C   IN ANY SPHERICAL TRIANGLE', OR THE NATIONAL BUREAU OF STANDARDS
C   'HANDBOOK OF MATHEMATICAL FUNCTIONS' UNDER SECTION 4.3.149,
C   FORMULAS FOR SOLUTION OF SPHERICAL TRIANGLES.
C         DOUBLE PRECISION IS USED INTERNALLY BECAUSE OF THE WIDE
C   RANGE OF GEOGRAPHIC VALUES THAT MAY BE USED.
C
C ATTRIBUTES:
C   LANGUAGE: CRAY CFT77 FORTRAN
C   MACHINE:  CRAY Y-MP8/832
C
C$$$
C
C *** IMPLICIT TYPE DEFAULTS.....
C
      IMPLICIT REAL  (A-H,O-Z)
C
C *** CONSTANTS......
C
      REAL   PI
      REAL   HALFPI
      REAL   DR
      REAL   RD
      REAL   TDEG, TRAD, TPOD, TFLT
      REAL   EARTHR
      REAL   WHOLCD, HALFCD, QUARCD
C
C *** VARIABLES......
C
      REAL   RLAT1,  RLAT2,  COSLA1, COSLA2, SINLA1, SINLA2
      REAL   DLOND,  RLOND,  COSLO,  SINLO,  SANGG,  ABEAR
      REAL   YDISP,  XDISP,  DDLAT1, DDLAT2, DBANG
      REAL   DLAT1,  DLAT2,  DLON1,  DLON2,  BEARD,  GCDKM
C
C *** CONVERT LATITUDES AND LONGITUDE DIFFERENCE TO RADIANS.
C
      DATA   PI    /3.141592653589793238462643/
      DATA   HALFPI/1.570796326794896619231322/
      DATA   DR    /0.017453292519943295769237/
      DATA   RD    /57.295779513082320876798155/
      DATA   TDEG  /1E-10/, TRAD/1E-10/, TPOD/1E-6/, TFLT/1E-6/
      DATA   EARTHR/6371.2/
      DATA   WHOLCD/360.0/, HALFCD/180.0/, QUARCD/90.0/

      DDLAT1 = DLAT1
      DDLAT2 = DLAT2
      RLAT1  = DR    * DDLAT1
      RLAT2  = DR    * DDLAT2
      DLOND  = DLON2 - DLON1
      IF (DLOND .GT. HALFCD)  DLOND = DLOND - WHOLCD
      IF (DLOND .LT. -HALFCD) DLOND = DLOND + WHOLCD
      RLOND = DR * DLOND
C
C *** FIRST WE ATTACK THE CASES WHERE POINT 1 IS VERY CLOSE TO THE
C *** NORTH OR SOUTH POLES.
C *** HERE WE USE CONVENTIONAL VALUE FOR BEARING.. - LONG (2) AT THE
C *** NORTH POLE, AND + LONG (2) AT THE SOUTH POLE.
C
      IF (ABS(DDLAT1-QUARCD) .LT. TDEG) THEN
        IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
          DBANG = 0.0
          SANGG = 0.0
        ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
          DBANG = 0.0
          SANGG = PI
        ELSE
          DBANG = -DLON2
          SANGG = HALFPI - RLAT2
        ENDIF
      ELSE IF (ABS(DDLAT1+QUARCD) .LT. TDEG) THEN
        IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
          DBANG = 0.0
          SANGG = PI
        ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
          DBANG = 0.0
          SANGG = 0.0
        ELSE
          DBANG = +DLON2
          SANGG = HALFPI + RLAT2
        ENDIF
C
C *** NEXT WE ATTACK THE CASES WHERE POINT 2 IS VERY CLOSE TO THE
C *** NORTH OR SOUTH POLES.
C *** HERE BEARING IS SIMPLY 0 OR 180 DEGREES.
C
      ELSE IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
        DBANG = 0.0
        SANGG = HALFPI - RLAT1
      ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
        DBANG = HALFCD
        SANGG = HALFPI + RLAT1
C
C *** THE CASE REMAINS THAT NEITHER POINT IS AT EITHER POLE.
C *** FIND COSINE AND SINE OF LATITUDES AND LONGITUDE DIFFERENCE
C *** SINCE THEY ARE USED IN MORE THAN ONE FORMULA.
C
      ELSE
        COSLA1 = COS(RLAT1)
        SINLA1 = SIN(RLAT1)
        COSLA2 = COS(RLAT2)
        SINLA2 = SIN(RLAT2)
        COSLO  = COS(RLOND)
        SINLO  = SIN(RLOND)
C
C *** FOLLOWING IS FORMULA FOR GREAT CIRCLE SUBTENDED ANGLE BETWEEN
C *** POINTS IN RADIAN MEASURE.
C
        SANGG = ACOS(SINLA1*SINLA2 + COSLA1*COSLA2*COSLO)
C
C *** IF THE GREAT CIRCLE SUBTENDED ANGLE IS VERY SMALL, FORCE BOTH
C *** BEARING AND DISTANCE TO BE ZERO.
C
        IF (ABS(SANGG) .LT. TRAD) THEN
          DBANG = 0.0
          SANGG = 0.0
C
C *** IF THE GREAT CIRCLE SUBTENDED ANGLE IS JUST SMALL, ASSUME A
C *** FLAT EARTH AND CALCULATE Y- AND X-DISPLACEMENTS.  THEN FIND
C *** BEARING USING THE ARCTANGENT FUNCTION AND DISTANCE USING THE
C *** SQUARE ROOT.
C
        ELSE IF (ABS(SANGG) .LT. TFLT) THEN
          YDISP = RLAT2-RLAT1
          XDISP = RLOND*COSLA2
          ABEAR = ATAN2(XDISP, YDISP)
          DBANG = RD*ABEAR
          SANGG = SQRT(YDISP**2 + XDISP**2)
C
C *** IF THE ANGLE IS RATHER CLOSE TO PI RADIANS, FORCE BEARING TO
C *** BE ZERO AND DISTANCE TO BE PI.
C *** THE TEST FOR 'CLOSE TO PI' IS MORE RELAXED THAN THE TEST FOR
C *** 'CLOSE TO ZERO' TO ALLOW FOR GREATER RELATIVE ERROR.
C
        ELSE IF (ABS(SANGG-PI) .LT. TPOD) THEN
          DBANG = 0.0
          SANGG = PI
C
C *** OTHERWISE COMPUTE THE PRINCIPAL VALUE OF THE BEARING ANGLE
C *** USING THE LAW OF SINES.  THE DIVISION BY THE SINE FORCES US TO
C *** LIMIT THE DOMAIN OF THE ARCSINE TO (-1,1).
C
        ELSE
          ABEAR = ASIN(AMAX1(-1.0,AMIN1(+1.0,COSLA2*SINLO/
     &            SIN(SANGG))))
C
C *** IF THE LONGITUDE DIFFERENCE IS LESS THAN PI/2 IT IS NECESSARY
C *** TO CHECK WHETHER THE BEARING ANGLE IS ACUTE OR OBTUSE BY
C *** COMPARING LATITUDE (1) WITH THE LATITUDE OF THE GREAT CIRCLE
C *** THROUGH POINT (2) NORMAL TO MERIDIAN OF LONGITUDE (1).  IF
C *** LATITUDE (1) IS GREATER, BEARING IS OBTUSE AND THE ACTUAL
C *** BEARING ANGLE IS THE SUPPLEMENT OF THE ANGLE CALCULATED ABOVE.
C
          IF (0.0 .LE. COSLA1*SINLA2 .AND. COSLA1*SINLA2 .LE.
     &    COSLA2*SINLA1*COSLO .OR. COSLA1*SINLA2 .LE. 0.0 .AND.
     &    COSLA2*SINLA1*COSLO .GE. COSLA1*SINLA2) ABEAR =
     &    SIGN(PI,ABEAR) - ABEAR
          DBANG = RD * ABEAR
        ENDIF
      ENDIF
C
C *** THIS FINISHES THE CASE WHERE POINTS ARE NOT AT THE POLES.
C *** NOW CONVERT BEARING TO DEGREES IN RANGE -180 TO +180 AND FIND
C *** GREAT CIRCLE DISTANCE IN KILOMETERS.
C
      IF (DBANG .LE. -HALFCD) DBANG = DBANG + WHOLCD
      IF (DBANG .GT. HALFCD)  DBANG = DBANG - WHOLCD
      GCDKM = EARTHR * SANGG
      BEARD = DBANG
      RETURN
      END