C CHANGES TO THE ROUTINES BY A. STOHL
C XI,XI0,ETA,ETA0 ARE DOUBLE PRECISION VARIABLES TO AVOID PROBLEMS
C AT POLES

	SUBROUTINE CC2GLL (STRCMP, XLAT,XLONG, UE,VN, UG,VG)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	DOUBLE PRECISION XPOLG,YPOLG,ALONG,SLONG,CLONG,ROT
	REAL STRCMP(9)
	ALONG = CSPANF( XLONG - STRCMP(2), -180., 180.)
	IF (XLAT.GT.89.985) THEN
C*  NORTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = - STRCMP(1) * ALONG + XLONG - 180.
	ELSEIF (XLAT.LT.-89.985) THEN
C*  SOUTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = - STRCMP(1) * ALONG - XLONG
	ELSE
	  ROT = - STRCMP(1) * ALONG
	ENDIF
	SLONG = SIN( RADPDG * ROT )
	CLONG = COS( RADPDG * ROT )
	XPOLG = SLONG * STRCMP(5) + CLONG * STRCMP(6)
	YPOLG = CLONG * STRCMP(5) - SLONG * STRCMP(6)
	UG = YPOLG * UE + XPOLG * VN
	VG = YPOLG * VN - XPOLG * UE
	RETURN
	END

	SUBROUTINE CCRVLL (STRCMP, XLAT,XLONG, GX,GY)
C*  WRITTEN ON 9/20/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	DOUBLE PRECISION XPOLG,YPOLG,TEMP,ALONG,SLONG,CLONG,CTEMP
	REAL STRCMP(9)
	ALONG = CSPANF( XLONG - STRCMP(2), -180., 180.)
	SLONG = SIN( RADPDG * STRCMP(1) * ALONG)
	CLONG = COS( RADPDG * STRCMP(1) * ALONG)
	XPOLG = - SLONG * STRCMP(5) + CLONG * STRCMP(6)
	YPOLG = CLONG * STRCMP(5) + SLONG * STRCMP(6)
        TEMP = SIN(RADPDG * XLAT)
        CTEMP = COS(RADPDG * XLAT)
        CURV = (STRCMP(1) - TEMP) / CTEMP / REARTH
        GX = CURV * XPOLG
        GY = CURV * YPOLG
	RETURN
	END

	SUBROUTINE CCRVXY (STRCMP, X,Y, GX,GY)
C*  WRITTEN ON 9/20/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	REAL STRCMP(9)
	DOUBLE PRECISION XPOLG,YPOLG,TEMP,YMERC,EFACT,CURV
	TEMP = STRCMP(1) * STRCMP(7) /REARTH
	XPOLG = STRCMP(6) + TEMP * (STRCMP(3) - X)
	YPOLG = STRCMP(5) + TEMP * (STRCMP(4) - Y)
	TEMP = SQRT ( XPOLG ** 2 + YPOLG ** 2 )
	IF (TEMP.GT.0.) THEN
	  YMERC = - LOG( TEMP) /STRCMP(1)
	  EFACT = EXP(YMERC)
 	  CURV = ( (STRCMP(1) - 1.D0) * EFACT +
     A            (STRCMP(1) + 1.D0) / EFACT )
     B           * .5D0 / REARTH
	  GX = XPOLG * CURV / TEMP
	  GY = YPOLG * CURV / TEMP
	ELSE
	  IF (ABS(STRCMP(1)) .EQ. 1.) THEN
	    GX = 0.
	    GY = 0.
	  ELSE
	    GX = 1./REARTH
	    GY = 1./REARTH
	  ENDIF
	ENDIF
	RETURN
	END

	SUBROUTINE CG2CLL (STRCMP, XLAT,XLONG, UG,VG, UE,VN)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	DOUBLE PRECISION XPOLG,YPOLG,ALONG,SLONG,CLONG,ROT
	REAL STRCMP(9)
	ALONG = CSPANF( XLONG - STRCMP(2), -180., 180.)
	IF (XLAT.GT.89.985) THEN
C*  NORTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = - STRCMP(1) * ALONG + XLONG - 180.
	ELSEIF (XLAT.LT.-89.985) THEN
C*  SOUTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = - STRCMP(1) * ALONG - XLONG
	ELSE
	  ROT = - STRCMP(1) * ALONG
	ENDIF
	SLONG = SIN( RADPDG * ROT )
	CLONG = COS( RADPDG * ROT )
	XPOLG = SLONG * STRCMP(5) + CLONG * STRCMP(6)
	YPOLG = CLONG * STRCMP(5) - SLONG * STRCMP(6)
	UE = YPOLG * UG - XPOLG * VG
	VN = YPOLG * VG + XPOLG * UG
	RETURN
	END

	SUBROUTINE CG2CXY (STRCMP, X,Y, UG,VG, UE,VN)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
	REAL STRCMP(9)
	DOUBLE PRECISION XPOLG,YPOLG,TEMP,XI0,ETA0,XI,ETA
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	XI0 = ( X - STRCMP(3) ) * STRCMP(7) / REARTH
	ETA0 = ( Y - STRCMP(4) ) * STRCMP(7) /REARTH
	XI = XI0 * STRCMP(5) - ETA0 * STRCMP(6)
	ETA = ETA0 * STRCMP(5) + XI0 * STRCMP(6)
	RADIAL = 2. * ETA - STRCMP(1) * (XI*XI + ETA*ETA)
	IF (RADIAL.GT.STRCMP(8)) THEN
