      SUBROUTINE SUNANG(GCLAT,LON,YEAR,DOFYR,MSEC,ZENANG,AZMUTH,DIST)
      IMPLICIT NONE
C
C**********************************************************************
C
C   NAME- SUNANG
C
C   LANGUAGE- FORTRAN 77     TYPE- SUBROUTINE    COMPUTER- CRAY
C
C   VERSION- 1.0    DATE- 12/29/88   PROGRAMMER- KARL W. COX (SMSRC)
C   VERSION- 2.0    DATE- 05/10/96   PROGRAMMER- ALEX POZDNYAK (FDC)
C      MODIFIED TO RUN ON THE CRAY.  REPLACED I*4 WITH INT AND
C      R*4/8 WITH REAL.
C   VERSION- 3.0    DATE- 07/07/97   PROGRAMMER- LEAH WILSON (HSTX)
C      OUTPUTS NODE.
C   VERSION- 3.1    DATE- 06/26/01   PROGRAMMER- FRANK TILLEY (RAYTHEON)
C      REMOVED COMPUTATION OF ORBIT NODE UNDER WR #164.
C   VERSION- 3,2    DATE- 09/04/03   PROGRAMMER- RICKEY ROLLINS (CSC)
C      CHANGED MOD TO AMOD
C  
C   FUNCTION- COMPUTES THE SOLAR ZENITH ANGLE AND THE SOLAR
C             AZIMUTH ANGLE FOR A LOCATION ON THE EARTH
C             FOR YEAR,DAY OF YEAR AND MILLISECONDS AFTER
C             MIDNIGHT IN COORDINATED UNIVERSAL TIME (UTC).
C
C   CALLING ARGUMENTS-
C      VARIABLE                TYPE I/O FUNCTION
C      --------                ---- --- --------
C      GCLAT                   REAL I   LATITUDE (IN DEGREES) OF EARTH
C                                       LOCATION, + IS NORTH, - IS SOUTH
C      LON                     REAL I   LONGITUDE (IN DEGREES) OF EARTH
C                                       LOCATION, + IS EAST, - IS WEST
C      YEAR                    INT  I   YEAR
C      DOFYR                   INT  I   DAY OF YEAR
C      MSEC                    INT  I   MILLISECONDS SINCE MIDNIGHT
C                                       UTC (I.E. 00Z) OF DAY OF YEAR
C      ZENANG                  REAL O   ANGLE OF SUN, MEASURED FROM
C                                       ZENITH (IN DEGREES)
C      AZMUTH                  REAL O   AZIMUTH ANGLE (IN DEGREES) OF
C                                       THE SUN MEASURED FROM THE NORTH
C                                       THROUGH EAST IN THE PLANE OF THE
C                                       HORIZON
C      DIST                    REAL O   DISTANCE FROM THE EARTH TO THE
C                                       SUN (IN AU)
C
C   LOCAL PARAMETERS-
C      VARIABLE                TYPE  FUNCTION
C      --------                ----  --------
C      DELTA                   REAL  CHECK FOR ZERO APPROXIMATION
C      DBL360                  REAL  DOUBLE PRECISION 360
C      M1                      INT   MONTH ONE
C      D1                      INT   DAY ONE
C      MJ2000                  INT   CORRECTION FACTOR FOR THE JULIAN
C                                    DAY CORRESPONDING TO THE YEAR 2000
C      YEARL                   INT   YEAR FOR WHICH THE SLOWLY
C                                    VARYING PARAMETERS WERE CALCULATED
C      DOFYRL                  INT   DAY OF YEAR FOR WHICH THE SLOWLY
C                                    VARYING PARAMETERS WERE CALCULATED
C      HOURL                   REAL  NEAREST HOUR OF THE LAST TIME THE
C                                    SLOWLY VARYING PARAMETERS WERE
C                                    CALCULATED
C      HOURS                   REAL  TIME OF DAY IN HOURS
C      NDAYS                   INT   NUMBER OF WHOLE DAYS FROM 1/1/2000
C      N                       REAL  DAYS + FRACTION FROM 00Z 1/1/2000
C      L                       REAL  MEAN LONGITUDE OF THE SUN
C      G                       REAL  MEAN ANOMALY OF THE SUN
C      LAMDA                   REAL  ECLIPTIC LONGITUDE OF THE SUN
C      EPSLON                  REAL  OBLIQUITY OF THE ECLIPTIC
C      ALPHA                   REAL  ANGLE OF RIGHT ASCENSION
C      TALPHA                  REAL  TANGENT OF ANGLE OF RIGHT ASCENSION
C      RAD                     REAL  DISTANCE FROM EARTH TO SUN
C      SINDEC                  REAL  SINE OF SOLAR DECLINATION
C      COSDEC                  REAL  COSINE OFF SOLAR DECLINATION
C      E                       REAL  CORRECTION FACTOR FOR GMT FROM
C                                    EQUATION OF TIME
C      GHA                     REAL  GMT HOUR ANGLE
C      LHA                     REAL  LOCAL HOUR ANGLE
C      COSLHA                  REAL  COSINE OF HOUR ANGLE
C      SINLAT                  REAL  SINE OF LATITUDE
C      COSLAT                  REAL  COSINE OF LATITUDE
C      COSZA                   REAL  COSINE OF THE SOLAR ZENITH ANGLE
C      SINZA                   REAL  SINE OF THE SOLAR ZENITH ANGLE
C      COSAZ                   REAL  COSINE OF THE SOLAR AZIMUTH ANGLE
C
C
C   SUBPROGRAMS CALLED- NONE
C
C   FUNCTIONS USED - REAL, ABS, MOD, COS, TAN,
C                    ATAN, SIN, ANINT, SIGN, ACOS, SQRT, MJDAY
C
C   COMMON AREAS- NONE
C
C   EXIT STATES- NONE
C
C   REFERENCE
C   ---------
C   "THE ASTRONOMICAL ALMANAC, 1987", EQUATIONS LISTED ON P. C24, EXCEPT
C   TIME-ARGUMENT IS DETERMINED USING THE MODIFIED JULIAN DAY INSTEAD OF
C   THE JULIAN DAY.
C
C   LIMITATIONS
C   -----------
C   EMPLOYS FORMULAE GIVEN IN "THE ASTRONOMICAL ALMANAC" FOR THE SUN'S
C   COORDINATES WITH STATED PRECISION OF 0.01 DEGREE, EQ. OF TIME
C   TO 0.1 MINUTE, AND THE PRECISION OF THE EARTH-SUN DISTANCE IS
C   0.0001 AU.  THE ACCURACY, HOWEVER, IS TYPICALLY 0.001 DEGREES IN
C   POSITION AND 0.05 MINUTES IN TIME.  NO ADJUSTMENTS ARE MADE FOR
C   PRECESSION, NUTATION, OR REFRACTION.
C
C   CAUTIONS
C   -------
C   1)  NO REASONABLENESS CHECKS ARE PERFORMED ON THE INPUT QUANTITIES
C   2)  THIS ROUTINE ASSUMES A SPHERICAL SURFACE FOR DETERMINING THE
C       SOLAR ZENITH AND AZIMUTH ANGLES:  I.E., THAT THE INPUT LATITUDE
C       IS THE GEOCENTRIC LATITUDE.  GEODETIC LATITUDE VALUES, WHICH
C       INCORPORATE THE EFFECTS OF THE EARTH'S OBLATENESS, MUST BE
C       CONVERTED TO GEOCENTRIC LATITUDE VALUES PRIOR TO CALLING SUNANG,
C       OR ERRORS UP TO 10 MINUTES OF ARC IN MID-LATITUDES CAN OCCUR.
C
C**********************************************************************
C
C  SUBROUTINE ARGUMENTS
C
      REAL   GCLAT,LON
      INTEGER   YEAR,DOFYR,MSEC,MJDAY
      REAL   ZENANG,AZMUTH,DIST
