SUBROUTINE gfluxquad(oz, sz, pr, rs, lat, nquad, quad, success)
c
c PURPOSE
c   calculates weighted quadrature (integral) of the global flux at the
c   earth's surface, given the parameters (oz, sz, pr, rs, lat), the
c   solar spectrum (in common /solspec/), and the weighting functions
c   described in common /weight/.
c
c CALLING PARAMETERS
c   INPUT
c     oz	ozone (DU)
c     sz	solar zenith angle (deg)
c     pr        terrain pressure (atm)
c     rs        surface reflectivity (fraction)
c     lat       latitude (deg)
c     nquad     dimension of quad
c
c OUTPUT
c     quad      is the array (must be of dimension .GE. the number 
c               of weighting functions) that will catch the output.
c     success   flags whether the subroutine executed successfully.

      INCLUDE "params.inc"
c---------------------------------------------------------------------
c Calling parameters
c
      REAL*8 oz, sz, pr, rs, lat, quad(0:nquad-1)
      INTEGER*4 nquad
      LOGICAL success


c---------------------------------------------------------------------
c Common blocks
c
      INCLUDE "weights.cmn"
      INCLUDE "solspec.cmn"
      INCLUDE "rttbls_cs.cmn"
      INCLUDE "cloud.cmn"

c---------------------------------------------------------------------
c Local variables
c
      INTEGER*4 ifn, ilam, ilam0, ilam1, ix
      REAL*8 flux(0:n_waveln-1), sum

      IF (nquad .LT. n_w_fns) GOTO 900

      DO ifn= 0, n_w_fns-1

        npts= w_ix0(ifn+1)-w_ix0(ifn)
        ilam0= w_ilam0(ifn)
        ilam1= ilam0 + npts - 1 

d       WRITE(9,92) oz, sz, pr, rs, lat
d  92   FORMAT('oz: ',f10.3,'sz: ',f10.3,'pr: ',f10.3,'rs: ',f10.3,
d    &       'lat: ', f10.3) 

        CALL gfluxinterp(oz, sz, pr, rs, lat, ilam0, ilam1,
     &         flux, success)

        IF (.NOT. success) RETURN

d       WRITE(9,90) ifn, w_ix0(ifn), npts, ilam0, ilam1,
d    &    w_zero0(ifn), waveln(ilam0), waveln(ilam1), w_zero1(ifn)
d  90 format(/'ifn:',i3,'  w_ix0:',i4,'  npts:',i5,'  ilam0,1:',2i5,
d    &       '  lams:', 4f8.3)

        DO ix=0, npts-1

d         WRITE(9,91)ix, waveln(ilam0+ix),flux(ix),solspec(ilam0+ix),
d    &      weight(w_ix0(ifn)+ix), cldfac(w_ix0(ifn)+ix)
d  91 format(3x,i5,f8.3,5(1pe14.6))

        END DO

        ix=0
        flux(ix)= flux(ix) *
     &            solspec(ilam0+ix) *
     &            weight(w_ix0(ifn)+ix) *
     &            cldfac(w_ix0(ifn)+ix)
      
        sum= (waveln(ilam0) - w_zero0(ifn))*flux(0)

        DO ix=1, npts-1
          flux(ix)= flux(ix) *
     &              solspec(ilam0+ix) *
     &              weight(w_ix0(ifn)+ix) *
     &              cldfac(w_ix0(ifn)+ix)

          sum= sum + (waveln(ilam0+ix)-waveln(ilam0+ix-1))*
     &             (flux(ix-1) + flux(ix))
        END DO

        sum= sum + (w_zero1(ifn) - waveln(ilam1))*flux(npts-1)
        quad(ifn)= sum/2.d0

d       WRITE(9,93) quad(ifn)
d  93 FORMAT(' *** Integrated flux: ',e14.6)

      END DO

      success= .TRUE.

      RETURN

  900 CONTINUE
      WRITE(6,1) nquad, n_w_fns
      success= .FALSE.
      RETURN

    1 FORMAT('!!! gfluxquad: supplied output array has dimension, ',
     &    I4,'.  It must have dimension >= ',I4)

      END