C*  CASE NORTH OF 89 DEGREES.  METEOROLOGICAL WIND DIRECTION DEFINITION
C*      CHANGES.
	  CALL CNXYLL(STRCMP, XI,ETA, XLAT,XLONG) 
C*  NORTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = STRCMP(1) * (XLONG - STRCMP(2)) - XLONG - 180.
	  SLONG = - SIN( RADPDG * ROT )
	  CLONG = COS( RADPDG * ROT )
	  XPOLG = SLONG * STRCMP(5) + CLONG * STRCMP(6)
	  YPOLG = CLONG * STRCMP(5) - SLONG * STRCMP(6)
	ELSE IF (RADIAL.LT.STRCMP(9)) THEN
C*  CASE SOUTH OF -89 DEGREES.  METEOROLOGICAL WIND DIRECTION DEFINITION
C*      CHANGES.
	  CALL CNXYLL(STRCMP, XI,ETA, XLAT,XLONG) 
C*  SOUTH POLAR METEOROLOGICAL ORIENTATION: "NORTH" ALONG PRIME MERIDIAN
	  ROT = STRCMP(1) * (XLONG - STRCMP(2)) + XLONG
	  SLONG = - SIN( RADPDG * ROT )
	  CLONG = COS( RADPDG * ROT )
	  XPOLG = SLONG * STRCMP(5) + CLONG * STRCMP(6)
	  YPOLG = CLONG * STRCMP(5) - SLONG * STRCMP(6)
	ELSE
C* NORMAL CASE.  METEOROLOGICAL DIRECTION OF WIND RELATED TO TRUE NORTH.
	  XPOLG = STRCMP(6) - STRCMP(1) * XI0
	  YPOLG = STRCMP(5) - STRCMP(1) * ETA0
	  TEMP = SQRT ( XPOLG ** 2 + YPOLG ** 2 )
	  XPOLG = XPOLG / TEMP
	  YPOLG = YPOLG / TEMP
	END IF
	UE = ( YPOLG * UG - XPOLG * VG )
	VN = ( YPOLG * VG + XPOLG * UG )
	RETURN
	END

	REAL FUNCTION CGSZLL (STRCMP, XLAT,XLONG)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180D0/PI)
	REAL STRCMP(9)
	DOUBLE PRECISION SLAT,YMERC,EFACT
	IF (XLAT .GT. 89.985) THEN
C* CLOSE TO NORTH POLE
	  IF (STRCMP(1) .GT. 0.9999) THEN
C* AND TO GAMMA == 1.
	    CGSZLL = 2. * STRCMP(7)
	    RETURN
	  ENDIF
	  EFACT = COS(RADPDG * XLAT)
	  IF (EFACT .LE. 0.) THEN
	    CGSZLL = 0.
	    RETURN
	  ELSE
	    YMERC = - LOG( EFACT /(1. + SIN(RADPDG * XLAT)))
	  ENDIF
	ELSE IF (XLAT .LT. -89.985) THEN
C* CLOSE TO SOUTH POLE
	  IF (STRCMP(1) .LT. -0.9999) THEN
C* AND TO GAMMA == -1.0
	    CGSZLL = 2. * STRCMP(7)
	    RETURN
	  ENDIF
	  EFACT = COS(RADPDG * XLAT)
	  IF (EFACT .LE. 0.) THEN
	    CGSZLL = 0.
	    RETURN
	  ELSE
	    YMERC = LOG( EFACT /(1. - SIN(RADPDG * XLAT)))
	  ENDIF
	ELSE
	SLAT = SIN(RADPDG * XLAT)
	YMERC = LOG((1. + SLAT) / (1. - SLAT))/2.
