``` 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
```