SUBROUTINE ALBTRN( ALBEDO, AMB, APB, ARRAY, B, BDR, CBAND, CC, 1,15
     &                   CMU, CWT, DTAUCP, EVAL, EVECC, GL, GC, GU,
     &                   IPVT, KK, LL, NLYR, NN, NSTR, NUMU, PRNT,
     &                   TAUCPR, UMU, U0U, WK, YLMC, YLMU, Z, AAD,
     &                   EVALD, EVECCD, WKD, MI, MI9M2, MAXULV, MAXUMU,
     &                   MXCMU, MXUMU, NNLYRI, SQT, ALBMED, TRNMED )

c    DISORT special case to get only albedo and transmissivity
c    of entire medium as a function of incident beam angle
c    (many simplifications because boundary condition is just
c    isotropic illumination, there are no thermal sources, and
c    particular solutions do not need to be computed).  See
c    Ref. S2 and references therein for details.

c    The basic idea is as follows.  The reciprocity principle leads to
c    the following relationships for a plane-parallel, vertically
c    inhomogeneous medium lacking thermal (or other internal) sources:
c
c       albedo(theta) = u_0(theta) for unit-intensity isotropic
c                       illumination at *top* boundary
c
c       trans(theta) =  u_0(theta) for unit-intensity isotropic
c                       illumination at *bottom* boundary
c
c    where 
c
c       albedo(theta) = albedo for beam incidence at angle theta
c       trans(theta) = transmissivity for beam incidence at angle theta
c       u_0(theta) = upward azim-avg intensity at top boundary
c                    at angle theta


c   O U T P U T    V A R I A B L E S:
c
c       ALBMED(IU)   Albedo of the medium as a function of incident
c                    beam angle cosine UMU(IU)
c     
c       TRNMED(IU)   Transmissivity of the medium as a function of
c                    incident beam angle cosine UMU(IU)


c    I N T E R N A L   V A R I A B L E S:

c       NCD         number of diagonals below/above main diagonal

c       RCOND       estimate of the reciprocal condition of matrix
c                   CBAND; for system  CBAND*X = B, relative 
c                   perturbations in CBAND and B of size epsilon may
c                   cause relative perturbations in X of size 
c                   epsilon/RCOND.  If RCOND is so small that 
c                          1.0 + RCOND .EQ. 1.0
c                   is true, then CBAND may be singular to working
c                   precision.

c       CBAND       Left-hand side matrix of linear system Eq. SC(5),
c                   scaled by Eq. SC(12); in banded form required
c                   by LINPACK solution routines

c       NCOL        number of columns in CBAND matrix

c       IPVT        INTEGER vector of pivot indices

c       (most others documented in DISORT)

c   Called by- DISORT
c   Calls- LEPOLY, ZEROIT, SGBCO, SOLEIG, TERPEV, SETMTX, SOLVE1,
c          ALTRIN, SPALTR, PRALTR
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

      INTEGER   MAXULV, MAXUMU, MI, MI9M2, MXCMU, MXUMU, NLYR, NN,
     &          NNLYRI, NSTR, NUMU
      REAL      ALBEDO
c     ..
c     .. Array Arguments ..

      LOGICAL   PRNT( * )
      INTEGER   IPVT( * )
      REAL      ALBMED( MAXUMU ), AMB( MI, MI ), APB( MI, MI ),
     &          ARRAY( MXCMU, MXCMU ), B( NNLYRI ), BDR( MI, 0:MI ),
     &          CBAND( MI9M2, NNLYRI ), CC( MXCMU, MXCMU ),
     &          CMU( MXCMU ), CWT( MXCMU ), DTAUCP( * ), EVAL( MI ),
     &          EVECC( MXCMU, MXCMU ), GC( MXCMU, MXCMU, * ),
     &          GL( 0:MXCMU, * ), GU( MXUMU, MXCMU, * ), KK( MXCMU, * ),
     &          LL( MXCMU, * ), TAUCPR( 0:* ), TRNMED( MAXUMU ),
     &          U0U( MAXUMU, MAXULV ), UMU( MAXUMU ), WK( MXCMU ),
     &          YLMC( 0:MXCMU, MXCMU ), YLMU( 0:MXCMU, * ), Z( NNLYRI ),
     &          SQT( * )

      DOUBLE PRECISION AAD( MI, MI ), EVALD( MI ), EVECCD( MI, MI ),
     &                 WKD( MXCMU )