C	EFACT = EXP(YMERC)
C	CGSZLL = 2. * STRCMP(7) * EXP (STRCMP(1) * YMERC)
C     C			 / (EFACT + 1./EFACT)
	ENDIF
	CGSZLL = STRCMP(7) * COS(RADPDG * XLAT) * EXP(STRCMP(1) *YMERC)
	RETURN
	END

	REAL FUNCTION CGSZXY (STRCMP, X,Y)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2,ALMST1=.9999999)
	REAL STRCMP(9)
	DOUBLE PRECISION YMERC,EFACT
	DOUBLE PRECISION XI0,ETA0,XI,ETA
	XI0 = ( X - STRCMP(3) ) * STRCMP(7) / REARTH
	ETA0 = ( Y - STRCMP(4) ) * STRCMP(7) /REARTH
	XI = XI0 * STRCMP(5) - ETA0 * STRCMP(6)
	ETA = ETA0 * STRCMP(5) + XI0 * STRCMP(6)
	RADIAL = 2. * ETA - STRCMP(1) * (XI*XI + ETA*ETA)
	EFACT = STRCMP(1) * RADIAL
	IF (EFACT .GT. ALMST1) THEN
	  IF (STRCMP(1).GT.ALMST1) THEN
	    CGSZXY = 2. * STRCMP(7)
	  ELSE
	    CGSZXY = 0.
	  ENDIF
	  RETURN
	ENDIF
	IF (ABS(EFACT) .LT. 1.E-2) THEN
	  TEMP = (EFACT / (2. - EFACT) )**2
	  YMERC = RADIAL / (2. - EFACT) * (1.    + TEMP *
     C				          (1./3. + TEMP *
     C				          (1./5. + TEMP *
     C				          (1./7. ))))
	ELSE
	  YMERC = - LOG( 1. - EFACT ) /2. /STRCMP(1)
	ENDIF
	IF (YMERC .GT. 6.) THEN
	  IF (STRCMP(1) .GT. ALMST1) THEN
	    CGSZXY = 2. * STRCMP(7)
	  ELSE
	    CGSZXY = 0.
	  ENDIF
	ELSE IF (YMERC .LT. -6.) THEN
	  IF (STRCMP(1) .LT. -ALMST1) THEN
	    CGSZXY = 2. * STRCMP(7)
	  ELSE
	    CGSZXY = 0.
	  ENDIF
	ELSE
	  EFACT = EXP(YMERC)
	  CGSZXY = 2. * STRCMP(7) * EXP (STRCMP(1) * YMERC)
     C				 / (EFACT + 1./EFACT)
	ENDIF
	RETURN
	END

	SUBROUTINE CLL2XY (STRCMP, XLAT,XLONG, X,Y)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
	REAL STRCMP(9)
	CALL CNLLXY(STRCMP, XLAT,XLONG, XI,ETA)
	X = STRCMP(3) + REARTH/STRCMP(7) *
     C		 (XI * STRCMP(5) + ETA * STRCMP(6) )
	Y = STRCMP(4) + REARTH/STRCMP(7) *
     C		 (ETA * STRCMP(5) - XI * STRCMP(6) )
	RETURN
	END

	SUBROUTINE CNLLXY (STRCMP, XLAT,XLONG, XI,ETA)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
C  MAIN TRANSFORMATION ROUTINE FROM LATITUDE-LONGITUDE TO
C  CANONICAL (EQUATOR-CENTERED, RADIAN UNIT) COORDINATES
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	PARAMETER (ALMST1=.9999999)
	REAL STRCMP(9)
 	DOUBLE PRECISION GAMMA
	DOUBLE PRECISION DLONG,DLAT,SLAT,MERCY,GMERCY
	GAMMA = STRCMP(1)
	DLAT = XLAT
	DLONG = CSPANF(XLONG - STRCMP(2), -180., 180.)
	DLONG = DLONG * RADPDG
	GDLONG = GAMMA * DLONG
	IF (ABS(GDLONG) .LT. .01) THEN
C  CODE FOR GAMMA SMALL OR ZERO.  THIS AVOIDS ROUND-OFF ERROR OR DIVIDE-
C  BY ZERO IN THE CASE OF MERCATOR OR NEAR-MERCATOR PROJECTIONS.
	  GDLONG = GDLONG * GDLONG
	  SNDGAM = DLONG * (1. - 1./6. * GDLONG *
     C			   (1. - 1./20. * GDLONG *
     C			   (1. - 1./42. * GDLONG )))
	  CSDGAM = DLONG * DLONG * .5 *
     C			   (1. - 1./12. * GDLONG *
     C			   (1. - 1./30. * GDLONG * 
     C			   (1. - 1./56. * GDLONG )))
	ELSE
C CODE FOR MODERATE VALUES OF GAMMA
	  SNDGAM = SIN (GDLONG) /GAMMA
	  CSDGAM = (1. - COS(GDLONG) )/GAMMA /GAMMA
	ENDIF
	SLAT = SIN(RADPDG * DLAT)
	IF ((SLAT .GE. ALMST1) .OR. (SLAT .LE. -ALMST1)) THEN
	  ETA = 1./STRCMP(1)
	  XI = 0.
	  RETURN
	ENDIF
	MERCY = .5 * LOG( (1. + SLAT) / (1. - SLAT) )
	GMERCY = GAMMA * MERCY
	IF (ABS(GMERCY) .LT. .001) THEN
C  CODE FOR GAMMA SMALL OR ZERO.  THIS AVOIDS ROUND-OFF ERROR OR DIVIDE-
C  BY ZERO IN THE CASE OF MERCATOR OR NEAR-MERCATOR PROJECTIONS.
	  RHOG1 = MERCY * (1. - .5 * GMERCY *
     C			  (1. - 1./3. * GMERCY *
     C			  (1. - 1./4. * GMERCY ) ) )
	ELSE
C CODE FOR MODERATE VALUES OF GAMMA
	  RHOG1 = (1. - EXP(-GMERCY)) / GAMMA
	ENDIF
	ETA = RHOG1 + (1. - GAMMA * RHOG1) * GAMMA * CSDGAM
	XI = (1. - GAMMA * RHOG1 ) * SNDGAM
	END

	SUBROUTINE CNXYLL (STRCMP, XI,ETA, XLAT,XLONG)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
