c$Author: fredfox $
c$Date: 2002-05-08 18:47:45 $
c$Revision: 1.2 $
c$Source: /weru/cvs/weps/weps.src/util/solar/hourangle.for,v $

      real function   hourangle(dlat, dec, riseangle)

c     + + + PURPOSE + + +
c     This function calculates the hour angle (degrees) 
c     of sunrise (-), sunset (+) based on the declination of the earth

c     + + + KEYWORDS + + +
c     sunrise sunset hourangle

c     + + + ARGUMENT DECLARATIONS + + +
      real dlat
      real dec
      real riseangle

c     + + + ARGUMENT DEFINITIONS + + +
c     dlat   - Latitude of the site, degrees (north > 0, south < 0)
c     dec    - declination of earth wrt the sun (degrees)
c     riseangle - angle of earths rotation where sunrise occurs
c                 this varies depending on whether you are calculating
c                 direct beam, civil twilight, nautical twilight or
c                 astronomical twilight hourangle

c     + + + LOCAL VARIABLES + + +
      real h, coshr
      real dlat_rad, dec_rad
      real dlat_rad_lim

      parameter( dlat_rad_lim = 1.57079 ) ! pi/2 minus a small bit
c      parameter( dlat_rad_lim = 1.570796327 ) ! pi/2

c     + + + LOCAL DEFINITIONS + + +
c     coshr   - Cosine of hour angle at sunrise
c     dlat_rad - latitude of site, converted to radians
c     dec_rad - declination of earth wrt the sun (radians)

c     + + + COMMON BLOCKS + + +
      include 'p1unconv.inc'

c     + + + END SPECIFICATIONS + + +

c     convert to radians
      dlat_rad = dlat * degtorad
      dec_rad = dec * degtorad

c     Calculate the cosine of hour angle (h) at sunset
c     To get the sunrise hour angle, take the negative.
c     Using the equation from "Solar Thermal Energy Systems,
c     Howell, Bannerot, Vliet, 1982, page 51 equation 3-4)
c     modified to account for atmospheric refraction as in
c     NOAA document (it just indicates that the sun is seen
c     before it physically is above the horizon)
c     ie. not at 90 degrees, but 90.833 degrees
c     THis expression is undefined at 90 and -90 degrees. If
c     roundoff error pushes it beyond the answer flips. Limit
c     set here to get correct answer at 90 and -90 degrees.
      dlat_rad = max( -dlat_rad_lim, min(dlat_rad_lim, dlat_rad))
      coshr = cos(riseangle*degtorad)/(cos(dlat_rad)*cos(dec_rad))
     &      - tan(dlat_rad)*tan(dec_rad)

c     check for artic circle conditions
      if( coshr.ge.1.0) then  
          hourangle = 0.0          !sunrise occurs at solar noon
      else if( coshr.le.-1.0) then
          hourangle = 180.0        !the sun is always above the horizon
      else
          hourangle = acos(coshr) * radtodeg
      end if

      return
      end