```
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
```