subroutine zensun(day,time,lat,lon,pi,sun_zenith,sun_azimuth)
!$$$ subprogram documentation block
! . . . .
! subprogram: zensun make sun zenith and sun azimuth angle
!
! prgmmr: Paul Ricchiazzi org: Earth Space Research Group,UCSB date: 1992-10-23
!
! abstract:
! Compute solar position information as a function of
! geographic coordinates, date and time.
!
! program history log:
! 2005-10-21 kazumori - reformatted for GSI
!
! input argument list:
! day - Julian day (positive scalar or vector)
! (spring equinox = 80)
! (summer solstice= 171)
! (fall equinox = 266)
! (winter solstice= 356)
! time - Universal Time in hours (scalar or vector)
! lat - geographic latitude of point on earth's surface (degrees)
! lon - geographic longitude of point on earth's surface (degrees)
!
! output argument list:
! sun_zenith - solar zenith angle
! sun_azimuth - solar azimuth angle
!
! comments:
!
!
! PROCEDURE:
!
! 1. Calculate the subsolar point latitude and longitude, based on
! DAY and TIME. Since each year is 365.25 days long the exact
! value of the declination angle changes from year to year. For
! precise values consult THE AMERICAN EPHEMERIS AND NAUTICAL
! ALMANAC published yearly by the U.S. govt. printing office. The
! subsolar coordinates used in this code were provided by a
! program written by Jeff Dozier.
!
! 2. Given the subsolar latitude and longitude, spherical geometry is
! used to find the solar zenith, azimuth and flux multiplier.
!
! eqt = equation of time (minutes) ! solar longitude correction = -15*eqt
! dec = declination angle (degrees) = solar latitude
!
! LOWTRAN v7 data (25 points)
! The LOWTRAN solar position data is characterized by only 25 points.
! This should predict the subsolar angles within one degree. For
! increased accuracy add more data points.
!
!nday=[ 1., 9., 21., 32., 44., 60., 91., 121., 141., 152.,$
! 160., 172., 182., 190., 202., 213., 244., 274., 305., 309.,$
! 325., 335., 343., 355., 366.]
!
!eqt=[ -3.23, -6.83,-11.17,-13.57,-14.33,-12.63, -4.2, 2.83, 3.57, 2.45,$
! 1.10, -1.42, -3.52, -4.93, -6.25, -6.28,-0.25, 10.02, 16.35, 16.38,$
! 14.3, 11.27, 8.02, 2.32, -3.23]
!
!dec=[-23.07,-22.22,-20.08,-17.32,-13.62, -7.88, 4.23, 14.83, 20.03, 21.95,$
! 22.87, 23.45, 23.17, 22.47, 20.63, 18.23, 8.58, -2.88,-14.18,-15.45,$
! -19.75,-21.68,-22.75,-23.43,-23.07]
!
! Analemma information from Jeff Dozier
! This data is characterized by 74 points
!
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
!
use kinds, only: r_kind,i_kind
implicit none
integer(i_kind),intent(in):: day
real(r_kind),intent(in) :: pi,time,lat,lon
integer(i_kind) di
real(r_kind) z,gaz
real(r_kind) ut,noon
real(r_kind) y(5),y2(5),x(2,5),Tx(5,2),xTx(2,2),aTx(5,2),det
real(r_kind) tt,eqtime,decang,latsun,lonsun
real(r_kind) nday(74),eqt(74),dec(74)
real(r_kind) beta(2), beta2(2), a(2,2)
real(r_kind) t0,t1,p0,p1,zz,xx,yy
real(r_kind),intent(out) :: sun_zenith,sun_azimuth
data nday/1.0, 6.0, 11.0, 16.0, 21.0, 26.0, 31.0, 36.0, 41.0, 46.0,&
51.0, 56.0, 61.0, 66.0, 71.0, 76.0, 81.0, 86.0, 91.0, 96.0,&
101.0, 106.0, 111.0, 116.0, 121.0, 126.0, 131.0, 136.0, 141.0, 146.0,&
151.0, 156.0, 161.0, 166.0, 171.0, 176.0, 181.0, 186.0, 191.0, 196.0,&
201.0, 206.0, 211.0, 216.0, 221.0, 226.0, 231.0, 236.0, 241.0, 246.0,&
251.0, 256.0, 261.0, 266.0, 271.0, 276.0, 281.0, 286.0, 291.0, 296.0,&
301.0, 306.0, 311.0, 316.0, 321.0, 326.0, 331.0, 336.0, 341.0, 346.0,&
351.0, 356.0, 361.0, 366.0/
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/
!
! compute the subsolar coordinates
!
tt= mod(real((int(day)+time/24.-1.)),365.25) +1. ! fractional day number
! with 12am 1jan = 1.
do di = 1, 73
if ((tt .ge. nday(di)) .and. (tt .le. nday(di+1))) exit
end do
!============== Perform a least squares regression on doy**3 ==============
x(1,:) = 1.0
if ((di .ge. 3) .and. (di .le. 72)) then
y(:) = eqt(di-2:di+2)
y2(:) = dec(di-2:di+2)
x(2,:) = nday(di-2:di+2)**3
end if
if (di .eq. 2) then
y(1) = eqt(73)
y(2:5) = eqt(di-1:di+2)
y2(1) = dec(73)
y2(2:5) = dec(di-1:di+2)
x(2,1) = nday(73)**3
x(2,2:5) = (365.+nday(di-1:di+2))**3
end if
if (di .eq. 1) then
y(1:2) = eqt(72:73)
y(3:5) = eqt(di:di+2)
y2(1:2) = dec(72:73)
y2(3:5) = dec(di:di+2)
x(2,1:2) = nday(72:73)**3
x(2,3:5) = (365.+nday(di:di+2))**3
end if
if (di .eq. 73) then
y(1:4) = eqt(di-2:di+1)
y(5) = eqt(2)
y2(1:4) = dec(di-2:di+1)
y2(5) = dec(2)
x(2,1:4) = nday(di-2:di+1)**3
x(2,5) = (365.+nday(2))**3
end if
if (di .eq. 74) then
y(1:3) = eqt(di-2:di)
y(4:5) = eqt(2:3)
y2(1:3) = dec(di-2:di)
y2(4:5) = dec(2:3)
x(2,1:3) = nday(di-2:di)**3
x(2,4:5) = (365.+nday(2:3))**3
end if
! Tx = transpose(x)
Tx(1:5,1)=x(1,1:5)
Tx(1:5,2)=x(2,1:5)
! xTx = MATMUL(x,Tx)
xTx(1,1)=x(1,1)*Tx(1,1)+x(1,2)*Tx(2,1)+x(1,3)*Tx(3,1)+x(1,4)*Tx(4,1)+x(1,5)*Tx(5,1)
xTx(1,2)=x(1,1)*Tx(1,2)+x(1,2)*Tx(2,2)+x(1,3)*Tx(3,2)+x(1,4)*Tx(4,2)+x(1,5)*Tx(5,2)
xTx(2,1)=x(2,1)*Tx(1,1)+x(2,2)*Tx(2,1)+x(2,3)*Tx(3,1)+x(2,4)*Tx(4,1)+x(2,5)*Tx(5,1)
xTx(2,2)=x(2,1)*Tx(1,2)+x(2,2)*Tx(2,2)+x(2,3)*Tx(3,2)+x(2,4)*Tx(4,2)+x(2,5)*Tx(5,2)
det = xTx(1,1)*xTx(2,2) - xTx(1,2)*xTx(2,1)
a(1,1) = xTx(2,2)/det
a(1,2) = -xTx(1,2)/det
a(2,1) = -xTx(2,1)/det
a(2,2) = xTx(1,1)/det
! aTx = MATMUL(Tx,a)
aTx(1,1)=Tx(1,1)*a(1,1)+Tx(1,2)*a(2,1)
aTx(2,1)=Tx(2,1)*a(1,1)+Tx(2,2)*a(2,1)
aTx(3,1)=Tx(3,1)*a(1,1)+Tx(3,2)*a(2,1)
aTx(4,1)=Tx(4,1)*a(1,1)+Tx(4,2)*a(2,1)
aTx(5,1)=Tx(5,1)*a(1,1)+Tx(5,2)*a(2,1)
aTx(1,2)=Tx(1,1)*a(1,2)+Tx(1,2)*a(2,2)
aTx(2,2)=Tx(2,1)*a(1,2)+Tx(2,2)*a(2,2)
aTx(3,2)=Tx(3,1)*a(1,2)+Tx(3,2)*a(2,2)
aTx(4,2)=Tx(4,1)*a(1,2)+Tx(4,2)*a(2,2)
aTx(5,2)=Tx(5,1)*a(1,2)+Tx(5,2)*a(2,2)
! beta = MATMUL(y,aTx)
beta(1) = y(1)*aTx(1,1)+y(2)*aTx(2,1)+y(3)*aTx(3,1)+y(4)*aTx(4,1)+y(5)*aTx(5,1)
beta(2) = y(1)*aTx(1,2)+y(2)*aTx(2,2)+y(3)*aTx(3,2)+y(4)*aTx(4,2)+y(5)*aTx(5,2)
! beta2 = MATMUL(y2,aTx)
beta2(1) = y2(1)*aTx(1,1)+y2(2)*aTx(2,1)+y2(3)*aTx(3,1)+y2(4)*aTx(4,1)+y2(5)*aTx(5,1)
beta2(2) = y2(1)*aTx(1,2)+y2(2)*aTx(2,2)+y2(3)*aTx(3,2)+y2(4)*aTx(4,2)+y2(5)*aTx(5,2)
!============== finished least squares regression on doy**3 ==============
if ((di .lt. 3) .or. (di .gt. 72)) tt = tt + 365.
eqtime=(beta(1) + beta(2)*tt**3)/60.
decang=beta2(1) + beta2(2)*tt**3
latsun=decang
ut=time
noon=12.-lon/15. ! universal time of noon
lonsun=-15.*(ut-12.+eqtime)
t0=(90.-lat)*pi/180.
t1=(90.-latsun)*pi/180.
p0=lon*pi/180.
p1=lonsun*pi/180.
zz=cos(t0)*cos(t1)+sin(t0)*sin(t1)*cos(p1-p0)
! zz2=sin(t0)*sin(t1)+cos(t0)*cos(t1)*cos(p1-p0)
xx=sin(t1)*sin(p1-p0)
yy=sin(t0)*cos(t1)-cos(t0)*sin(t1)*cos(p1-p0)
sun_zenith=90-acos(zz)/(pi/180)
sun_azimuth=atan2(xx,yy)/(pi/180)
return
end subroutine zensun