c     ..
c     .. Local Scalars ..

      LOGICAL   LAMBER, LYRCUT
      INTEGER   IQ, IU, L, LC, MAZIM, NCD, NCOL, NCUT
      REAL      DELM0, FISOT, RCOND, SGN, SPHALB, SPHTRN
c     ..
c     .. External Subroutines ..

      EXTERNAL  ALTRIN, LEPOLY, PRALTR, SETMTX, SGBCO, SOLEIG, SOLVE1,
     &          SPALTR, TERPEV, ZEROIT
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC EXP
c     ..

      MAZIM  = 0
      DELM0  = 1.0
c                    ** Set DISORT variables that are ignored in this
c                    ** special case but are needed below in argument
c                    ** lists of subroutines shared with general case
      NCUT   = NLYR
      LYRCUT = .FALSE.
      FISOT  = 1.0
      LAMBER = .TRUE.
c                          ** Get Legendre polynomials for computational
c                          ** and user polar angle cosines

      CALL LEPOLY( NUMU, MAZIM, MXCMU, NSTR-1, UMU, SQT, YLMU )

      CALL LEPOLY( NN, MAZIM, MXCMU, NSTR-1, CMU, SQT, YLMC )

c                       ** Evaluate Legendre polynomials with negative
c                       ** arguments from those with positive arguments;
c                       ** Dave/Armstrong Eq. (15)
      SGN  = - 1.0

      DO 20 L = MAZIM, NSTR - 1

         SGN  = - SGN

         DO 10 IQ = NN + 1, NSTR
            YLMC( L, IQ ) = SGN*YLMC( L, IQ - NN )
   10    CONTINUE

   20 CONTINUE
c                                  ** Zero out bottom reflectivity
c                                  ** (ALBEDO is used only in analytic
c                                  ** formulae involving ALBEDO = 0
c                                  ** solutions; Eqs 16-17 of Ref S2)

      CALL ZEROIT( BDR, MI*( MI + 1 ) )


c ===================  BEGIN LOOP ON COMPUTATIONAL LAYERS  =============

      DO 30 LC = 1, NLYR

c                        ** Solve eigenfunction problem in Eq. STWJ(8b)

         CALL SOLEIG( AMB, APB, ARRAY, CMU, CWT, GL( 0,LC ), MI, MAZIM,
     &                MXCMU, NN, NSTR, YLMC, CC, EVECC, EVAL,
     &                KK( 1,LC ), GC( 1,1,LC ), AAD, EVECCD, EVALD,
     &                WKD )

c                          ** Interpolate eigenvectors to user angles

         CALL TERPEV( CWT, EVECC, GL( 0,LC ), GU( 1,1,LC ), MAZIM,
     &                MXCMU, MXUMU, NN, NSTR, NUMU, WK, YLMC, YLMU )

   30 CONTINUE

c ===================  END LOOP ON COMPUTATIONAL LAYERS  ===============


c                      ** Set coefficient matrix (CBAND) of equations
c                      ** combining boundary and layer interface 
c                      ** conditions (in band-storage mode required by
c                      ** LINPACK routines)

      CALL SETMTX( BDR, CBAND, CMU, CWT, DELM0, DTAUCP, GC, KK,
     &             LAMBER, LYRCUT, MI, MI9M2, MXCMU, NCOL, NCUT,
     &             NNLYRI, NN, NSTR, TAUCPR, WK )

c                      ** LU-decompose the coeff. matrix (LINPACK)

      NCD  = 3*NN - 1
      CALL SGBCO( CBAND, MI9M2, NCOL, NCD, NCD, IPVT, RCOND, Z )
      IF( 1.0+RCOND .EQ. 1.0 )
     &    CALL ERRMSG('ALBTRN--SGBCO says matrix near singular',.FALSE.)

