! new2.f90
real pi, radian, L, julian, csoff, Ds, arg, Etime, tciv
real zenith, azimuth, hourset, tsunset, tsunrise
pi=3.141593
radian=180./pi
write(6,*)'Enter latitude (deg. N),Julian day, and offset'
write(6,*)' of civil time from local solar time as decimal hours'
read(5,*)L,julian,csoff
L=L/57.29
! write(6,'("Latitude = ",f7.3," deg., JD = ",i3,", t offset = ",f5.2," h")')&
! & L, julian, csoff
write(6,'("t(civil) zenith azimuth")')
! Compute some quantities approximately fixed over the day
Ds=asin(sin(23.44/radian)*cos(2.*pi*(julian-172.)/365.)) ! solar declin., in radians
! Correct for the equation of time
! t(solar) = t(mean solar) + E
! sine-wave argument is omega(annual)*t(as Julian day) = 2 pi jd / 365
arg=2.*pi*julian/365.
! Eorbital=-7.8 min * sin(arg), Etilt = -10 min * sin(2*arg)
Etime=-0.13*sin(arg)-0.167*sin(2.*arg)
do tciv=4.,20.,0.5
call solar_subr(julian,tciv,L,csoff,zenith,azimuth,Ds,Etime)
write(6,'(f6.2,3xf8.2,2x,f8.2)')tciv,zenith,azimuth
enddo
arg=-tan(L)*tan(Ds)
hourset=acos(arg)
tsunset =12*(1+hourset/pi)-csoff-Etime
! Above is possibly not corrected for the eq. of time
tsunrise=24.-tsunset
write(6,'("Sunrise is at ",f6.2,", sunset at ",f6.2)')tsunrise,tsunset
stop
end
!
!
subroutine solar_subr(julian,tciv,L,csoff,zenithdeg,azimuthdeg,Ds,Etime)
! Subroutine to determine cosine of solar zenith angle (same as sine of elevation angle)
! at any time of day, plus times of sunrise and sunset
! This subroutine is used in the analysis of our weather data
! It was modified from solar.angles.total.3.f written in 2013, which was, in turn
! modified from a program written originally in 1999
! mdt, tsunrise, and tsunset are all given in decimal time from midnight
real julian, tciv, L, csoff, zenithdeg, azimuthdeg, Ds, Etime
real pi, radian, sinL, cosL, time, hangle, sinalphas, alphas, zenith
real cosDx, sinha, cosha, cosa, sinphis, phis, dsinphisdh
pi=3.1415926
radian=180./pi
! Latitude of Las Cruces (our home) = 32.272o, from UTM 13S E=0337110, N=3571896
! and longitude = -106.730
sinL=sin(L)
cosL=cos(L)
! csoff =-0.06 ! Civil time is 0.06 hours earlier than solar time in Las Cruces
! Convert local civil time (MDT) to local standard time and then to
! local solar time (ignoring 'equation of time'
! For Las Cruces, at -106.73o, this seems to be -0.06 (-1.73/15)
time=tciv+Etime+csoff
! Convert to hour angle
hangle=(time-12.)*pi/12.
sinalphas=sinL*sin(Ds)+cosL*cos(Ds)*cos(hangle)
alphas=asin(sinalphas)
zenith=pi/2.-alphas
zenithdeg=zenith*radian
! Compute the azimuth
! I didn't need to compute solar azimuth (though the Preston building shades
! the radiation sensor in summer at late hours; phis would be a bit useful then
cosDs=cos(Ds)
sinha=sin(hangle)
cosha=cos(hangle)
cosa=cos(alphas)
sinphis=-cosDs*sinha/cosa
phis=asin(sinphis)
! Get it into the correct quadrant
dsinphisdh=(cosa**2*cosha-sinha**2*sinalphas*cosL*cosDs)/cosa**3
if(sinalphas.gt.0)then !We're in daytime, only time the correc. is relevant
if(dsinphisdh.lt.0.)then
if(hangle.lt.0.)phis=pi-phis
if(hangle.gt.0.)phis=-pi-phis
endif
endif
azimuthdeg=phis*radian
return
end