FUNCTION gfluxdiurn(oz, pr, rs, par1, par2, parmode, 
     &                                           taucl, success)
c
c PURPOSE
c  Calculate the diurnally- and spectrally-integrated, weighted irradiances
c
c CALLING PARAMETERS
c  oz		i:R8	Profile ozone value, in DU
c  pr		i:R8	Terrain pressure, in atm
c  rs		i:R8	Surface reflectivity (fractional)
c  par1		i:R8	First parameter (see below)
c  par2		i:R8	Second parameter (see below)
c  parmode	i:I4	Flag that indicates what the parameters are (vide infra)
c  taucl	i:R8	Cloud optical thickness
c  success	o:L	Success flag
c
c Notes: Parameter modes
c  This routine can be called using solar declination angle and latitude,
c  or two equivalent parameterizations (a,b) or (p,q).  The following 
c  table indicates the identities of the parameters for the three modes:
c
c 	parmode		par1		par2
c	-------		---------	--------
c	  1		soldecang	latitude	(both in degrees)
c	  2		a		b		(both dimensionless)
c	  3		p		q		(both dimensionless)
c
c  By definition, 
c	a= sin(soldecang)*sin(latitude)
c	b= cos(soldecang)*cos(latitude)
c	p= b + a = cos(soldecang - latitude)
c	q= b - a = cos(soldecang + latitude)
c
      IMPLICIT NONE

      INCLUDE "params.inc"

c----------------------------------------------------------------------
c Calling parameters
c
      REAL*8 gfluxdiurn
      REAL*8 oz, pr, rs, par1, par2, taucl
      INTEGER*4 parmode
      LOGICAL success

c----------------------------------------------------------------------
c Common Blocks
c
      INCLUDE "gauscheb.cmn"

c----------------------------------------------------------------------
c Local variables
c
      INTEGER*4 llim
      INTEGER*4 nquad /1/, i
      REAL*8    sd, sl, cd, cl, a, b, p, q, soldecang, lat
      REAL*8    z, sum, quad

      REAL*8 radeg      /57.295779513082323/    !deg-rad conversion

c----------------------------------------------------------------------
c Whatever the parmode, find a and b...
c
      IF(parmode .EQ. 1) THEN
        soldecang= par1/radeg
        sd= sin(soldecang)
        cd= cos(soldecang)

        lat= par2/radeg
        sl= sin(lat)
        cl= cos(lat)

        a= sd*sl
        b= cd*cl

      ELSE IF(parmode .EQ. 2) THEN
        lat= 0.d0
        a= par1
        b= par2

      ELSE IF(parmode .EQ. 3) THEN
        lat= 0.d0
        p= par1
        q= par2

        a= (p - q)/2.d0
        b= (p + q)/2.d0

      ELSE
        WRITE(6,*)' !!! gfluxdiurn !!! - called with invalid parmode:',
     &               parmode
        success=.FALSE.
        gfluxdiurn= -999.
        RETURN

      END IF

c------------------------------------------------------------------------
c On with the calculation...

      llim= int(nscl*acos(-a/b)-.5)

      sum= 0.d0
      DO i=0,llim
        z= acos(b*y(i)+a)*radeg
        CALL ld_cloud_ccf(taucl, rs, pr, z, success)
           IF(.NOT. success) GOTO 901
        CALL gfluxquad(oz, z, pr, rs, lat, nquad, quad, success)
           IF(.NOT. success) GOTO 902
        sum= sum + quad
      END DO

      gfluxdiurn= t_fac*sum/nscl
      RETURN

  901 CONTINUE
      WRITE(6,*)' !!! gfluxdiurn !!! call to ld_cloud_ccf failed.'
      gfluxdiurn= 0.d0
      RETURN

  902 CONTINUE
      WRITE(6,*)' !!! gfluxdiurn !!! call to gfluxquad failed.'
      gfluxdiurn= 0.d0
      RETURN

      END