C
C  CONSTANTS
C
      INCLUDE'RMPHSCON'
C
      REAL   DELTA
      PARAMETER (DELTA = 1.0E-10)
      REAL   DBL360
      PARAMETER (DBL360 = 360.0)
      INTEGER   M1,D1,MJ2000
      PARAMETER (M1 = 1)
      PARAMETER (D1 = 1)
      PARAMETER (MJ2000 = 51545)
C
C  FOR DETERMINING SOLAR DECLINATION, EQUATION OF TIME AND DISTANCE TO
C  SUN
C
      INTEGER   NDAYS
      REAL   HOURS,N
      REAL   L,G,LAMDA,EPSLON,ALPHA,TALPHA,RAD
C
C  FOR COMPUTING SOLAR ZENITH ANGLE AND AZIMUTH
C
      REAL   GHA,LHA,COSLHA,COSLAT,SINLAT,COSZA,SINZA,COSAZ
C
C  YEAR, DAY AND HOUR FOR WHICH THE SLOWLY CHANGING
C  PARAMETERS HAVE BEEN DETERMINED.  NEW VALUES ARE CALCULATED IF THE
C  INPUT DATA AND TIME DIFFERS BY MORE THAN +/- 0.5 HOURS.
C
      INTEGER   YEARL,DOFYRL
      REAL   HOURL
C
C  SLOWLY VARYING PARAMETERS
C
      REAL   SINDEC,COSDEC,E
C
      SAVE YEARL,DOFYRL,HOURL,E,RAD,SINDEC,COSDEC
      DATA YEARL /-1/, DOFYRL /-1/, HOURL /-1.0/
C
      HOURS = REAL(MSEC)/3600000.0
      IF((YEAR .NE. YEARL) .OR. (DOFYR .NE. DOFYRL) .OR.
     1   (ABS(HOURS - HOURL) .GE. 0.5)) THEN
