SUBROUTINE UPISOT( ARRAY, CC, CMU, IPVT, MXCMU, NN, NSTR, OPRIM, 1,3
     &                   WK, XR0, XR1, Z0, Z1, ZPLK0, ZPLK1 )

c       Finds the particular solution of thermal radiation of SS(15)
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       OPRIM  :  Delta-M scaled single scattering albedo
c
c       XR0    :  Expansion of thermal source function
c
c       XR1    :  Expansion of thermal source function Eqs. SS(14-16)
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       Z0     :  Solution vectors Z-sub-zero of Eq. SS(16)
c
c       Z1     :  Solution vectors Z-sub-one  of Eq. SS(16)
c
c       ZPLK0, :  Permanent storage for Z0,Z1, but re-ordered
c        ZPLK1
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(16)
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   MXCMU, NN, NSTR
      REAL      OPRIM, XR0, XR1
c     ..
c     .. Array Arguments ..

      INTEGER   IPVT( * )
      REAL      ARRAY( MXCMU, MXCMU ), CC( MXCMU, MXCMU ), CMU( MXCMU ),
     &          WK( MXCMU ), Z0( MXCMU ), Z1( MXCMU ), ZPLK0( MXCMU ),
     &          ZPLK1( MXCMU )
c     ..
c     .. Local Scalars ..

      INTEGER   IQ, JQ
      REAL      RCOND
c     ..
c     .. External Subroutines ..

      EXTERNAL  ERRMSG, SGECO, SGESL
c     ..


      DO 20 IQ = 1, NSTR

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

         ARRAY( IQ, IQ ) = 1.0 + ARRAY( IQ, IQ )

         Z1( IQ ) = XR1
         Z0( IQ ) = ( 1.- OPRIM )*XR0 + CMU( IQ )*Z1( IQ )

   20 CONTINUE
c                       ** Solve linear equations: same as in UPBEAM,
c                       ** except ZJ replaced by Z0
      RCOND  = 0.0

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

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

      CALL SGESL( ARRAY, MXCMU, NSTR, IPVT, Z0, 0 )

CDIR$ IVDEP
      DO 30 IQ = 1, NN
         ZPLK0( IQ + NN ) = Z0( IQ )
         ZPLK1( IQ + NN ) = Z1( IQ )
         ZPLK0( NN + 1 - IQ ) = Z0( IQ + NN )
         ZPLK1( NN + 1 - IQ ) = Z1( IQ + NN )
   30 CONTINUE


      RETURN
      END