SUBROUTINE UPBEAM( ARRAY, CC, CMU, DELM0, FBEAM, GL, IPVT, MAZIM, 1,3
     &                   MXCMU, NN, NSTR, PI, UMU0, WK, YLM0, YLMC, ZJ,
     &                   ZZ )

c         Finds the incident-beam particular solution of SS(18)
c
c
c   I N P U T    V A R I A B L E S:
c
c       CC     :  C-sub-ij in Eq. SS(5)
c
c       CMU    :  Abscissae for Gauss quadrature over angle cosine
c
c       DELM0  :  Kronecker delta, delta-sub-m0
c
c       GL     :  Delta-M scaled Legendre coefficients of phase function
c                 (including factors 2L+1 and single-scatter albedo)
c
c       MAZIM  :  Order of azimuthal component
c
c       YLM0   :  Normalized associated Legendre polynomial
c                 at the beam angle
c
c       YLMC   :  Normalized associated Legendre polynomial
c                 at the quadrature angles
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       ZJ     :  Right-hand side vector X-sub-zero in SS(19); also the
c                 solution vector Z-sub-zero after solving that system
c
c       ZZ     :  Permanent storage for ZJ, but re-ordered
c
c
c   I N T E R N A L    V A R I A B L E S:
c
c       ARRAY  :  Coefficient matrix in left-hand side of Eq. SS(19)
c       IPVT   :  Integer vector of pivot indices required by LINPACK
c       WK     :  Scratch array required by LINPACK
c
c   Called by- DISORT
c   Calls- SGECO, ERRMSG, SGESL
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

      INTEGER   MAZIM, MXCMU, NN, NSTR
      REAL      DELM0, FBEAM, PI, UMU0
c     ..
c     .. Array Arguments ..

      INTEGER   IPVT( * )
      REAL      ARRAY( MXCMU, MXCMU ), CC( MXCMU, MXCMU ), CMU( MXCMU ),
     &          GL( 0:MXCMU ), WK( MXCMU ), YLM0( 0:MXCMU ),
     &          YLMC( 0:MXCMU, * ), ZJ( MXCMU ), ZZ( MXCMU )
c     ..
c     .. Local Scalars ..

      INTEGER   IQ, JOB, JQ, K
      REAL      RCOND, SUM
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, SGECO, SGESL
c     ..


      DO 30 IQ = 1, NSTR

         DO 10 JQ = 1, NSTR
            ARRAY( IQ, JQ ) = -CC( IQ, JQ )
   10    CONTINUE

         ARRAY( IQ, IQ ) = 1.+ CMU( IQ ) / UMU0 + ARRAY( IQ, IQ )

         SUM  = 0.
         DO 20 K = MAZIM, NSTR - 1
            SUM  = SUM + GL( K )*YLMC( K, IQ )*YLM0( K )
   20    CONTINUE

         ZJ( IQ ) = ( 2.- DELM0 )*FBEAM*SUM / ( 4.*PI )
   30 CONTINUE

c                  ** Find L-U (lower/upper triangular) decomposition
c                  ** of ARRAY and see if it is nearly singular
c                  ** (NOTE:  ARRAY is altered)
      RCOND  = 0.0

      CALL SGECO( ARRAY, MXCMU, NSTR, IPVT, RCOND, WK )

      IF( 1.0 + RCOND.EQ.1.0 )
     &    CALL ERRMSG('UPBEAM--SGECO says matrix near singular',.FALSE.)

c                ** Solve linear system with coeff matrix ARRAY
c                ** (assumed already L-U decomposed) and R.H. side(s)
c                ** ZJ;  return solution(s) in ZJ
      JOB  = 0

      CALL SGESL( ARRAY, MXCMU, NSTR, IPVT, ZJ, JOB )

CDIR$ IVDEP
      DO 40 IQ = 1, NN
         ZZ( IQ + NN )     = ZJ( IQ )
         ZZ( NN + 1 - IQ ) = ZJ( IQ + NN )
   40 CONTINUE


      RETURN
      END