c                             ** First, illuminate from top; if only
c                             ** one layer, this will give us everything

c                             ** Solve for constants of integration in
c                             ** homogeneous solution

      CALL SOLVE1( B, CBAND, FISOT, 1, IPVT, LL, MI9M2, MXCMU,
     &             NCOL, NLYR, NN, NNLYRI, NSTR )

c                             ** Compute azimuthally-averaged intensity
c                             ** at user angles; gives albedo if multi-
c                             ** layer (Eq. 9 of Ref S2); gives both
c                             ** albedo and transmissivity if single
c                             ** layer (Eqs. 3-4 of Ref S2)

      CALL ALTRIN( GU, KK, LL, MXCMU, MXUMU, MAXUMU, NLYR, NN, NSTR,
     &             NUMU, TAUCPR, UMU, U0U, WK )

c                               ** Get beam-incidence albedos from
c                               ** reciprocity principle
      DO 40 IU = 1, NUMU / 2
         ALBMED( IU ) = U0U( IU + NUMU/2, 1 )
   40 CONTINUE


      IF( NLYR.EQ.1 ) THEN

         DO 50 IU = 1, NUMU / 2
c                               ** Get beam-incidence transmissivities
c                               ** from reciprocity principle (1 layer);
c                               ** flip them end over end to correspond
c                               ** to positive UMU instead of negative

            TRNMED( IU ) = U0U( NUMU/2 + 1 - IU, 2 )
     &                    + EXP( -TAUCPR( NLYR ) / UMU( IU + NUMU/2 ) )

   50    CONTINUE

      ELSE
c                             ** Second, illuminate from bottom
c                             ** (if multiple layers)

         CALL SOLVE1( B, CBAND, FISOT, 2, IPVT, LL, MI9M2, MXCMU,
     &                NCOL, NLYR, NN, NNLYRI, NSTR )

         CALL ALTRIN( GU, KK, LL, MXCMU, MXUMU, MAXUMU, NLYR, NN, NSTR,
     &                NUMU, TAUCPR, UMU, U0U, WK )

c                               ** Get beam-incidence transmissivities
c                               ** from reciprocity principle
         DO 60 IU = 1, NUMU / 2
            TRNMED(IU) = U0U( IU + NUMU/2, 1 )
     &                   + EXP( - TAUCPR(NLYR) / UMU(IU+NUMU/2) )
   60    CONTINUE

      END IF


      IF( ALBEDO.GT.0.0 ) THEN

c                             ** Get spherical albedo and transmissivity
         IF( NLYR.EQ.1 ) THEN

            CALL SPALTR( CMU, CWT, GC, KK, LL, MXCMU, NLYR,
     &                    NN, NSTR, TAUCPR, SPHALB, SPHTRN )
         ELSE

            CALL SPALTR( CMU, CWT, GC, KK, LL, MXCMU, NLYR,
     &                    NN, NSTR, TAUCPR, SPHTRN, SPHALB )
         END IF

c                                ** Ref. S2, Eqs. 16-17 (these eqs. have
c                                ** a simple physical interpretation
c                                ** like that of adding-doubling eqs.)
         DO 70 IU = 1, NUMU

            ALBMED(IU) = ALBMED(IU) + ( ALBEDO / (1.-ALBEDO*SPHALB) )
     &                                * SPHTRN * TRNMED(IU)

            TRNMED(IU) = TRNMED(IU) + ( ALBEDO / (1.-ALBEDO*SPHALB) )
     &                                * SPHALB * TRNMED(IU)
   70    CONTINUE

      END IF
c                          ** Return UMU to all positive values, to
c                          ** agree with ordering in ALBMED, TRNMED
      NUMU  = NUMU / 2
      DO 80 IU = 1, NUMU
         UMU( IU ) = UMU( IU + NUMU )
   80 CONTINUE

      IF( PRNT(6) ) CALL PRALTR( UMU, NUMU, ALBMED, TRNMED )

      RETURN
      END