C  MAIN TRANSFORMATION ROUTINE FROM CANONICAL (EQUATOR-CENTERED,
C  RADIAN UNIT) COORDINATES
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	PARAMETER (ALMST1=.9999999)
	REAL STRCMP(9)
	DOUBLE PRECISION GAMMA,TEMP,ARG1,ARG2,YMERC,ALONG,GXI,CGETA
        DOUBLE PRECISION XI,ETA
	GAMMA = STRCMP(1)
C  CALCULATE EQUIVALENT MERCATOR COORDINATE
	ODIST = XI*XI + ETA*ETA
	ARG2 = 2. * ETA - GAMMA * (XI*XI + ETA*ETA)
	ARG1 = GAMMA * ARG2
C Change by A. Stohl to avoid problems close to the poles
C IF (ARG1 .GE. ALMST1) THEN
C  DISTANCE TO NORTH (OR SOUTH) POLE IS ZERO (OR IMAGINARY ;) )
C XLAT = SIGN(90.,STRCMP(1))
C XLONG = STRCMP(2)
C RETURN
C ENDIF
	IF (ABS(ARG1) .LT. .01) THEN
C  CODE FOR GAMMA SMALL OR ZERO.  THIS AVOIDS ROUND-OFF ERROR OR DIVIDE-
C  BY ZERO IN THE CASE OF MERCATOR OR NEAR-MERCATOR PROJECTIONS.
	  TEMP = (ARG1 / (2. - ARG1) )**2
	  YMERC = ARG2 / (2. - ARG1) * (1.    + TEMP *
     C				       (1./3. + TEMP *
     C				       (1./5. + TEMP *
     C				       (1./7. ))))
	ELSE
C CODE FOR MODERATE VALUES OF GAMMA
	  YMERC = - LOG ( 1. - ARG1 ) /2. / GAMMA
	ENDIF
C  CONVERT YMERC TO LATITUDE
	TEMP = EXP( - ABS(YMERC) )
	XLAT = SIGN(ATAN2((1. - TEMP) * (1. + TEMP), 2. * TEMP), YMERC)
C  FIND LONGITUDES
	GXI = GAMMA*XI
	CGETA = 1. - GAMMA * ETA
	IF ( ABS(GXI) .LT. .01*CGETA ) THEN
C  CODE FOR GAMMA SMALL OR ZERO.  THIS AVOIDS ROUND-OFF ERROR OR DIVIDE-
C  BY ZERO IN THE CASE OF MERCATOR OR NEAR-MERCATOR PROJECTIONS.
	  TEMP = ( GXI /CGETA )**2
	  ALONG = XI / CGETA * (1.    - TEMP *
     C			       (1./3. - TEMP *
     C			       (1./5. - TEMP *
     C			       (1./7.   ))))
	ELSE
C CODE FOR MODERATE VALUES OF GAMMA
	  ALONG = ATAN2( GXI, CGETA) / GAMMA
	ENDIF
	XLONG = SNGL(STRCMP(2) + DGPRAD * ALONG)
	XLAT = XLAT * DGPRAD
	RETURN
	END

	SUBROUTINE CPOLLL (STRCMP, XLAT,XLONG, ENX,ENY,ENZ)
C*  WRITTEN ON 11/23/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180.,DGPRAD=180./PI)
	DOUBLE PRECISION XPOLG,YPOLG,ALONG,SLONG,CLONG,ROT
	REAL STRCMP(9)
	ALONG = CSPANF( XLONG - STRCMP(2), -180., 180.)
	  ROT = - STRCMP(1) * ALONG
	SLONG = SIN( RADPDG * ROT )
	CLONG = COS( RADPDG * ROT )
	XPOLG = SLONG * STRCMP(5) + CLONG * STRCMP(6)
	YPOLG = CLONG * STRCMP(5) - SLONG * STRCMP(6)
	CLAT = COS(RADPDG * XLAT)
	ENX = CLAT * XPOLG
	ENY = CLAT * YPOLG
	ENZ = SIN(RADPDG * XLAT)
	RETURN
	END

	SUBROUTINE CPOLXY (STRCMP, X,Y, ENX,ENY,ENZ)
C*  WRITTEN ON 11/26/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	REAL STRCMP(9)
	DOUBLE PRECISION XPOL,YPOL,TEMP,XI0,ETA0,XI,ETA,RADIAL
	DOUBLE PRECISION TEMP2,YMERC,ARG,OARG,CLAT
	XI0 = ( X - STRCMP(3) ) * STRCMP(7) / REARTH
	ETA0 = ( Y - STRCMP(4) ) * STRCMP(7) /REARTH
	XI = XI0 * STRCMP(5) - ETA0 * STRCMP(6)
	ETA = ETA0 * STRCMP(5) + XI0 * STRCMP(6)
	RADIAL = 2. * ETA -  STRCMP(1) * (XI*XI + ETA*ETA)
	TEMP = STRCMP(1) * RADIAL
	IF (TEMP .GE. 1.) THEN
	  ENX = 0.
	  ENY = 0.
	  ENZ = SIGN(1.,STRCMP(1))
	  RETURN
	ENDIF
	IF (ABS(TEMP).LT.1.E-2) THEN
	  TEMP2 = (TEMP / (2. - TEMP))**2
	  YMERC = RADIAL / (2. - TEMP) * (1. + TEMP2 *
     C					 (1./3. + TEMP2 *
     C					 (1./5. + TEMP2 *
     C					 (1./7.))))
	ELSE
	  YMERC = -.5 * LOG(1. - TEMP) / STRCMP(1)
	ENDIF
	ARG = EXP( YMERC )
	OARG = 1./ARG
	CLAT = 2./(ARG + OARG)
	ENZ = (ARG - OARG) * CLAT /2.
	TEMP = CLAT / SQRT(1. - TEMP)
	XPOL = - XI * STRCMP(1) * TEMP
	YPOL = (1. - ETA * STRCMP(1) ) * TEMP
	ENX = XPOL * STRCMP(5) + YPOL * STRCMP(6)
	ENY = YPOL * STRCMP(5) - XPOL * STRCMP(6)
	RETURN
	END

	REAL FUNCTION CSPANF (VALUE, BEGIN, END)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
