parameter (mxtim=100)
      real tary(mxtim)
      real*8 cozena, dlat, dlon, dday, dtime, val
      REAL lon, lat, musun
      real a,e,az el 
C test program for solar position codes

      year=2003.0
      day=182.0
      time=10.5 
      lon=20.0
      lat=50.0
      
      pi=4.*atan(1.)
      dtor=pi/180.
      
      
      print*, ' year day time lat lon ',year,day,time,lat,lon
      
      call sunae1(year,day,time,lat,lon,az,el)
      print*, az,el


      STOP
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine sunae1(year,day,hour,lat,long,az,el)
      implicit real(a-z)
C Purpose:
C Calculates azimuth and elevation of sun
C
C References:
C (a) Michalsky, J. J., 1988, The Astronomical Almanac's algorithm for
C approximate solar position (1950-2050), Solar Energy, 227---235, 1988
C
C (b) Spencer, J. W., 1989, Comments on The Astronomical
C Almanac's algorithm for approximate solar position (1950-2050)
C Solar Energy, 42, 353
C
C Input:
C year - the year number (e.g. 1977)
C day  - the day number of the year starting with 1 for
C        January 1
C time - decimal time. E.g. 22.89 (8.30am eastern daylight time is
C        equal to 8.5+5(hours west of Greenwich) -1 (for daylight savings
C        time correction
C lat -  local latitude in degrees (north is positive)
C lon -  local longitude (east of Greenwich is positive.
C                         i.e. Honolulu is 15.3, -157.8)
C
C Output:
C a - azimuth angle of the sun (measured east from north 0 to 360)
C e - elevation of the sun
C
C Spencer correction introduced and 3 lines of Michalsky code
C commented out (after calculation of az)
C
      pi=4.*atan(1.)
      twopi=2.*pi
      rad=pi/180.
C get the current Julian date
      delta=year-1949.
      leap=aint(delta/4.)
      jd=32916.5+delta*365.+leap+day+hour/24.
C calculate ecliptic coordinates
       time=jd-51545.0
C force mean longitude between 0 and 360 degs
       mnlong=280.460+0.9856474*time
       mnlong=mod(mnlong,360.)
       if(mnlong.lt.0) mnlong=mnlong+360.
C mean anomaly in radians between 0, 2*pi
       mnanom=357.528+0.9856003*time
       mnanom=mod(mnanom,360.)
       if(mnanom.lt.0.) mnanom=mnanom+360.
       mnanom=mnanom*rad
C compute ecliptic longitude and obliquity of ecliptic
C       eclong=mnlong+1.915*(mnanom)+0.20*sin(2.*mnanom)
       eclong=mnlong+1.915*sin(mnanom)+0.020*sin(2.*mnanom)
       eclong=mod(eclong, 360.)
       if(eclong.lt.0) eclong=eclong+360.
       oblqec=23.429-0.0000004*time
       eclong=eclong*rad
       oblqec=oblqec*rad
C calculate right ascention and declination
       num=cos(oblqec)*sin(eclong)
       den=cos(eclong)
       ra=atan(num/den)
C force ra between 0 and 2*pi
       if(den.lt.0.) then
        ra=ra+pi
       elseif(num.lt.0.) then
        ra=ra+twopi
       endif
C dec in radians
       dec=asin(sin(oblqec)*sin(eclong))
C calculate Greenwich mean sidereal time in hours
       gmst=6.697375+0.0657098242*time+hour
C hour not changed to sidereal sine "time" includes the fractional day
       gmst=mod(gmst,24.)
       if(gmst.lt.0.) gmst=gmst+24.
C calculate local mean sidereal time in radians
       lmst=gmst+long/15.
       lmst=mod(lmst,24.)
       if(lmst.lt.0.) lmst=lmst+24.
       lmst=lmst*15.*rad
C calculate hour angle in radians between -pi, pi
        ha = lmst -ra
        if(ha.lt.-pi) ha=ha+twopi
        if(ha.gt.pi) ha=ha-twopi
        lat=lat*rad
C calculate azimuth and elevation
        el=asin(sin(dec)*sin(lat)+cos(dec)*cos(lat)*cos(ha))
        az=asin(-cos(dec)*sin(ha)/cos(el))
C add J. W. Spencer code (next 5 lines)
        if(sin(dec) - sin(el)*sin(lat) .ge. 0. ) then
         if(sin(az) .lt. 0.) az=az+twopi
        else
         az=pi-az
        endif
C end Spencer's corrections
ccc        elc=asin(sin(dec)/sin(lat))
ccc        if(el.ge.elc) az=pi-az
ccc        if(el.le.elc  .and. ha.gt.0.) az=twopi+az

C this puts azimuth between 0 and 2*pi radians
C calculate refraction correction for US stand. atm.
        el=el/rad
        if(el.gt.-0.56) then
           refrac=3.51561*(0.1594+0.0196*el+0.00002*el**2)/
     $            (1.+0.505*el+0.0845*el**2)
        else
           refrac=0.56
        endif
        el=el+refrac
        az=az/rad
        lat=lat/rad
C
        return
        end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc