c  
c----------------------------------------------------------------------

      subroutine zensun(iday,time,alat,alon,zenith,azimuth,solfac) 1
c
c   routine:      zensun
c  
c   purpose:  compute the solar zenith and azimuth angles and solar flux
c             multiplier for a given location, time and day.
c
c   input:
c     iday    day of year (used to fraction into a 364.75 day year)
c  
c     time    universal time in decimal hours
c  
c     alat    geographic latitude of point on earth's surface (degrees)
c  
c     alon    geographic longitude of point on earth's surface (degrees)
c  
c   output:
c  
c     zenith  solar zenith angle (degrees)
c  
c     azimuth solar azimuth measured clockwise from due north (degrees)
c
c     solfac  solar flux multiplier. solfac=1./rsun**2 
c             where rsun is the current earth-sun distance in 
c             astronomical units.  
c  
      parameter (pi=3.141592654)
      parameter (dtor=pi/180.)


      dimension nday(74),eqt(74),dec(74)
      data nday/
     &   1,   6,  11,  16,  21,  26,  31,  36,  41,
     &  46,  51,  56,  61,  66,  71,  76,  81,  86,
     &  91,  96, 101, 106, 111, 116, 121, 126, 131,
     & 136, 141, 146, 151, 156, 161, 166, 171, 176,
     & 181, 186, 191, 196, 201, 206, 211, 216, 221,
     & 226, 231, 236, 241, 246, 251, 256, 261, 266,
     & 271, 276, 281, 286, 291, 296, 301, 306, 311,
     & 316, 321, 326, 331, 336, 341, 346, 351, 356,
     & 361, 366/

      data eqt/
     &    -3.23, -5.49, -7.60, -9.48,-11.09,-12.39,-13.34,-13.95,
     &   -14.23,-14.19,-13.85,-13.22,-12.35,-11.26,-10.01, -8.64,
     &    -7.18, -5.67, -4.16, -2.69, -1.29, -0.02,  1.10,  2.05,
     &     2.80,  3.33,  3.63,  3.68,  3.49,  3.09,  2.48,  1.71,
     &     0.79, -0.24, -1.33, -2.41, -3.45, -4.39, -5.20, -5.84,
     &    -6.28, -6.49, -6.44, -6.15, -5.60, -4.82, -3.81, -2.60,
     &    -1.19,  0.36,  2.03,  3.76,  5.54,  7.31,  9.04, 10.69,
     &    12.20, 13.53, 14.65, 15.52, 16.12, 16.41, 16.36, 15.95,
     &    15.19, 14.09, 12.67, 10.93,  8.93,  6.70,  4.32,  1.86,
     &    -0.62, -3.23/

      data dec/
     &   -23.06,-22.57,-21.91,-21.06,-20.05,-18.88,-17.57,-16.13,
     &   -14.57,-12.91,-11.16, -9.34, -7.46, -5.54, -3.59, -1.62,
     &     0.36,  2.33,  4.28,  6.19,  8.06,  9.88, 11.62, 13.29,
     &    14.87, 16.34, 17.70, 18.94, 20.04, 21.00, 21.81, 22.47,
     &    22.95, 23.28, 23.43, 23.40, 23.21, 22.85, 22.32, 21.63,
     &    20.79, 19.80, 18.67, 17.42, 16.05, 14.57, 13.00, 11.33,
     &     9.60,  7.80,  5.95,  4.06,  2.13,  0.19, -1.75, -3.69,
     &    -5.62, -7.51, -9.36,-11.16,-12.88,-14.53,-16.07,-17.50,
     &   -18.81,-19.98,-20.99,-21.85,-22.52,-23.02,-23.33,-23.44,
     &   -23.35,-23.06/

c
c compute the subsolar coordinates
c
      dd=mod(iday-1, 365)+1
      do 10 i=1,74
        if(nday(i) .gt. dd ) goto 20
 10   continue
 20   continue
      frac=(dd-nday(i-1))/(nday(i)-nday(i-1))
      eqtime=eqt(i-1)*(1.-frac)+frac*eqt(i)
      decang=dec(i-1)*(1.-frac)+frac*dec(i)
      sunlat=decang
      sunlon=-15.*(time-12.+eqtime/60.)
      
c      write(12,'(2i4,1p9e11.3)') i,nday(i),frac,eqtime,sunlat,sunlon
c      write(12,'(1p9e11.3)') eqt(i-1),eqt(i),eqtime
c      write(12,'(1p9e11.3)') dec(i-1),dec(i),decang
c
c compute the solar zenith, azimuth and flux multiplier
c
      t0=(90.-alat)*dtor                            
      t1=(90.-sunlat)*dtor                         
      p0=alon*dtor                                  
      p1=sunlon*dtor                               
      zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0) 
c      write(12,'(1p9e11.3)') t0,t1,p0,p1,zz
c      write(12,'(1p9e11.3)') cos(t0),cos(t1),sin(t0),sin(t1),cos(p1-p0)
      xx=sin(t1)*sin(p1-p0)
      yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
      azimuth=atan2(xx,yy)/dtor
      zenith=acos(zz)/dtor                         
c
      rsun=1.-0.01673*cos(.9856*(dd-2.)*dtor)      
      solfac=1./rsun**2                              
      return
      end