SUBROUTINE SPOINTS
     A  (FLAT1, ELON1, FLAT2, ELON2, N, FLAT, ELON)
	 
C	S/R to determine the latitude and longitude of N equally spaced
C	points between the points (FLAT1, ELON1) and (FLAT2, ELON2).
C	All values specified in degrees.  End points are excluded from
C       the set of returned points.

C	Written by G. IRLAM 8-Apr-87
C	Last modified 23-mar-88

C	INPUT DATA  FLAT1   North latitude of first end point
C		    ELON1   East longitude of first end point
C		    FLAT2   North latitude of second end point
C		    ELON2   East longitude of second end point
C		    N	    No. of intermediate points
C	OUTPUT DATA FLAT    N values of latitude of intermediate points
C		    FLON    N values of longitude of intermediate points

	DIMENSION FLAT (N), ELON (N)

	CALL DISTBEAR (FLAT1, ELON1, FLAT2, ELON2, DIST, BEAR)
	D = DIST / (N + 1)

	DO 100 I = 1, N
	  CALL LATLON (FLAT1, ELON1, I * D, BEAR, FLAT (I), ELON (I))
100	continue

	RETURN

        END


	SUBROUTINE DISTBEAR
     A  (X1, Y1, X2, Y2, D, B)
	 
C	S/R to determine the distance and bearing from point
C	(X1, Y1) to point (X2, Y2).  All values specified in degrees.
C	All distances in kilometres.

C	Written by G. IRLAM 8-Apr-87
c       Last modified 15-june-87 david barker

c	formulae courtesy of: "Fundamental Formulas of Physics"
c	  ed: Donald H. Menzel    p. 8, section 2.11
c         and considering the spherical triangle formed by
c         the two given vertices and the North Pole.

C	INPUT DATA  X1	    North latitude of first end point
C		    Y1	    East longitude of first end point
C		    X2	    North latitude of second end point
C		    Y2	    East longitude of second end point
C	OUTPUT DATA D	    Shortest distance between points
C		    B	    Bearing of second point from first point
C			    along along shortest path.  Degrees East
C			    of North , 0 < B <= 360.

	PI = 4.0 * ATAN (1.0)
	RAD = PI / 180
	DEG = 180 / PI
	R0 = 6372.1
c	modified 16-oct-87 to agree with QPAR
	r0 = 6370.

	YY = Y2 - Y1
	IF (YY .LE. -180.0) YY = YY + 360.0
	IF (YY .GT. 180.0) YY = YY - 360.0
	SX1 = SIN (X1 * RAD)
	CX1 = COS (X1 * RAD)
	SX2 = SIN (X2 * RAD)
	CX2 = COS (X2 * RAD)
	CYY = COS (YY * RAD)
	CD = SX1 * SX2 + CX1 * CX2 * CYY
	IF (ABS (CD) .GT. 1.0) CD = SIGN (1.0, CD)
	D = R0 * ACOS (CD)
	SD = SQRT (1.0 - CD * CD)
c	modified 23-mar-88 to account for co-located points, or point 1
c	being N or S pole - in this case, bearing is set to zero
	if(cx1*sd.eq.0.) then
	    b = 0.
	else
	CBN = (SX2 - SX1 * CD) / (CX1 * SD)
	if (abs(cbn) .gt. 1.0) cbn = sign(1.0, cbn)
	bn = acos(cbn) *deg
	if (yy .lt. 0.0) bn = 360. -bn
	b = bn
	endif


	RETURN

	END


	SUBROUTINE LATLON
     A  (X1, Y1, D, B, X2, Y2)
	 
C	S/R to determine the point (X2, Y2) at a distance D and bearing
C	B from the point (X1, Y1).
C	All values specified in degrees.  All distances in kilometres.

C	Written by G. IRLAM 8-Apr-87
c       Last modified 15-june-87 david barker

c	Formulae courtesy of: "Fundamental Formulas of Physics."
c	  ed: Donald H. Menzel     p. 8  section 2.11 
c	  and considering the spherical triangle formed by
c 	  the bearing point, the North Pole, and the side 
c         described by bearing B and distance D.

C	INPUT DATA  X1	    North latitude of bearing point
C		    Y1	    East longitude of bearing point
C		    D	    Shortest distance to new points
C		    B	    Bearing of new point from bearing point
C			    along along shortest path.  Degrees East
C			    of North, 0 < B <= 360.
C	OUTPUT DATA X2	    North latitude of new point
C		    Y2	    East longitude of new point

	PI = 4.0 * ATAN (1.0)
	RAD = PI / 180
	DEG = 180 / PI
	R0 = 6372.1
c	modified 16-oct-87 to agree with QPAR
c	modified 7-jun-88 to have longitudes between 0 and 359.999999
	r0 = 6370.

	BN = B
 	SX1 = SIN (X1 * RAD)
	CX1 = COS (X1 * RAD)
	SD = SIN (D / R0)
	CD = COS (D / R0)
	CBN = COS (BN * RAD)
        sx2 = sx1*cd + cx1*sd*cbn
	IF (ABS (sX2) .GT. 1.0) sX2 = SIGN (1.0, sX2)
	x2 = asin(sx2)*deg

	CX2 = SQRT (1.0 - sX2 * sX2)
	CYY = (CD - SX1 * SX2) / (CX1 * CX2)
	IF (ABS (CYY) .GT. 1.0) CYY = SIGN (1.0, CYY)
	YY = ACOS (CYY) * DEG
	if (bn .gt. 180.0) yy = -yy
	y2 = y1 + yy
	if (y2 .lt. 0.0) y2 = y2 + 360.0

	RETURN

	END