SUBROUTINE USRINT( BPLANK, CMU, CWT, DELM0, DTAUCP, EMU, EXPBEA, 1
     &                   FBEAM, FISOT, GC, GU, KK, LAMBER, LAYRU, LL,
     &                   LYRCUT, MAZIM, MXCMU, MXULV, MXUMU, NCUT, NLYR,
     &                   NN, NSTR, PLANK, NUMU, NTAU, PI, RMU, TAUCPR,
     &                   TPLANK, UMU, UMU0, UTAUPR, WK, ZBEAM, Z0U, Z1U,
     &                   ZZ, ZPLK0, ZPLK1, UUM )

c       Computes intensity components at user output angles
c       for azimuthal expansion terms in Eq. SD(2)
c
c
c   I N P U T    V A R I A B L E S:
c
c       BPLANK :  Integrated Planck function for emission from
c                 bottom boundary
c
c       CMU    :  Abscissae for Gauss quadrature over angle cosine
c
c       CWT    :  Weights for Gauss quadrature over angle cosine
c
c       DELM0  :  Kronecker delta, delta-sub-M0
c
c       EMU    :  Surface directional emissivity (user angles)
c
c       EXPBEA :  Transmission of incident beam, EXP(-TAUCPR/UMU0)
c
c       GC     :  Eigenvectors at polar quadrature angles, SC(1)
c
c       GU     :  Eigenvectors interpolated to user polar angles
c                    (i.e., G in Eq. SC(1) )
c
c       KK     :  Eigenvalues of coeff. matrix in Eq. SS(7)
c
c       LAYRU  :  Layer number of user level UTAU
c
c       LL     :  Constants of integration in Eq. SC(1), obtained
c                 by solving scaled version of Eq. SC(5);
c                 exponential term of Eq. SC(12) not included
c
c       LYRCUT :  Logical flag for truncation of computational layer
c
c       MAZIM  :  Order of azimuthal component
c
c       NCUT   :  Total number of computational layers considered
c
c       NN     :  Order of double-Gauss quadrature (NSTR/2)
c
c       RMU    :  Surface bidirectional reflectivity (user angles)
c
c       TAUCPR :  Cumulative optical depth (delta-M-Scaled)
c
c       TPLANK :  Integrated Planck function for emission from
c                 top boundary
c
c       UTAUPR :  Optical depths of user output levels in delta-M
c                 coordinates;  equal to UTAU if no delta-M
c
c       Z0U    :  Z-sub-zero in Eq. SS(16) interpolated to user
c                 angles from an equation derived from SS(16)
c
c       Z1U    :  Z-sub-one in Eq. SS(16) interpolated to user
c                 angles from an equation derived from SS(16)
c
c       ZZ     :  Beam source vectors in Eq. SS(19)
c
c       ZPLK0  :  Thermal source vectors Z0, by solving Eq. SS(16)
c
c       ZPLK1  :  Thermal source vectors Z1, by solving Eq. SS(16)
c
c       ZBEAM  :  Incident-beam source vectors
c
c       nobeam :  1 => exclude atmospheric single scattering contribution
c                      to radiance. (direct beam contribution to surface
c                      reflectance is still included).
c                 0 => normal disort treatment of direct beam contribution
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       UUM    :  Azimuthal components of the intensity in EQ. STWJ(5)
c
c
c    I N T E R N A L    V A R I A B L E S:
c
c       BNDDIR :  Direct intensity down at the bottom boundary
c       BNDDFU :  Diffuse intensity down at the bottom boundary
c       BNDINT :  Intensity attenuated at both boundaries, STWJ(25-6)
c       DTAU   :  Optical depth of a computational layer
c       LYREND :  End layer of integration
c       LYRSTR :  Start layer of integration
c       PALINT :  Intensity component from parallel beam
c       PLKINT :  Intensity component from planck source
c       WK     :  Scratch vector for saving EXP evaluations
c
c       All the exponential factors ( EXP1, EXPN,... etc.)
c       come from the substitution of constants of integration in
c       Eq. SC(12) into Eqs. S1(8-9).  They all have negative
c       arguments so there should never be overflow problems.
c
c   Called by- DISORT
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

      LOGICAL   LAMBER, LYRCUT, PLANK
      INTEGER   MAZIM, MXCMU, MXULV, MXUMU, NCUT, NLYR, NN, NSTR, NTAU,
     &          NUMU
      REAL      BPLANK, DELM0, FBEAM, FISOT, PI, TPLANK, UMU0
