SUBROUTINE SETDIS( CMU, CWT, DELTAM, DTAUC, DTAUCP, EXPBEA, FBEAM, 1,2
     &                   FLYR, GL, HL, HLPR, IBCND, LAMBER, LAYRU,
     &                   LYRCUT, MAXUMU, MAXCMU, MXCMU, NCUT, NLYR,
     &                   NTAU, NN, NSTR, PLANK, NUMU, ONLYFL, OPRIM,
     &                   PMOM, SSALB, TAUC, TAUCPR, UTAU, UTAUPR, UMU,
     &                   UMU0, USRTAU, USRANG )

c          Perform miscellaneous setting-up operations
c
c    INPUT :  all are DISORT input variables (see DOC file)
c
c
c    O U T P U T     V A R I A B L E S:
c
c       NTAU,UTAU   if USRTAU = FALSE (defined in DISORT.doc)
c       NUMU,UMU    if USRANG = FALSE (defined in DISORT.doc)
c
c       CMU,CWT     computational polar angles and
c                   corresponding quadrature weights
c
c       EXPBEA      transmission of direct beam
c
c       FLYR        truncated fraction in delta-M method
c
c       GL          phase function Legendre coefficients multiplied
c                   by (2L+1) and single-scatter albedo
c
c       HLPR        Legendre moments of surface bidirectional
c                   reflectivity, times 2K+1
c
c       LAYRU       Computational layer in which UTAU falls
c
c       LYRCUT      flag as to whether radiation will be zeroed
c                   below layer NCUT
c
c       NCUT        computational layer where absorption
c                   optical depth first exceeds  ABSCUT
c
c       NN          NSTR / 2
c
c       OPRIM       delta-M-scaled single-scatter albedo
c
c       TAUCPR      delta-M-scaled optical depth
c
c       UTAUPR      delta-M-scaled version of  UTAU
c
c   Called by- DISORT
c   Calls- QGAUSN, ERRMSG
c ---------------------------------------------------------------------

c     .. Scalar Arguments ..

      LOGICAL   DELTAM, LAMBER, LYRCUT, ONLYFL, PLANK, USRANG, USRTAU
      INTEGER   IBCND, MAXCMU, MAXUMU, MXCMU, NCUT, NLYR, NN, NSTR,
     &          NTAU, NUMU
      REAL      FBEAM, UMU0
c     ..
c     .. Array Arguments ..

      INTEGER   LAYRU( * )
      REAL      CMU( MXCMU ), CWT( MXCMU ), DTAUC( * ), DTAUCP( * ),
     &          EXPBEA( 0:* ), FLYR( * ), GL( 0:MXCMU, * ),
     &          HL( 0:MAXCMU ), HLPR( 0:MXCMU ), OPRIM( * ),
     &          PMOM( 0:MAXCMU, * ), SSALB( * ), TAUC( 0:* ),
     &          TAUCPR( 0:* ), UMU( MAXUMU ), UTAU( * ), UTAUPR( * )
c     ..
c     .. Local Scalars ..

      INTEGER   IQ, IU, K, LC, LU
      REAL      ABSCUT, ABSTAU, F
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, QGAUSN
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, EXP
c     ..
      DATA      ABSCUT / 10. /


      IF( .NOT.USRTAU ) THEN
c                              ** Set output levels at computational
c                              ** layer boundaries
         NTAU  = NLYR + 1

         DO 10 LC = 0, NTAU - 1
            UTAU( LC + 1 ) = TAUC( LC )
   10    CONTINUE

      END IF
c                        ** Apply delta-M scaling and move description
c                        ** of computational layers to local variables
      EXPBEA( 0 ) = 1.0
      TAUCPR( 0 ) = 0.0
      ABSTAU      = 0.0

      DO 40 LC = 1, NLYR

         PMOM( 0, LC ) = 1.0

         IF( ABSTAU.LT.ABSCUT ) NCUT  = LC

         ABSTAU = ABSTAU + ( 1.- SSALB( LC ) )*DTAUC( LC )

         IF( .NOT.DELTAM ) THEN

            OPRIM( LC )  = SSALB( LC )
            DTAUCP( LC ) = DTAUC( LC )
            TAUCPR( LC ) = TAUC( LC )

            DO 20 K = 0, NSTR - 1
               GL( K, LC ) = ( 2*K + 1 )*OPRIM( LC )*PMOM( K, LC )
   20       CONTINUE

            F  = 0.0


         ELSE
