SUBROUTINE SURFAC( ALBEDO, DELM0, FBEAM, HLPR, LAMBER, MI, MAZIM, 1,7
     &                   MXCMU, MXUMU, NN, NUMU, NSTR, ONLYFL, UMU,
     &                   USRANG, YLM0, YLMC, YLMU, BDR, EMU, BEM, RMU,
     &                   SQT )

c       Specifies user's surface bidirectional properties, STWJ(21)
c
c
c   I N P U T     V A R I A B L E S:
c
c       DELM0  :  Kronecker delta, delta-sub-m0
c
c       HLPR   :  Legendre moments of surface bidirectional reflectivity
c                 (with 2K+1 factor included)
c
c       MAZIM  :  Order of azimuthal component
c
c       NN     :  Order of double-Gauss quadrature (NSTR/2)
c
c       YLM0   :  Normalized associated Legendre polynomial
c                 at the beam angle
c
c       YLMC   :  Normalized associated Legendre polynomials
c                 at the quadrature angles
c
c       YLMU   :  Normalized associated Legendre polynomials
c                 at the user angles
c
c       SQT(k) :  Square root of k
c
c       (remainder are DISORT input variables)
c
c
c    O U T P U T     V A R I A B L E S:
c
c       BDR :  Surface bidirectional reflectivity (computational angles)
c
c       RMU :  Surface bidirectional reflectivity (user angles)
c
c       BEM :  Surface directional emissivity (computational angles)
c
c       EMU :  Surface directional emissivity (user angles)
c
c
c    I N T E R N A L     V A R I A B L E S:
c
c       DREF      Directional reflectivity
c
c       NMUG   :  Number of angle cosine quadrature points on (0,1) for
c                 integrating bidirectional reflectivity to get
c                 directional emissivity (it is necessary to use a
c                 quadrature set distinct from the computational
c                 angles, because the computational angles may not be
c                 dense enough--NSTR may be too small--to give an
c                 accurate approximation for the integration).
c
c       GMU    :  The NMUG angle cosine quadrature points on (0,1)
c       GWT    :  The NMUG angle cosine quadrature weights on (0,1)
c
c       YLMG   :  Normalized associated Legendre polynomials
c                 at the NMUG quadrature angles
c
c   Called by- DISORT
c   Calls- QGAUSN, LEPOLY, ZEROIT, ERRMSG
c +-------------------------------------------------------------------+

c     .. Parameters ..

      INTEGER   NMUG, MAXSTR
c                             ** CAUTION:  Do not increase MAXSTR
c                             **           without checking if this
c                             **           would require a larger
c                             **           dimension for SQT
      PARAMETER ( NMUG = 10, MAXSTR = 100 )
c     ..
c     .. Scalar Arguments ..

      LOGICAL   LAMBER, ONLYFL, USRANG
      INTEGER   MAZIM, MI, MXCMU, MXUMU, NN, NSTR, NUMU
      REAL      ALBEDO, DELM0, FBEAM
c     ..
c     .. Array Arguments ..

      REAL      BDR( MI, 0:MI ), BEM( MI ), EMU( MXUMU ),
     &          HLPR( 0:MXCMU ), RMU( MXUMU, 0:MI ), UMU( * ),
     &          YLM0( 0:MXCMU ), YLMC( 0:MXCMU, MXCMU ),
     &          YLMU( 0:MXCMU, MXUMU ), SQT( * )
c     ..
c     .. Local Scalars ..

      LOGICAL   PASS1
      INTEGER   IQ, IU, JG, JQ, K
      REAL      DREF, SGN, SUM
c     ..
c     .. Local Arrays ..

      REAL      GMU( NMUG ), GWT( NMUG ), YLMG( 0:MAXSTR, NMUG )
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, LEPOLY, QGAUSN, ZEROIT
c     ..
      SAVE      PASS1, GMU, GWT, YLMG
      DATA      PASS1 / .True. /


      IF( PASS1 ) THEN

         PASS1  = .FALSE.

         CALL QGAUSN( NMUG, GMU, GWT )

         CALL LEPOLY( NMUG, 0, MAXSTR, MAXSTR, GMU, SQT, YLMG )