c     ..
c     .. Array Arguments ..

      INTEGER   LAYRU( * )
      REAL      CMU( MXCMU ), CWT( MXCMU ), DTAUCP( * ), EMU( MXUMU ),
     &          EXPBEA( 0:* ), GC( MXCMU, MXCMU, * ),
     &          GU( MXUMU, MXCMU, * ), KK( MXCMU, * ), LL( MXCMU, * ),
     &          RMU( MXUMU, 0:* ), TAUCPR( 0:* ), UMU( * ),
     &          UTAUPR( MXULV ), UUM( MXUMU, MXULV ), WK( MXCMU ),
     &          Z0U( MXUMU, * ), Z1U( MXUMU, * ), ZBEAM( MXUMU, * ),
     &          ZPLK0( MXCMU, * ), ZPLK1( MXCMU, * ), ZZ( MXCMU, * )
c     ..
c     .. Local Scalars ..

      LOGICAL   NEGUMU

c pjr-- beamon --   New logical variable indicates when 
c                   to include parallel beam contribution

      INTEGER   IQ, IU, JQ, LC, LU, LYREND, LYRSTR, LYU
      REAL      BNDDFU, BNDDIR, BNDINT, DENOM, DFUINT, DTAU, DTAU1,
     &          DTAU2, EXP0, EXP1, EXP2, EXPN, FACT, PALINT, PLKINT, SGN
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, EXP
c     ..

c                          ** Incorporate constants of integration into
c                          ** interpolated eigenvectors



      DO 30 LC = 1, NCUT

         DO 20 IQ = 1, NSTR

            DO 10 IU = 1, NUMU
               GU( IU, IQ, LC ) = GU( IU, IQ, LC )*LL( IQ, LC )
   10       CONTINUE

   20    CONTINUE

   30 CONTINUE
c                           ** Loop over levels at which intensities
c                           ** are desired ('user output levels')
      DO 160 LU = 1, NTAU

         IF( fbeam.gt.0.0 ) EXP0 = EXP( - UTAUPR(LU) / UMU0 )
         LYU  = LAYRU( LU )
c                              ** Loop over polar angles at which
c                              ** intensities are desired
         DO 150 IU = 1, NUMU

            IF( LYRCUT .AND. LYU.GT.NCUT ) GO TO  150

            NEGUMU = UMU( IU ).LT.0.0

            IF( NEGUMU ) THEN

               LYRSTR = 1
               LYREND = LYU - 1
               SGN    = -1.0

            ELSE

               LYRSTR = LYU + 1
               LYREND = NCUT
               SGN    = 1.0

            END IF