c                                    ** Do delta-M transformation

            F  = PMOM( NSTR, LC )
            OPRIM(LC) = SSALB(LC) * ( 1.- F ) / ( 1.- F * SSALB(LC) )
            DTAUCP( LC ) = ( 1.- F*SSALB( LC ) )*DTAUC( LC )
            TAUCPR( LC ) = TAUCPR( LC-1 ) + DTAUCP( LC )

            DO 30 K = 0, NSTR - 1
               GL( K, LC ) = ( 2*K + 1 ) * OPRIM( LC ) *
     &                       ( PMOM( K,LC ) - F ) / ( 1.- F )
   30       CONTINUE

         END IF

         FLYR( LC )   = F
         EXPBEA( LC ) = 0.0

         IF( FBEAM.GT.0.0 ) EXPBEA( LC ) = EXP( -TAUCPR( LC ) / UMU0 )

   40 CONTINUE
c                      ** If no thermal emission, cut off medium below
c                      ** absorption optical depth = ABSCUT ( note that
c                      ** delta-M transformation leaves absorption
c                      ** optical depth invariant ).  Not worth the
c                      ** trouble for one-layer problems, though.
      LYRCUT = .FALSE.

      IF( ABSTAU.GE.ABSCUT .AND. .NOT.PLANK .AND. IBCND.NE.1 .AND.
     &    NLYR.GT.1 ) LYRCUT = .TRUE.

      IF( .NOT.LYRCUT ) NCUT   = NLYR

c                             ** Set arrays defining location of user
c                             ** output levels within delta-M-scaled
c                             ** computational mesh
      DO 70 LU = 1, NTAU

         DO 50 LC = 1, NLYR

            IF( UTAU( LU ).GE.TAUC( LC - 1 ) .AND.
     &          UTAU( LU ).LE.TAUC( LC ) ) GO TO  60

   50    CONTINUE
         LC   = NLYR

   60    CONTINUE
         UTAUPR( LU ) = UTAU( LU )
         IF( DELTAM ) UTAUPR( LU ) = TAUCPR( LC - 1 ) +
     &                               ( 1.- SSALB( LC )*FLYR( LC ) )*
     &                               ( UTAU( LU ) - TAUC( LC-1 ) )
         LAYRU( LU ) = LC

   70 CONTINUE
c                      ** Calculate computational polar angle cosines
c                      ** and associated quadrature weights for Gaussian
c                      ** quadrature on the interval (0,1) (upward)
      NN   = NSTR / 2

      CALL QGAUSN( NN, CMU, CWT )
c                                  ** Downward (neg) angles and weights
      DO 80 IQ = 1, NN
         CMU( IQ + NN ) = - CMU( IQ )
         CWT( IQ + NN ) = CWT( IQ )
   80 CONTINUE


      IF( FBEAM.GT.0.0 ) THEN
c                               ** Compare beam angle to comput. angles
         DO 90 IQ = 1, NN

            IF( ABS( UMU0 - CMU( IQ ) ) / UMU0.LT.1.E-4 ) CALL ERRMSG(
     &          'SETDIS--beam angle=computational angle; change NSTR',
     &          .True. )

   90    CONTINUE

      END IF


      IF( .NOT.USRANG .OR. ( ONLYFL .AND. MAXUMU.GE.NSTR ) ) THEN

c                                   ** Set output polar angles to
c                                   ** computational polar angles
         NUMU   = NSTR

         DO 100 IU = 1, NN
            UMU( IU ) = - CMU( NN + 1 - IU )
  100    CONTINUE

         DO 110 IU = NN + 1, NSTR
            UMU( IU ) = CMU( IU - NN )
  110    CONTINUE

      END IF


      IF( USRANG .AND. IBCND.EQ.1 ) THEN

c                               ** Shift positive user angle cosines to
c                               ** upper locations and put negatives
c                               ** in lower locations
         DO 120 IU = 1, NUMU
            UMU( IU + NUMU ) = UMU( IU )
  120    CONTINUE

         DO 130 IU = 1, NUMU
            UMU( IU ) = -UMU( 2*NUMU + 1 - IU )
  130    CONTINUE

         NUMU   = 2*NUMU

      END IF


      IF( .NOT.LYRCUT .AND. .NOT.LAMBER ) THEN

         DO 140 K = 0, NSTR
            HLPR( K ) = ( 2*K + 1 )*HL( K )
  140    CONTINUE

      END IF


      RETURN
      END