! 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