subroutine sunmod(radeg,iday,iyr,gmt,i,up,no,ea,sunz,rs)
 
!  Given year, day of year, time in hours (GMT) and latitude and 
!  longitude, returns an accurate solar zenith and azimuth angle.  
!  Based on IAU 1976 Earth ellipsoid.  Method for computing solar 
!  vector and local vertical from Patt and Gregg, 1994, Int. J. 
!  Remote Sensing.  Only outputs solar zenith angle.  This version
!  utilizes a pre-calculation of the local up, north, and east 
!  vectors, since the locations where the solar zenith angle are 
!  calculated in the model are fixed.
 
!  Subroutines required: sun2000
!                        gha2000
!                        jd
 
#include "gloparam.F.h"
      real up(len,3),no(len,3),ea(len,3)
      real suni(3),sung(3)
 
!  Compute sun vector
!   Compute unit sun vector in geocentric inertial coordinates
      sec = gmt*3600.0D0
      call sun2000(radeg,iyr,iday,sec,suni,rs)

!   Get Greenwich mean sidereal angle
      day = float(iday) + sec/86400.0D0
      call gha2000(radeg,iyr,day,gha)
      ghar = gha/radeg

!   Transform Sun vector into geocentric rotating frame
      sung(1) = suni(1)*cos(ghar) + suni(2)*sin(ghar)
      sung(2) = suni(2)*cos(ghar) - suni(1)*sin(ghar)
      sung(3) = suni(3)
 
!  Compute components of spacecraft and sun vector in the
!  vertical (up), North (no), and East (ea) vectors frame
      sunv = 0.0
      sunn = 0.0
      sune = 0.0
      do j = 1,3
       sunv = sunv + sung(j)*up(i,j)
       sunn = sunn + sung(j)*no(i,j)
       sune = sune + sung(j)*ea(i,j)
      enddo
 
!  Compute the solar zenith 
      sunz = radeg*atan2(sqrt(sunn*sunn+sune*sune),sunv)
 
      return
      end
 
! *****************************************************************
      subroutine sun2000(radeg,iyr,iday,sec,sunvec,rs)
 
!  This subroutine computes the Sun vector in geocentric inertial 
!  (equatorial) coodinates.  It uses the model referenced in The 
!  Astronomical Almanac for 1984, Section S (Supplement) and documented
!  in Exact closed-form geolocation algorithm for Earth survey
!  sensors, by F.S. Patt and W.W. Gregg, Int. Journal of Remote
!  Sensing, 1993.  The accuracy of the Sun vector is approximately 0.1 
!  arcminute.
 
!	Arguments:
 
!	Name	Type	I/O	Description
!	--------------------------------------------------------
!	IYR	I*4	 I	Year, four digits (i.e, 1993)
!	IDAY	I*4	 I	Day of year (1-366)
!	SEC	R*8	 I	Seconds of day 
!	SUN(3)	R*8	 O	Unit Sun vector in geocentric inertial 
!				 coordinates of date
!	RS	R*8	 O	Magnitude of the Sun vector (AU)
 
!	Subprograms referenced:
 
!	JD		Computes Julian day from calendar date
!	EPHPARMS	Computes mean solar longitude and anomaly and
!			 mean lunar lontitude and ascending node
!	NUTATE		Compute nutation corrections to lontitude and 
!			 obliquity
 
!	Coded by:  Frederick S. Patt, GSC, November 2, 1992
!	Modified to include Earth constants subroutine by W. Gregg,
!		May 11, 1993.


      real sunvec(3)
      common /nutcm/dpsi,eps,nutime
      parameter (xk=0.0056932)    !Constant of aberration 
      parameter (imon=1)

!   Compute floating point days since Jan 1.5, 2000 
!    Note that the Julian day starts at noon on the specified date
       rjd = float(jd(iyr,imon,iday))
       t = rjd - 2451545.0D0 + (sec-43200.0D0)/86400.0D0

!  Compute solar ephemeris parameters
       call ephparms (t, xls, gs, xlm, omega)

!  Check if need to compute nutation corrections for this day
      nt = int(t)
      if (nt.ne.nutime) then
        nutime = nt
        call nutate (radeg, t, xls, gs, xlm, omega, dpsi, eps)
      end if

!  Compute planet mean anomalies
!   Venus Mean Anomaly 	
      g2 = 50.40828D0 + 1.60213022D0*t
      g2 = mod(g2,360.0)

!   Mars Mean Anomaly 		
      g4 = 19.38816D0 + 0.52402078D0*t
      g4 = mod(g4,360.0)

