c$Author: fredfox $
c$Date: 2002-05-09 14:08:42 $
c$Revision: 1.2 $
c$Source: /weru/cvs/weps/weps.src/util/solar/radext.for,v $

      real function radext(idoy, bmalat)

c     + + + purpose + + +
c     this subroutine estimates the incoming extraterrestial radiation 
c     for a given location (Mj/m^2/day)

c     + + + key words + + +
c     radiation, solar, extraterrestrial

c     + + + argument declarations + + +
      integer idoy
      real bmalat

c     + + + argument definitions + + +
c     idoy     - julian day of year, 1-366
c     bmalat   - latitude of the site, degrees

c     + + + COMMON BLOCK + + +
      include 'p1const.inc'
      include 'p1unconv.inc'
      include 'p1solar.inc'

c     + + + local variables + + +
      real rlat
      real dec
      real rdec
      real dr
      real ws
      real ra1
      real ra2

c     + + + local definitions + + +
c     rlat   - latitude (radians)
c     rdec   - declination (radians)
c     dr     - direct radiation variation with distance from sun of earth orbit
c     ws     - sunset hour angle (radians)
c     ra1, ra2 - intermediate calculations

c     + + + parameters + + +

      real gsc

      parameter (gsc = 0.08202)

c     gsc - solar_constant in Mj/m^2-min (0.08202 Mj/m^2-min = 1367 W/m^2)

c     + + + data initialization + + +

c     + + + FUNCTION DECLARATIONS + + +
      real declination
      real hourangle

c     + + + end specifications + + +

c     convert to radians for trig functions
      rlat = bmalat * degtorad
      rdec  = declination(idoy) * degtorad

c     compute factor for variable distance from sun along orbital path
      dr = 1 + 0.033*cos(2*pi*idoy/365)							!h-21

      ws = hourangle(bmalat, dec, beamrise ) * degtorad
      ra1 = ((24.0*60.0)/pi)*gsc*dr									!h-20(a)
      ra2 = (ws*sin(rlat)*sin(rdec))+(cos(rlat)*cos(rdec)*sin(ws))	!h-20(b)
      radext = ra1*ra2												!h-20

      return
      end