SUBROUTINE TERPSO( CWT, DELM0, FBEAM, GL, MAZIM, MXCMU, PLANK, 1
     &                   NUMU, NSTR, OPRIM, PI, YLM0, YLMC, YLMU, PSI,
     &                   XR0, XR1, Z0, ZJ, ZBEAM, Z0U, Z1U )

c         Interpolates source functions to user angles
c
c
c    I N P U T      V A R I A B L E S:
c
c       CWT    :  Weights 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       OPRIM  :  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       YLM0   :  Normalized associated Legendre polynomial
c                 at the beam angle
c
c       YLMC   :  Normalized associated Legendre polynomial
c                 at the quadrature angles
c
c       YLMU   :  Normalized associated Legendre polynomial
c                 at the user angles
c
c       Z0     :  Solution vectors Z-sub-zero of Eq. SS(16)
c
c       ZJ     :  Solution vector Z-sub-zero after solving Eq. SS(19)
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       ZBEAM  :  Incident-beam source function at user angles
c
c       Z0U,Z1U:  Components of a linear-in-optical-depth-dependent
c                 source (approximating the Planck emission source)
c
c
c   I N T E R N A L    V A R I A B L E S:
c
c       PSI    :  Sum just after square bracket in  Eq. SD(9)
c
c   Called by- DISORT
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

      LOGICAL   PLANK
      INTEGER   MAZIM, MXCMU, NSTR, NUMU
      REAL      DELM0, FBEAM, OPRIM, PI, XR0, XR1
c     ..
c     .. Array Arguments ..

      REAL      CWT( MXCMU ), GL( 0:MXCMU ), PSI( MXCMU ),
     &          YLM0( 0:MXCMU ), YLMC( 0:MXCMU, MXCMU ),
     &          YLMU( 0:MXCMU, * ), Z0( MXCMU ), Z0U( * ), Z1U( * ),
     &          ZBEAM( * ), ZJ( MXCMU )
c     ..
c     .. Local Scalars ..

      INTEGER   IQ, IU, JQ
      REAL      FACT, PSUM, SUM
c     ..


      IF( FBEAM.GT.0.0 ) THEN
c                                  ** Beam source terms; Eq. SD(9)

         DO 20 IQ = MAZIM, NSTR - 1

            PSUM   = 0.
            DO 10 JQ = 1, NSTR
               PSUM  = PSUM + CWT( JQ )*YLMC( IQ, JQ )*ZJ( JQ )
   10       CONTINUE

            PSI( IQ + 1 ) = 0.5*GL( IQ )*PSUM

   20    CONTINUE

         FACT   = ( 2. - DELM0 )*FBEAM / ( 4.0*PI )

         DO 40 IU = 1, NUMU

            SUM  = 0.
            DO 30 IQ = MAZIM, NSTR - 1
               SUM  = SUM + YLMU( IQ, IU )*
     &                     ( PSI( IQ + 1 ) + FACT*GL( IQ )*YLM0( IQ ) )
   30       CONTINUE

            ZBEAM( IU ) = SUM

   40    CONTINUE

      END IF


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

c                                   ** Thermal source terms, STWJ(27c)
         DO 60 IQ = MAZIM, NSTR - 1

            PSUM   = 0.0
            DO 50 JQ = 1, NSTR
               PSUM  = PSUM + CWT( JQ )*YLMC( IQ, JQ )*Z0( JQ )
   50       CONTINUE

            PSI( IQ + 1 ) = 0.5*GL( IQ )*PSUM

   60    CONTINUE

         DO 80 IU = 1, NUMU

            SUM  = 0.0
            DO 70 IQ = MAZIM, NSTR - 1
               SUM  = SUM + YLMU( IQ, IU )*PSI( IQ + 1 )
   70       CONTINUE

            Z0U( IU ) = SUM + ( 1.- OPRIM )*XR0
            Z1U( IU ) = XR1

   80    CONTINUE

      END IF


      RETURN
      END