c########################################################
c
      subroutine sun(yyyy,doy,timeday,latit,long,
     &sinele,sunrise,sunset)
c
c    This subroutine calculates month (zmm), day of month (dddd),
c  local solar noon, sun rise, sun set, and solar declination 
c from year (yyyy), day of year (doy), and geographic location.
c 
c =======================Input Variables===================
c yyyy:   year
c doy:    day of year
c timeday: timeday in GMT
c latit:  latitude,  degrees, north possitive
c long: longitude, degrees, east positive
c
c =======================Output Variables==================
c
c zmm:     month
c dddd:    date (day of month)
c locnoon: local apparent noon time, hours, GMT
c soldec:  solar declination in radians
c sunrise: sun rise time hours GMT
c sunset:  sun set time hpurs GMT
c sinele: sine of solar elevation angle
c 
c =======================Subroutines Called=================
c subroutines: none
c functions: none
c
      implicit double precision(a-h,l,o-z)
c
      longit=-long
      ia=1889
      if (mod(idint(yyyy),4) .eq. 0) ia=1523
      ib=idint((doy+dble(ia)-122.1d0)/365.25d0)
      ic=idint(doy)+ia-idint(365.25d0*dble(ib))
      ie=idint(dble(ic)/30.6001d0)
      if (dble(ie) .lt. 13.5d0) then
        zmm=dble(ie-1)
      else
        zmm=dble(ie-13)
      end if
      dddd=dble(ic-idint(30.6001d0*dble(ie)))
      if (zmm .gt. 2.1d0) then
        y=yyyy
        zm=zmm
      else
        y=yyyy-1.0d0
        zm=zmm+12.0d0
      end if
      ia=idint(y/100.0d0)
      ib=2-ia+idint(dble(ia)/4.0d0)
      fjd=dble(idint(365.25d0*y)+idint(30.6001d0*
     &    (zm+1.0d0)))+dddd+1720994.5d0+dble(ib)
      capt=(fjd-2415020.0d0)/36525.0d0
      capl=279.69668d0+36000.76892d0*capt+
     &     0.0003025d0*capt*capt
      capm=358.47583d0+35999.04975d0*capt-0.00015*capt*
     &     capt-3.3d-6*capt*capt*capt
      smalle=0.01675104d0-4.18d-5*capt-1.26d-7*capt*capt
      epis=23.452294d0-0.0130125d0*capt-1.64d-6*capt*capt+
     &     5.03d-7*capt*capt*capt
      fact=3.14159265d0/180.0d0
      yepis=dtan(epis*fact/2.0d0)*dtan(epis*fact/2.0d0)
      cape=yepis*dsin(2.0d0*capl*fact)-2.0d0*smalle*dsin
     &     (capm*fact)+4.0d0*smalle*yepis*dsin(capm*fact)*
     &     dcos(2.0d0*capl*fact)-0.5d0*yepis*yepis*dsin(
     &     4.0d0*capl*fact)-1.25d0*smalle*smalle*
     &     dsin(2.0d0*capm*fact)
      eqtime=4.0d0*cape/fact
      capc=(1.91946d0-0.004789d0*capt-1.4d-5*capt*capt)*
     &      dsin(capm*fact)+(0.020094d0-1.0d-4*capt)*dsin(2.0d0*
     &      capm)+2.93d-4*dsin(3.0d0*capm)
      theta=capl+capc
      soldec=dasin(dsin(epis*fact)*dsin(theta*fact))

      locnoon=12.0d0-eqtime/60.0d0+4.0d0*longit/60.0d0
c      decmax=23.5d0*fact
c      sols=4141.0d0/24.0d0+dble(mod(idint(yyyy)+3,4))*0.25d0
c      season=(doy-sols)/365.2d0
c      soldec=decmax*dcos(2.0d0*3.14159265d0*season)
      hac = -dtan(latit*fact)*dtan(soldec)-
     &     0.01454d0/(dcos(latit*fact)*dcos(soldec))
      hac  = dmin1(hac,1.0d0)
      hac  = dmax1(hac,-1.0d0)
c                                                                               
c-----------------------------------------------------------------------        
c     h is the half-day length (in radians)                                     
c-----------------------------------------------------------------------        
c                                                                               
      h   = dacos(hac)                                                           
      sunrise = locnoon-(h/(15.0d0*fact))
      sunset  = locnoon+(h/(15.0d0*fact))

      if (sunset .gt. 24.0d0)then
        sunset=sunset-24.0d0
      end if
      if (sunrise.lt.0.0d0)then
        sunrise=24.0d0+sunrise
      end if

      fact=3.14159265d0/180.0d0
      hourang=15.0d0*(timeday-locnoon)*fact
      sinele=dsin(latit*fact)*dsin(soldec)+
     &       dcos(latit*fact)*dcos(soldec)*dcos(hourang)
      sinele=dmin1(sinele,1.0d0)
      sinele=dmax1(sinele,-1.0d0)
      return
      end