REAL FUNCTION DREF( MU, HL, NSTR ) 1,3

c        Exact flux albedo for given angle of incidence, given
c        a bidirectional reflectivity characterized by its
c        Legendre coefficients ( NOTE** these will only agree
c        with bottom-boundary albedos calculated by DISORT in
c        the limit as number of streams go to infinity, because
c        DISORT evaluates the integral 'CL' only approximately,
c        by quadrature, while this routine calculates it exactly.)
c
c  INPUT :   MU     Cosine of incidence angle
c
c            HL     Legendre coefficients of bidirectional reflectivity
c
c          NSTR     Number of elements of HL to consider
c
c
c  INTERNAL VARIABLES (P-sub-L is the L-th Legendre polynomial) :
c
c       CL      Integral from 0 to 1 of  MU * P-sub-L(MU)
c                   (vanishes for  L = 3, 5, 7, ... )
c       PL      P-sub-L
c       PLM1    P-sub-(L-1)
c       PLM2    P-sub-(L-2)
c
c   Called by- CHEKIN
c   Calls- ERRMSG
c +-------------------------------------------------------------------+

c     .. Parameters ..

      INTEGER   MAXTRM
      PARAMETER ( MAXTRM = 100 )
c     ..
c     .. Scalar Arguments ..

      INTEGER   NSTR
      REAL      MU
c     ..
c     .. Array Arguments ..

      REAL      HL( 0:NSTR )
c     ..
c     .. Local Scalars ..

      LOGICAL   PASS1
      INTEGER   L
      REAL      CL, PL, PLM1, PLM2
c     ..
c     .. Local Arrays ..

      REAL      C( MAXTRM )
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC MOD
c     ..
      SAVE      PASS1, C
      DATA      PASS1 / .True. /
c     ..


      IF( PASS1 ) THEN

         PASS1  = .FALSE.
         CL     = 0.125
         C( 2 ) = 10.*CL

         DO 10 L = 4, MAXTRM, 2
            CL     = - CL*( L - 3 ) / ( L + 2 )
            C( L ) = 2.*( 2*L + 1 )*CL
   10    CONTINUE

      END IF


      IF( NSTR.LT.2 .OR. ABS(MU).GT.1.0 )
     &    CALL ERRMSG( 'DREF--input argument error(s)',.True. )

      IF( NSTR.GT.MAXTRM )
     &    CALL ERRMSG( 'DREF--parameter MAXTRM too small',.True. )


      DREF  = HL( 0 ) - 2.*HL( 1 )*MU
      PLM2  = 1.0
      PLM1  = - MU

      DO 20 L = 2, NSTR - 1
c                                ** Legendre polynomial recurrence

         PL = (
         IF( MOD( L,2 ).EQ.0 ) DREF   = DREF + C( L )*HL( L )*PL

         PLM2  = PLM1
         PLM1  = PL

   20 CONTINUE

      IF( DREF.LT.0.0 .OR. DREF.GT.1.0 )
     &    CALL ERRMSG( 'DREF--albedo value not in (0,1)',.False. )

      RETURN
      END