!  Jupiter Mean Anomaly 
      g5 = 20.35116D0 + 0.08309121D0*t
      g5 = mod(g5,360.0)

!  Compute solar distance (AU)
      rs = 1.00014D0 - 0.01671D0*cos(gs/radeg) 
     *       - 0.00014D0*cos(2.0D0*gs/radeg)

!  Compute Geometric Solar Longitude 
      dls = (6893.0D0 - 4.6543463D-4*t)*sin(gs/radeg) 
     * + 72.0D0*sin(2.0D0*gs/radeg) 
     * - 7.0D0*cos((gs - g5)/radeg)
     * + 6.0D0*sin((xlm - xls)/radeg) 
     * + 5.0D0*sin((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg) 
     * - 5.0D0*cos((2.0D0*gs - 2.0D0*g2)/radeg)
     * - 4.0D0*sin((gs - g2)/radeg) 
     * + 4.0D0*cos((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg) 
     * + 3.0D0*sin((2.0D0*gs - 2.0D0*g2)/radeg)
     * - 3.0D0*sin(g5/radeg) 
     * - 3.0D0*sin((2.0D0*gs - 2.0D0*g5)/radeg)  !arcseconds

      xlsg = xls + dls/3600.0D0

!  Compute Apparent Solar Longitude; includes corrections for nutation 
!   in longitude and velocity aberration
      xlsa = xlsg + dpsi - xk/rs

!   Compute unit Sun vector 
      sunvec(1) = cos(xlsa/radeg)
      sunvec(2) = sin(xlsa/radeg)*cos(eps/radeg)
      sunvec(3) = sin(xlsa/radeg)*sin(eps/radeg)
!	type *,' Sunlon = ',xlsg,xlsa,eps

      return
      end
 
! *********************************************************************
      subroutine gha2000 (radeg, iyr, day, gha)

!  This subroutine computes the Greenwich hour angle in degrees for the
!  input time.  It uses the model referenced in The Astronomical Almanac
!  for 1984, Section S (Supplement) and documented in Exact 
!  closed-form geolocation algorithm for Earth survey sensors, by 
!  F.S. Patt and W.W. Gregg, Int. Journal of Remote Sensing, 1993.
!  It includes the correction to mean sideral time for nutation
!  as well as precession.

!  Calling Arguments

!  Name		Type 	I/O	Description
 
!  iyr		I*4	 I	Year (four digits)
!  day		R*8	 I	Day (time of day as fraction)
!  gha		R*8	 O	Greenwich hour angle (degrees)


!	Subprograms referenced:
 
!	JD		Computes Julian day from calendar date
!	EPHPARMS	Computes mean solar longitude and anomaly and
!			 mean lunar lontitude and ascending node
!	NUTATE		Compute nutation corrections to lontitude and 
!			 obliquity
  	
 
!	Program written by:	Frederick S. Patt
!				General Sciences Corporation
!				November 2, 1992
 
!	Modification History:
 
       common /nutcm/dpsi,eps,nutime
       parameter (imon=1)
       data nutime /-99999/

!  Compute days since J2000
       iday = int(day)
       fday = day - iday
       jday = jd(iyr,imon,iday)
       t = jday - 2451545.5D0 + fday
!  Compute Greenwich Mean Sidereal Time	(degrees)
       gmst = 100.4606184D0 + 0.9856473663D0*t + 2.908D-13*t*t

!  Check if need to compute nutation correction for this day
       nt = int(t)
       if (nt.ne.nutime) then
         nutime = nt
         call ephparms (t, xls, gs, xlm, omega)
         call nutate (radeg, t, xls, gs, xlm, omega, dpsi, eps)
       end if

!  Include apparent time correction and time-of-day
       gha = gmst + dpsi*cos(eps/radeg) + fday*360.0D0
       gha = mod(gha,360.0)
       if (gha.lt.0.0D0) gha = gha + 360.0D0

      return
      end

       subroutine ephparms (t, xls, gs, xlm, omega)

!  This subroutine computes ephemeris parameters used by other Mission
!  Operations routines:  the solar mean longitude and mean anomaly, and
!  the lunar mean longitude and mean ascending node.  It uses the model
!  referenced in The Astronomical Almanac for 1984, Section S 
!  (Supplement) and documented and documented in Exact closed-form 
!  geolocation algorithm for Earth survey sensors, by F.S. Patt and 
!  W.W. Gregg, Int. Journal of Remote Sensing, 1993.  These parameters 
!  are used to compute the solar longitude and the nutation in 
!  longitude and obliquity.

!  Calling Arguments

!  Name		Type 	I/O	Description
 
!  t		R*8	 I	Time in days since January 1, 2000 at 
!				 12 hours UT
!  xls		R*8	 O	Mean solar longitude (degrees)
!  gs		R*8	 O	Mean solar anomaly (degrees)
!  xlm		R*8	 O	Mean lunar longitude (degrees)
!  omega	R*8	 O	Ascending node of mean lunar orbit 
!				 (degrees)
 
 
!	Program written by:	Frederick S. Patt
!				General Sciences Corporation
!				November 2, 1992
 
!	Modification History:
 

!  Sun Mean Longitude 		
      xls = 280.46592D0 + 0.9856473516D0*t
      xls = mod(xls,360.0)
 
!  Sun Mean Anomaly		
      gs = 357.52772D0 + 0.9856002831D0*t 
      gs = mod(gs,360.0)

!  Moon Mean Longitude		
      xlm = 218.31643D0 + 13.17639648D0*t 
      xlm = mod(xlm,360.0)

!  Ascending Node of Moons Mean Orbit 	
      omega = 125.04452D0 - 0.0529537648D0*t 
      omega = mod(omega,360.0)

      return
      end

      subroutine nutate (radeg, t, xls, gs, xlm, omega, dpsi, eps)

!  This subroutine computes the nutation in longitude and the obliquity
!  of the ecliptic corrected for nutation.  It uses the model referenced
!  in The Astronomical Almanac for 1984, Section S (Supplement) and 
!  documented in Exact closed-form geolocation algorithm for Earth 
!  survey sensors, by F.S. Patt and W.W. Gregg, Int. Journal of 
!  Remote Sensing, 1993.  These parameters are used to compute the 
!  apparent time correction to the Greenwich Hour Angle and for the 
!  calculation of the geocentric Sun vector.  The input ephemeris 
!  parameters are computed using subroutine ephparms.  Terms are 
!  included to 0.1 arcsecond.

!  Calling Arguments

!  Name		Type 	I/O	Description
 
!  t		R*8	 I	Time in days since January 1, 2000 at 
!				 12 hours UT
!  xls		R*8	 I	Mean solar longitude (degrees)
!  gs		R*8	 I	Mean solar anomaly   (degrees)
!  xlm		R*8	 I	Mean lunar longitude (degrees)
!  Omega	R*8	 I	Ascending node of mean lunar orbit 
!  				 (degrees)
!  dPsi		R*8	 O	Nutation in longitude (degrees)
!  Eps		R*8	 O	Obliquity of the Ecliptic (degrees)
! 				 (includes nutation in obliquity)
 
 
!	Program written by:	Frederick S. Patt
!				General Sciences Corporation
!				October 21, 1992
 
!	Modification History:
 
!  Nutation in Longitude
      dpsi = - 17.1996D0*sin(omega/radeg) 
     * + 0.2062D0*sin(2.0D0*omega/radeg)
     * - 1.3187D0*sin(2.0D0*xls/radeg) 
     * + 0.1426D0*sin(gs/radeg) 
     * - 0.2274D0*sin(2.0D0*xlm/radeg) 

!  Mean Obliquity of the Ecliptic	
      epsm = 23.439291D0 - 3.560D-7*t 

!  Nutation in Obliquity 
      deps = 9.2025D0*cos(omega/radeg) + 0.5736D0*cos(2.0D0*xls/radeg)

!  True Obliquity of the Ecliptic 
      eps = epsm + deps/3600.0D0

      dpsi = dpsi/3600.0D0

      return
      end
 
! ************************************************************************
      function jd(i,j,k)
 
 
!    This function converts a calendar date to the corresponding Julian
!    day starting at noon on the calendar date.  The algorithm used is
!    from Van Flandern and Pulkkinen, Ap. J. Supplement Series 41, 
!    November 1979, p. 400.
 
 
!	Arguments
      
!     	Name    Type 	I/O 	Description
!     	----	---- 	--- 	-----------
!     	i	I*4  	 I 	Year - e.g. 1970
!     	j       I*4  	 I  	Month - (1-12)
!     	k       I*4  	 I  	Day  - (1-31)
!     	jd      I*4  	 O  	Julian day
 
!     external references
!     -------------------
!      none
 
 
!     Written by Frederick S. Patt, GSC, November 4, 1992
 
 
      jd = 367*i - 7*(i+(j+9)/12)/4 + 275*j/9 + k + 1721014

!  This additional calculation is needed only for dates outside of the 
!   period March 1, 1900 to February 28, 2100
!     	jd = jd + 15 - 3*((i+(j-9)/7)/100+1)/4
        return
        end