C* REAL FUNCTION CSPANF RETURNS A VALUE IN THE INTERVAL (BEGIN,END]
C* WHICH IS EQUIVALENT TO VALUE, MOD (END - BEGIN).  IT IS USED TO
C* REDUCE PERIODIC VARIABLES TO A STANDARD RANGE.  IT ADJUSTS FOR THE
C* BEHAVIOR OF THE MOD FUNCTION WHICH PROVIDES POSITIVE RESULTS FOR
c* POSITIVE INPUT, AND NEGATIVE RESULTS FOR NEGATIVE INPUT
C* INPUT:
C*       VALUE - REAL NUMBER TO BE REDUCED TO THE SPAN
C*       BEGIN - FIRST VALUE OF THE SPAN
C*       END   - LAST VALUE OF THE SPAN
C* RETURNS:
C*       THE REDUCED VALUE
C* EXAMPLES:
C*      ALONG = CSPANF(XLONG, -180., +180.)
C*      DIR  = CSPANF(ANGLE, 0., 360.)
	REAL FIRST,LAST
	FIRST = MIN(BEGIN,END)
	LAST = MAX(BEGIN,END)
	VAL = MOD( VALUE - FIRST , LAST - FIRST)
	IF ( VAL . LE. 0.) THEN
	  CSPANF = VAL + LAST
	ELSE
	  CSPANF = VAL + FIRST
	ENDIF
	RETURN
	END

	SUBROUTINE CXY2LL (STRCMP, X,Y, XLAT,XLONG)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (REARTH=6 371.2)
        DOUBLE PRECISION XI0,ETA0,XI,ETA
	REAL STRCMP(9)
	XI0 = ( X - STRCMP(3) ) * STRCMP(7) / REARTH
	ETA0 = ( Y - STRCMP(4) ) * STRCMP(7) /REARTH
	XI = XI0 * STRCMP(5) - ETA0 * STRCMP(6)
	ETA = ETA0 * STRCMP(5) + XI0 * STRCMP(6)
	CALL CNXYLL(STRCMP, XI,ETA, XLAT,XLONG)
	XLONG = CSPANF(XLONG, -180., 180.)
	RETURN
	END

	REAL FUNCTION EQVLAT (XLAT1,XLAT2)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	SIND(X) = SIN (RADPDG*X)
	SINL1 = SIND (XLAT1)
	SINL2 = SIND (XLAT2)
	IF (ABS(SINL1 - SINL2) .GT. .001) THEN
	  AL1 = LOG((1. - SINL1)/(1. - SINL2))
	  AL2 = LOG((1. + SINL1)/(1. + SINL2))
	ELSE
C  CASE LAT1 NEAR OR EQUAL TO LAT2
	  TAU = - (SINL1 - SINL2)/(2. - SINL1 - SINL2)
	  TAU = TAU*TAU
	  AL1  = 2. / (2. - SINL1 - SINL2) * (1.    + TAU *
     C 					     (1./3. + TAU *
     C					     (1./5. + TAU *
     C					     (1./7.))))
	  TAU =   (SINL1 - SINL2)/(2. + SINL1 + SINL2)
	  TAU = TAU*TAU
	  AL2  = -2. / (2. + SINL1 + SINL2) * (1.    + TAU *
     C 					      (1./3. + TAU *
     C					      (1./5. + TAU *
     C					      (1./7.))))
	ENDIF
	EQVLAT = ASIN((AL1 + AL2) / (AL1 - AL2))/RADPDG
	RETURN
	END

	SUBROUTINE STCM1P(STRCMP, X1,Y1, XLAT1,XLONG1,
     C XLATG,XLONGG, GRIDSZ, ORIENT)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	REAL STRCMP(9)
	DO K=3,4
	  STRCMP (K) = 0.
	ENDDO
        TURN = RADPDG * (ORIENT - STRCMP(1) *
     C            CSPANF(XLONGG - STRCMP(2), -180., 180.) )
	STRCMP (5) = COS (TURN)
	STRCMP (6) = - SIN (TURN)
	STRCMP (7) = 1.
	STRCMP (7) = GRIDSZ * STRCMP(7)
     C             / CGSZLL(STRCMP, XLATG, STRCMP(2))
	CALL CLL2XY (STRCMP, XLAT1,XLONG1, X1A,Y1A)
	STRCMP(3) = STRCMP(3) + X1 - X1A
	STRCMP(4) = STRCMP(4) + Y1 - Y1A
	RETURN
	END

	SUBROUTINE STCM2P(STRCMP, X1,Y1, XLAT1,XLONG1,
     C X2,Y2, XLAT2,XLONG2)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	REAL STRCMP(9)
	DO K=3,6
	  STRCMP (K) = 0.
	ENDDO
	STRCMP (5) = 1.
	STRCMP (7) = 1.
	CALL CLL2XY (STRCMP, XLAT1,XLONG1, X1A,Y1A)
	CALL CLL2XY (STRCMP, XLAT2,XLONG2, X2A,Y2A)
	DEN = SQRT( (X1 - X2)**2 + (Y1 - Y2)**2 )
	DENA = SQRT( (X1A - X2A)**2 + (Y1A - Y2A)**2 )
	STRCMP(5) = ((X1A - X2A)*(X1 - X2) + (Y1A - Y2A) * (Y1 - Y2))
     C  /DEN /DENA
	STRCMP(6) = ((Y1A - Y2A)*(X1 - X2) - (X1A - X2A) * (Y1 - Y2))
     C  /DEN /DENA
	STRCMP (7) = STRCMP(7) * DENA / DEN
	CALL CLL2XY (STRCMP, XLAT1,XLONG1, X1A,Y1A)
	STRCMP(3) = STRCMP(3) + X1 - X1A
	STRCMP(4) = STRCMP(4) + Y1 - Y1A
	RETURN
	END