c                       ** Convert Legendre polys. to negative GMU
         SGN  = - 1.0

         DO 20 K = 0, MAXSTR

            SGN  = - SGN

            DO 10 JG = 1, NMUG
               YLMG( K, JG ) = SGN*YLMG( K, JG )
   10       CONTINUE

   20    CONTINUE

      END IF


      CALL ZEROIT( BDR, MI*( MI + 1 ) )
      CALL ZEROIT( BEM, MI )

      IF( LAMBER .AND. MAZIM.EQ.0 ) THEN

         DO 40 IQ = 1, NN

            BEM( IQ ) = 1.- ALBEDO

            DO 30 JQ = 0, NN
               BDR( IQ, JQ ) = ALBEDO
   30       CONTINUE

   40    CONTINUE


      ELSE IF( .NOT.LAMBER ) THEN
c                                  ** Compute surface bidirectional
c                                  ** properties at computational angles
         DO 80 IQ = 1, NN

            DO 60 JQ = 1, NN

               SUM  = 0.0
               DO 50 K = MAZIM, NSTR - 1
                  SUM  = SUM + HLPR( K )*YLMC( K, IQ )*
     &                         YLMC( K, JQ + NN )
   50          CONTINUE

               BDR( IQ, JQ ) = ( 2.- DELM0 )*SUM

   60       CONTINUE


            IF( FBEAM.GT.0.0 ) THEN

               SUM  = 0.0
               DO 70 K = MAZIM, NSTR - 1
                  SUM  = SUM + HLPR( K )*YLMC( K, IQ )*YLM0( K )
   70          CONTINUE

               BDR( IQ, 0 ) = ( 2.- DELM0 )*SUM

            END IF

   80    CONTINUE


         IF( MAZIM.EQ.0 ) THEN

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

c                              ** Integrate bidirectional reflectivity
c                              ** at reflection polar angles CMU and
c                              ** incident angles GMU to get
c                              ** directional emissivity at
c                              ** computational angles CMU.
            DO 110 IQ = 1, NN

               DREF  = 0.0

               DO 100 JG = 1, NMUG

                  SUM  = 0.0
                  DO 90 K = 0, NSTR - 1
                     SUM  = SUM + HLPR( K )*YLMC( K, IQ )*
     &                            YLMG( K, JG )
   90             CONTINUE

                  DREF  = DREF + 2.*GWT( JG )*GMU( JG )*SUM

  100          CONTINUE

               BEM( IQ ) = 1.- DREF

  110       CONTINUE

         END IF

      END IF
c                                       ** Compute surface bidirectional
c                                       ** properties at user angles

      IF( .NOT.ONLYFL .AND. USRANG ) THEN

         CALL ZEROIT( EMU, MXUMU )
         CALL ZEROIT( RMU, MXUMU*( MI + 1 ) )

         DO 180 IU = 1, NUMU

            IF( UMU( IU ).GT.0.0 ) THEN

               IF( LAMBER .AND. MAZIM.EQ.0 ) THEN

                  DO 120 IQ = 0, NN
                     RMU( IU, IQ ) = ALBEDO
  120             CONTINUE

                  EMU( IU ) = 1.- ALBEDO


               ELSE IF( .NOT.LAMBER ) THEN

                  DO 140 IQ = 1, NN

                     SUM  = 0.0
                     DO 130 K = MAZIM, NSTR - 1
                        SUM  = SUM + HLPR( K )*YLMU( K, IU )*
     &                               YLMC( K, IQ + NN )
  130                CONTINUE

                     RMU( IU, IQ ) = ( 2.- DELM0 )*SUM

  140             CONTINUE


                  IF( FBEAM.GT.0.0 ) THEN

                     SUM  = 0.0
                     DO 150 K = MAZIM, NSTR - 1
                        SUM  = SUM + HLPR( K )*YLMU( K, IU )*YLM0( K )
  150                CONTINUE

                     RMU( IU, 0 ) = ( 2.- DELM0 )*SUM

                  END IF


                  IF( MAZIM.EQ.0 ) THEN

c                               ** Integrate bidirectional reflectivity
c                               ** at reflection angles UMU and
c                               ** incident angles GMU to get
c                               ** directional emissivity at
c                               ** user angles UMU.
                     DREF  = 0.0

                     DO 170 JG = 1, NMUG

                        SUM  = 0.0
                        DO 160 K = 0, NSTR - 1
                           SUM  = SUM + HLPR( K )*YLMU( K, IU )*
     &                                  YLMG( K, JG )
  160                   CONTINUE

                        DREF  = DREF + 2.*GWT( JG )*GMU( JG )*SUM

  170                CONTINUE

                     EMU( IU ) = 1.- DREF

                  END IF

               END IF

            END IF

  180    CONTINUE

      END IF


      RETURN
      END