c                          ** For downward intensity, integrate from top
c                          ** to LYU-1 in Eq. S1(8); for upward,
c                          ** integrate from bottom to LYU+1 in S1(9)
            PALINT = 0.0
            PLKINT = 0.0

            DO 60 LC = LYRSTR, LYREND

               DTAU = DTAUCP( LC )
               EXP1 = EXP(               EXP2 = EXP(
               IF ( PLANK .AND. MAZIM.EQ.0 )
     &           PLKINT = PLKINT + SGN * ( Z0U(IU,LC) * (EXP1 - EXP2) +
     &                    Z1U(IU,LC) * (     &                                   (TAUCPR(LC) + UMU(IU))*EXP2 ) )

               IF( fbeam.gt.0.0 ) THEN

                  DENOM  = 1.+ UMU( IU ) / UMU0

                  IF( ABS( DENOM ).LT.0.0001 ) THEN
c                                                   ** L'Hospital limit
                     EXPN  = ( DTAU / UMU0 )*EXP0

                  ELSE

                     EXPN  = ( EXP1*EXPBEA( LC-1 ) -
     &                         EXP2*EXPBEA( LC ) )*SGN / DENOM
                  END IF

                  PALINT = PALINT + ZBEAM( IU, LC )*EXPN

               END IF

c                                                   ** KK is negative
               DO 40 IQ = 1, NN

                  WK( IQ ) = EXP( KK( IQ,LC )*DTAU )
                  DENOM  = 1.0 + UMU( IU )*KK( IQ, LC )

                  IF( ABS( DENOM ).LT.0.0001 ) THEN
c                                                   ** L'Hospital limit
                     EXPN  = DTAU / UMU( IU )*EXP2

                  ELSE

                     EXPN  = SGN*( EXP1*WK( IQ ) - EXP2 ) / DENOM

                  END IF

                  PALINT = PALINT + GU( IU, IQ, LC )*EXPN

   40          CONTINUE

c                                                   ** KK is positive
               DO 50 IQ = NN + 1, NSTR

                  DENOM  = 1.0 + UMU( IU )*KK( IQ, LC )

                  IF( ABS( DENOM ).LT.0.0001 ) THEN
c                                                   ** L'Hospital limit
                     EXPN  = - DTAU / UMU( IU )*EXP1

                  ELSE

                     EXPN  = SGN*( EXP1 - EXP2*WK(NSTR+1 - IQ) ) /DENOM

                  END IF

                  PALINT = PALINT + GU( IU, IQ, LC )*EXPN

   50          CONTINUE


   60       CONTINUE
c                           ** Calculate contribution from user
c                           ** output level to next computational level

            DTAU1  = UTAUPR( LU ) - TAUCPR( LYU - 1 )
            DTAU2  = UTAUPR( LU ) - TAUCPR( LYU )

            IF( ABS( DTAU1 ).LT.1.E-6 .AND. NEGUMU ) GO TO  90
            IF( ABS( DTAU2 ).LT.1.E-6 .AND. (.NOT.NEGUMU) ) GO TO  90

            IF( NEGUMU )      EXP1   = EXP( DTAU1 / UMU( IU ) )
            IF( .NOT.NEGUMU ) EXP2   = EXP( DTAU2 / UMU( IU ) )

            IF( fbeam.gt.0.0 ) THEN

               DENOM  = 1.+ UMU( IU ) / UMU0

               IF( ABS( DENOM ).LT.0.0001 ) THEN

                  EXPN  = ( DTAU1 / UMU0 )*EXP0

               ELSE IF( NEGUMU ) THEN

                  EXPN  = ( EXP0 - EXPBEA( LYU - 1 )*EXP1 ) / DENOM

               ELSE

                  EXPN  = ( EXP0 - EXPBEA( LYU )*EXP2 ) / DENOM

               END IF

               PALINT = PALINT + ZBEAM( IU, LYU )*EXPN

            END IF

c                                                   ** KK is negative
            DTAU  = DTAUCP( LYU )

            DO 70 IQ = 1, NN

               DENOM  = 1.+ UMU( IU )*KK( IQ, LYU )

               IF( ABS( DENOM ).LT.0.0001 ) THEN

                  EXPN   = -DTAU2 / UMU( IU )*EXP2

               ELSE IF( NEGUMU ) THEN

                  EXPN   = ( EXP( -KK( IQ,LYU )*DTAU2 ) -
     &                       EXP(  KK( IQ,LYU )*DTAU  )*EXP1 ) / DENOM

               ELSE

                  EXPN   = ( EXP( -KK( IQ,LYU )*DTAU2 ) - EXP2 ) / DENOM

               END IF

               PALINT = PALINT + GU( IU, IQ, LYU )*EXPN

   70       CONTINUE

c                                                   ** KK is positive
            DO 80 IQ = NN + 1, NSTR

               DENOM  = 1.+ UMU( IU )*KK( IQ, LYU )

               IF( ABS( DENOM ).LT.0.0001 ) THEN

                  EXPN  = -DTAU1 / UMU( IU )*EXP1

               ELSE IF( NEGUMU ) THEN

                  EXPN  = ( EXP( -KK( IQ,LYU )*DTAU1 ) - EXP1 ) / DENOM

               ELSE

                  EXPN  = ( EXP( -KK( IQ,LYU )*DTAU1 ) -
     &                      EXP( -KK( IQ,LYU )*DTAU  )*EXP2 ) / DENOM
               END IF

               PALINT = PALINT + GU( IU, IQ, LYU )*EXPN

   80       CONTINUE


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

               IF( NEGUMU ) THEN

                  EXPN  = EXP1
                  FACT  = TAUCPR( LYU - 1 ) + UMU( IU )

               ELSE

                  EXPN  = EXP2
                  FACT  = TAUCPR( LYU ) + UMU( IU )

               END IF

               PLKINT = PLKINT + Z0U( IU, LYU )*( 1.- EXPN ) +
     &                  Z1U( IU, LYU )*( UTAUPR( LU ) + UMU( IU )
     &                                   - FACT*EXPN )
            END IF

c                            ** Calculate intensity components
c                            ** attenuated at both boundaries.
c                            ** NOTE: no azimuthal intensity
c                            ** component for isotropic surface
   90       CONTINUE
            BNDINT = 0.0

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

               BNDINT = ( FISOT + TPLANK )*
     &                  EXP( UTAUPR( LU ) / UMU( IU ) )


            ELSE IF( .NOT.NEGUMU ) THEN

               IF( LYRCUT .OR. ( LAMBER .AND. MAZIM.GT.0 ) ) GO TO 140

               DO 100 JQ = NN + 1, NSTR
                  WK( JQ ) = EXP( - KK( JQ,NLYR )*DTAUCP( NLYR ) )
  100          CONTINUE

               BNDDFU = 0.0

               DO 130 IQ = NN, 1, -1

                  DFUINT = 0.0
                  DO 110 JQ = 1, NN
                     DFUINT = DFUINT + GC( IQ, JQ, NLYR )*LL( JQ, NLYR )
  110             CONTINUE

                  DO 120 JQ = NN + 1, NSTR
                     DFUINT = DFUINT + GC( IQ, JQ, NLYR )*
     &                                 LL( JQ, NLYR )*WK( JQ )
  120             CONTINUE

                  IF( fbeam.gt.0.0 ) DFUINT = DFUINT +
     &                               ZZ( IQ, NLYR )*EXPBEA( NLYR )

                  DFUINT = DFUINT + DELM0*( ZPLK0( IQ,NLYR ) +
     &                              ZPLK1( IQ,NLYR )*TAUCPR( NLYR ) )
                 BNDDFU = BNDDFU + ( 1.+ DELM0 ) * RMU(IU,NN+1-IQ)
     &                           * CMU(NN+1-IQ) * CWT(NN+1-IQ) * DFUINT
  130          CONTINUE

               BNDDIR = 0.0
               IF( FBEAM.GT.0.0 ) BNDDIR = UMU0*FBEAM / PI*RMU( IU, 0 )*
     &                                     EXPBEA( NLYR )

               BNDINT = ( BNDDFU + BNDDIR + DELM0 * EMU(IU) * BPLANK )
     &                  * EXP(            END IF

  140       CONTINUE

            UUM( IU, LU ) = PALINT + PLKINT + BNDINT

  150    CONTINUE

  160 CONTINUE


      RETURN
      END