C*  GENERAL CONFORMAL MAP ROUTINES FOR METEOROLOGICAL MODELERS
C*  WRITTEN ON 3/31/94 BY

C* Dr. Albion Taylor
C* NOAA / OAR / ARL                  Phone: (301) 713-0295 x 132
C* Rm. 3151, 1315 East-West Highway  Fax:   (301) 713-0119
C* Silver Spring, MD 20910           E-mail: ADTaylor@arlrisc.ssmc.noaa.gov 

C*  SUBROUTINE STLMBR (STRCMP, TNGLAT, CLONG)
C*    THIS ROUTINE INITIALIZES THE MAP STRUCTURE ARRAY STRCMP TO
C*    THE FORM OF A SPECIFIC MAP PROJECTION
C*  INPUTS:
C*    TNGLAT - THE LATITUDE AT WHICH THE PROJECTION WILL BE TANGENT
C*             TO THE EARTH.  +90. FOR NORTH POLAR STEREOGRAPHIC,
C*             -90. FOR SOUTH POLAR STEREOGRAPHIC, 0. FOR MERCATOR,
C*             AND OTHER VALUES FOR LAMBERT CONFORMAL. 
C*             -90 <= TNGLAT <= 90.
C*    CLONG -  A LONGITUDE IN THE REGION UNDER CONSIDERATION.  LONGITUDES
C*             BETWEEN CLONG-180. AND CLONG+180.  WILL BE MAPPED IN ONE
C*             CONNECTED REGION
C*  OUTPUTS:
C*    STRCMP - A 9-VALUE MAP STRUCTURE ARRAY FOR USE WITH SUBSEQUENT
C*             CALLS TO THE COORDINATE TRANSFORM ROUTINES.
C*
C*  REAL FUNCTION EQVLAT (XLAT1,XLAT2)
C*    THIS FUNCTION IS PROVIDED TO ASSIST IN FINDING THE TANGENT LATITUDE
C*    EQUIVALENT TO THE 2-REFERENCE LATITUDE SPECIFICATION IN THE LEGEND
C*    OF MOST LAMBERT CONFORMAL MAPS.  IF THE MAP SPECIFIES "SCALE 
C*    1:XXXXX TRUE AT 40N AND 60N", THEN EQVLAT(40.,60.) WILL RETURN THE
C*    EQUIVALENT TANGENT LATITUDE.
C*  INPUTS:
C*    XLAT1,XLAT2:  THE TWO LATITUDES SPECIFIED IN THE MAP LEGEND
C*  RETURNS:
C*    THE EQUIVALENT TANGENT LATITUDE
C*  EXAMPLE:  CALL STLMBR(STRCMP, EQVLAT(40.,60.), 90.)