C
C**** COMPUTE THE SLOWLY VARYING QUANTITIES ONLY ONCE EVERY HOUR
C
C           COMPUTE TIME ARGUMENT
C
          NDAYS = MJDAY(M1,D1,YEAR) + DOFYR - MJ2000
          N = NDAYS + HOURS/24.0 - 0.5
C
C           COMPUTE MEAN LONGITUDE OF THE SUN (IN DEGREES)
C
          L = 280.460 + REAL(AMOD(0.9856474*N,DBL360))
          L = AMOD(L,360.0)
          IF(L .LT. 0.0) L = 360.0 + L
C
C           COMPUTE MEAN ANOMALY (IN DEGREES)
C
          G = 357.528 + REAL(AMOD(0.9856003*N,DBL360))
          IF(G .LT. 0.0) G = 360.0 + G
C
C           COMPUTE ECLIPTIC LONGITUDE OF THE SUN (IN DEGREES)
C
          LAMDA = L + 1.915*SIN(G*D2RAD) + 0.020*SIN(2.0*G*D2RAD)
          LAMDA = AMOD(LAMDA,360.0)
          IF(LAMDA .LT. 0.0) LAMDA = 360.0 + LAMDA
C
C           COMPUTE THE OBLIQUITY OF THE ECLIPTIC (IN DEGREES)
C
          EPSLON = 23.439 - REAL(0.0000004*N)
C
C           COMPUTE THE SINE AND COSINE OF THE SOLAR DECLINATION
C
          SINDEC=SIN(EPSLON*D2RAD)*SIN(LAMDA*D2RAD)
          COSDEC=SQRT(1.0 - SINDEC*SINDEC)
C
C           COMPUTE APPARENT RIGHT ASCENSION OF THE SUN (IN DEGREES)
C
          IF(ABS(COS(LAMDA*D2RAD)) .LT. DELTA) THEN
              ALPHA = LAMDA
            ELSE
              TALPHA = COS(EPSLON*D2RAD)*TAN(LAMDA*D2RAD)
              ALPHA = ATAN(TALPHA)/D2RAD + 180.0*ANINT(LAMDA/180.0)
            ENDIF
C
C           COMPUTE THE EQUATION OF TIME, I.E., THE APPARENT SOLAR TIME
C           MINUS THE MEAN SOLAR TIME (IN HOURS)
C
          E = (L - ALPHA)/15.0
C
C           COMPUTE THE DISTANCE FROM THE EARTH TO THE SUN
C
          RAD = 1.00014 - 0.01671*COS(G*D2RAD)
     1                     - 0.00014*COS(2.0*G*D2RAD)
C
C           SAVE THE YEAR, DAY OF YEAR, AND HOUR FOR WHICH THE SOLAR
C           DECLINATION, EQUATION OF TIME, AND EARTH-SUN DISTANCE HAVE
C           BEEN COMPUTED
C
          YEARL  = YEAR
          DOFYRL = DOFYR
          HOURL  = ANINT(REAL(HOURS))
C
        ENDIF
C
C***   THE FOLLOWING COMPUTATIONS ARE PERFORMED FOR EVERY CALL TO SUNANG
C
C           COMPUTE HOUR ANGLE (IN HOURS) AT GREENWICH, NEGATIVE VALUES
C           INDICATE HOURS BEFORE NOON, POSITIVE VALUES ARE HOURS AFTER
C           NOON
C
      GHA = REAL(HOURS) - 12.0 + E
C
C           COMPUTE THE COSINE OF THE LOCAL HOUR ANGLE
C
      LHA = AMOD(LON + GHA*15.0,360.0)
      IF(LHA .GT.  180.0) LHA = LHA - 360.0
      IF(LHA .LT. -180.0) LHA = LHA + 360.0
      COSLHA = COS(LHA*D2RAD)
C
C           COMPUTE THE SOLAR ZENITH ANGLE
C
      COSLAT=COS(GCLAT*D2RAD)
      SINLAT=SIN(GCLAT*D2RAD)
      COSZA=SINDEC*SINLAT+COSDEC*COSLAT*COSLHA
      IF(ABS(COSZA) .GT. 1.0) COSZA = SIGN(1.0,COSZA)
      ZENANG=ACOS(COSZA)/D2RAD
      SINZA=SQRT(1.-COSZA*COSZA)
C
C           COMPUTE THE SOLAR AZIMUTH ANGLE
C
      IF(SINZA .LT. DELTA) THEN
          AZMUTH = 0.0
        ELSE
          COSAZ = (SINDEC*COSLAT - COSDEC*COSLHA*SINLAT)/SINZA
          IF(ABS(COSAZ) .GT. 1.0) COSAZ = SIGN(1.0,COSAZ)
          AZMUTH = ACOS(COSAZ)/D2RAD
          IF(LHA .GT. 0.0) AZMUTH = 360.0 - AZMUTH
        ENDIF
C
C           THE EARTH TO SUN DISTANCE DOES NOT CHANGE APPRECIABLY IN
C           LESS  THAN AN HOUR
C
      DIST = RAD
C
      RETURN
      END