C*  SUBROUTINE STCM2P (STRCMP, X1,Y1, XLAT1,XLONG1,
C*          X2,Y2, XLAT2,XLONG2)
C*  SUBROUTINE STCM1P (STRCMP, X1,Y1, XLAT1,XLONG1,
C*          XLATG,XLONGG, GRIDSZ, ORIENT)
C*    THESE ROUTINES COMPLETE THE SPECIFICATION OF THE MAP STRUCTURE
C*    ARRAY BY CONFORMING THE MAP COORDINATES TO THE SPECIFICATIONS
C*    OF A PARTICULAR GRID.  EITHER STCM1P OR STCM2P MUST BE CALLED,
C*    BUT NOT BOTH
C*  INPUTS:
C*    STRCMP - A 9-VALUE MAP STRUCTURE ARRAY, SET TO A PARTICULAR MAP
C*             FORM BY A PREVIOUS CALL TO STLMBR
C*    FOR STCM2P:
C*      X1,Y1, X2,Y2 - THE MAP COORDINATES OF TWO POINTS ON THE GRID
C*      XLAT1,XLONG1, XLAT2,XLONG2 - THE GEOGRAPHIC COORDINATES OF THE
C*             SAME TWO POINTS
C*    FOR STCM1P:
C*      X1,Y1 - THE MAP COORDINATES OF ONE POINT ON THE GRID
C*      XLAT1,XLONG1 - THE GEOGRAPHIC COORDINATES OF THE SAME POINT
C*      XLATG,XLONGG - LATITUDE AND LONGITUDE OF REFERENCE POINT FOR
C*             GRIDSZ AND ORIENTATION SPECIFICATION.
C*      GRIDSZ - THE DESIRED GRID SIZE IN KILOMETERS, AT XLATG,XLONGG
C*      ORIENT - THE ANGLE, WITH RESPECT TO NORTH, OF A Y-GRID LINE, AT
C*             THE POINT XLATG,XLONGG
C*  OUTPUTS:
C*    STRCMP - A 9-VALUE MAP STRUCTURE ARRAY, FULLY SET FOR USE BY
C*             OTHER SUBROUTINES IN THIS SYSTEM

C*  SUBROUTINE CLL2XY (STRCMP, XLAT,XLONG, X,Y)
C*  SUBROUTINE CXY2LL (STRCMP, X,Y, XLAT,XLONG)
C*     THESE ROUTINES CONVERT BETWEEN MAP COORDINATES X,Y
C*     AND GEOGRAPHIC COORDINATES XLAT,XLONG
C*  INPUTS:
C*     STRCMP(9) - 9-VALUE MAP STRUCTURE ARRAY
C*     FOR CLL2XY:  XLAT,XLONG - GEOGRAPHIC COORDINATES
C*     FOR CXY2LL:  X,Y - MAP COORDINATES
C*  OUTPUTS:
C*     FOR CLL2XY:  X,Y - MAP COORDINATES
C*     FOR CXY2LL:  XLAT,XLONG - GEOGRAPHIC COORDINATES

C*  SUBROUTINE CC2GXY (STRCMP, X,Y, UE,VN, UG,VG)
C*  SUBROUTINE CG2CXY (STRCMP, X,Y, UG,VG, UE,VN)
C*  SUBROUTINE CC2GLL (STRCMP, XLAT,XLONG, UE,VN, UG,VG)
C*  SUBROUTINE CG2CLL (STRCMP, XLAT,XLONG, UG,VG, UE,VN)
C*     THESE SUBROUTINES CONVERT VECTOR WIND COMPONENTS FROM
C*     GEOGRAPHIC, OR COMPASS, COORDINATES, TO MAP OR
C*     GRID COORDINATES.  THE SITE OF THE WIND TO BE
C*     CONVERTED MAY BE GIVEN EITHER IN GEOGRAPHIC OR
C*     MAP COORDINATES.  WIND COMPONENTS ARE ALL IN KILOMETERS
C*     PER HOUR, WHETHER GEOGRAPHIC OR MAP COORDINATES.
C*  INPUTS:
C*    STRCMP(9) - 9-VALUE MAP STRUCTURE ARRAY
C*    FOR CC2GXY AND CG2CXY:  X,Y        -  MAP COORDINATES OF SITE
C*    FOR CC2GLL AND CG2CLL:  XLAT,XLONG -  GEOGRAPHIC COORDINATES OF SITE
C*    FOR CC2GXY AND CC2GLL:  UE,VN - EAST AND NORTH WIND COMPONENTS
C*    FOR CG2CXY AND CG2CLL:  UG,VG - X- AND Y- DIRECTION WIND COMPONENTS
C*  OUTPUTS:
C*    FOR CC2GXY AND CC2GLL:  UG,VG - X- AND Y- DIRECTION WIND COMPONENTS
C*    FOR CG2CXY AND CG2CLL:  UE,VN - EAST AND NORTH WIND COMPONENTS

C*  SUBROUTINE CCRVXY (STRCMP, X, Y,       GX,GY)
C*  SUBROUTINE CCRVLL (STRCMP, XLAT,XLONG, GX,GY)
C*    THESE SUBROUTINES RETURN THE CURVATURE VECTOR (GX,GY), AS REFERENCED
C*    TO MAP COORDINATES, INDUCED BY THE MAP TRANSFORMATION.  WHEN
C*    NON-LINEAR TERMS IN WIND SPEED ARE IMPORTANT, A "GEODESIC" FORCE
C*    SHOULD BE INCLUDED IN THE VECTOR FORM [ (U,U) G - (U,G) U ] WHERE THE
C*    INNER PRODUCT (U,G) IS DEFINED AS UX*GX + UY*GY.
C*  INPUTS:
C*    STRCMP(9) - 9-VALUE MAP STRUCTURE ARRAY
C*    FOR CCRVXY:  X,Y        -  MAP COORDINATES OF SITE
C*    FOR CCRVLL:  XLAT,XLONG -  GEOGRAPHIC COORDINATES OF SITE
C*  OUTPUTS:
C*    GX,GY       - VECTOR COEFFICIENTS OF CURVATURE, IN UNITS RADIANS
C*                  PER KILOMETER

C*  REAL FUNCTION CGSZLL (STRCMP, XLAT,XLONG)
C*  REAL FUNCTION CGSZXY (STRCMP, X,Y)
C*    THESE FUNCTIONS RETURN THE SIZE, IN KILOMETERS, OF EACH UNIT OF
C*    MOTION IN MAP COORDINATES (GRID SIZE).  THE GRID SIZE AT ANY
C*    LOCATION DEPENDS ON THAT LOCATION; THE POSITION MAY BE GIVEN IN
C*    EITHER MAP OR GEOGRAPHIC COORDINATES.
C*  INPUTS:
C*    STRCMP(9) - 9-VALUE MAP STRUCTURE ARRAY
C*    FOR CGSZXY:  X,Y        -  MAP COORDINATES OF SITE
C*    FOR CGSZLL:  XLAT,XLONG -  GEOGRAPHIC COORDINATES OF SITE
C*  RETURNS:
C*    GRIDSIZE IN KILOMETERS AT GIVEN SITE.

C*  SUBROUTINE CPOLXY (STRCMP, X,Y, ENX,ENY,ENZ)
C*  SUBROUTINE CPOLLL (STRCMP, XLAT,XLONG, ENX,ENY,ENZ)
C*    THESE SUBROUTINES PROVIDE 3-D VECTOR COMPONENTS OF A UNIT VECTOR
C*    IN THE DIRECTION OF THE NORTH POLAR AXIS.  WHEN MULTIPLIED
C*    BY TWICE THE ROTATION RATE OF THE EARTH (2 * PI/24 HR), THE
C*    VERTICAL COMPONENT YIELDS THE CORIOLIS FACTOR.
C*  INPUTS:
C*    STRCMP(9) - 9-VALUE MAP STRUCTURE ARRAY
C*    FOR CPOLXY:  X,Y        -  MAP COORDINATES OF SITE
C*    FOR CPOLLL:  XLAT,XLONG -  GEOGRAPHIC COORDINATES OF SITE
C*  RETURNS:
C*    ENX,ENY,ENZ THE DIRECTION COSINES OF A UNIT VECTOR IN THE
C*    DIRECTION OF THE ROTATION AXIS OF THE EARTH

C*  SUBROUTINE CNLLXY (STRCMP, XLAT,XLONG, XI,ETA)
C*  SUBROUTINE CNXYLL (STRCMP, XI,ETA, XLAT,XLONG)
C*    THESE SUBROUTINES PERFORM THE UNDERLYING TRANSFORMATIONS FROM
C*    GEOGRAPHIC COORDINATES TO AND FROM CANONICAL (EQUATOR CENTERED)
C*    COORDINATES.  THEY ARE CALLED BY CXY2LL AND CLL2XY, BUT ARE NOT
C*    INTENDED TO BE CALLED DIRECTLY

C*  REAL FUNCTION CSPANF (VALUE, BEGIN, END)
C*    THIS FUNCTION ASSISTS OTHER ROUTINES IN PROVIDING A LONGITUDE IN
C*    THE PROPER RANGE.  IT ADDS TO VALUE WHATEVER MULTIPLE OF 
C*    (END - BEGIN) IS NEEDED TO RETURN A NUMBER BEGIN < CSPANF <= END

	SUBROUTINE STLMBR(STRCMP, TNGLAT, XLONG)
C*  WRITTEN ON 3/31/94 BY Dr. Albion Taylor  NOAA / OAR / ARL
	PARAMETER (PI=3.14159265358979,RADPDG=PI/180,DGPRAD=180/PI)
	PARAMETER (REARTH=6 371.2)
	REAL STRCMP(9)
	STRCMP(1) = SIN(RADPDG * TNGLAT)
C*  GAMMA = SINE OF THE TANGENT LATITUDE
	STRCMP(2) = CSPANF( XLONG, -180., +180.)
C* LAMBDA_0 = REFERENCE LONGITUDE
	STRCMP(3) = 0.
C* X_0 = X- GRID COORDINATE OF ORIGIN (XI,ETA) = (0.,0.)
	STRCMP(4) = 0.
C* y_0 = Y-GRID COORDINATE OF ORIGIN (XI,ETA) = (0.,0.)
	STRCMP(5) = 1.
C* COSINE OF ROTATION ANGLE FROM XI,ETA TO X,Y
	STRCMP(6) = 0.
C* SINE OF ROTATION ANGLE FROM XI,ETA TO X,Y
	STRCMP(7) = REARTH
C* GRIDSIZE IN KILOMETERS AT THE EQUATOR
	CALL CNLLXY(STRCMP, 89.,XLONG, XI,ETA)
	STRCMP(8) = 2. * ETA - STRCMP(1) * ETA * ETA
C* RADIAL COORDINATE FOR 1 DEGREE FROM NORTH POLE
	CALL CNLLXY(STRCMP, -89.,XLONG, XI,ETA)
   	STRCMP(9) = 2. * ETA - STRCMP(1) * ETA * ETA
C* RADIAL COORDINATE FOR 1 DEGREE FROM SOUTH POLE
	RETURN
	END