c=======================================================================






c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c RCS version control information:
c $Header: /home/paul/rt/sbdart/RCS/sbdart.f,v 1.4 2000/07/13 23:43:37 paul Exp $
c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

c  Note: CDIR$ IVDEP comment lines are specific to Cray computers.  
c        They cause better optimization of loops immediately 
c        following, by ignoring potential dependencies.



      SUBROUTINE DISORT( NLYR, DTAUC, SSALB, PMOM, TEMPER, WVNMLO, 1,30
     &                   WVNMHI, USRTAU, NTAU, UTAU, NSTR, USRANG, NUMU,
     &                   UMU, NPHI, PHI, IBCND, FBEAM, UMU0, PHI0,
     &                   FISOT, LAMBER, ALBEDO, HL, BTEMP, TTEMP, TEMIS,
     &                   DELTAM, PLANK, ONLYFL, ACCUR, PRNT, HEADER,
     &                   MAXCLY, MAXULV, MAXUMU, MAXCMU, MAXPHI, RFLDIR,
     &                   RFLDN, FLUP, DFDT, UAVG, UU, U0U, ALBMED,
     &                   TRNMED)

c *******************************************************************
c       Plane-parallel discrete ordinates radiative transfer program
c             ( see DISORT.doc for complete documentation )
c *******************************************************************
c
c +------------------------------------------------------------------+
c  Calling Tree (omitting calls to ERRMSG):
c  (routines in parentheses are not in this file)
c
c  DISORT-+-(R1MACH)
c         +-SLFTST-+-(TSTBAD)
c         +-ZEROIT
c         +-CHEKIN-+-(WRTBAD)
c         |        +-(WRTDIM)
c         |        +-DREF
c         +-ZEROAL
c         +-SETDIS-+-QGAUSN (1)-+-(D1MACH)
c         +-PRTINP
c         +-ALBTRN-+-LEPOLY (2)
c         |        +-ZEROIT
c         |        +-SOLEIG (3)-+-ASYMTX-+-(D1MACH)
c         |        +-TERPEV
c         |        +-SETMTX (4)--ZEROIT
c         |        +-(SGBCO)
c         |        +-SOLVE1-+-ZEROIT
c         |        |        +-(SGBSL)
c         |        +-ALTRIN
c         |        +-SPALTR
c         |        +-PRALTR
c         +-PLKAVG-+-(R1MACH)
c         +-LEPOLY see 2
c         +-SURFAC-+-QGAUSN see 1
c         |        +-LEPOLY see 2
c         |        +-ZEROIT
c         +-SOLEIG see 3
c         +-UPBEAM-+-(SGECO)
c         |        +-(SGESL)
c         +-UPISOT-+-(SGECO)
c         |        +-(SGESL)
c         +-TERPEV
c         +-TERPSO
c         +-SETMTX see 4
c         +-SOLVE0-+-ZEROIT
c         |        +-(SGBCO)
c         |        +-(SGBSL)
c         +-FLUXES--ZEROIT
c         +-USRINT
c         +-CMPINT
c         +-PRAVIN
c         +-RATIO--(R1MACH)
c         +-PRTINT
c
c *** Intrinsic Functions used in DISORT package which take
c     non-negligible amount of time:
c
c    EXP :  Called by- ALBTRN, ALTRIN, CMPINT, FLUXES, SETDIS,
c                      SETMTX, SPALTR, USRINT, PLKAVG
c
c    SQRT : Called by- ASYMTX, SOLEIG
c
c +-------------------------------------------------------------------+
c
c  Index conventions (for all DO-loops and all variable descriptions):
c
c     IU     :  for user polar angles
c
c  IQ,JQ,KQ  :  for computational polar angles ('quadrature angles')
c
c   IQ/2     :  for half the computational polar angles (just the ones
c               in either 0-90 degrees, or 90-180 degrees)
c
c     J      :  for user azimuthal angles
c
c     K,L    :  for Legendre expansion coefficients or, alternatively,
c               subscripts of associated Legendre polynomials
c
c     LU     :  for user levels
c
c     LC     :  for computational layers (each having a different
c               single-scatter albedo and/or phase function)
c
c    LEV     :  for computational levels
c
c    MAZIM   :  for azimuthal components in Fourier cosine expansion
c               of intensity and phase function
c
c +------------------------------------------------------------------+
c
c               I N T E R N A L    V A R I A B L E S
c
c   AMB(IQ/2,IQ/2)    First matrix factor in reduced eigenvalue problem
c                     of Eqs. SS(12), STWJ(8E)  (used only in SOLEIG)
c
c   APB(IQ/2,IQ/2)    Second matrix factor in reduced eigenvalue problem
c                     of Eqs. SS(12), STWJ(8E)  (used only in SOLEIG)
c
c   ARRAY(IQ,IQ)      Scratch matrix for SOLEIG, UPBEAM and UPISOT
c                     (see each subroutine for definition)
c
c   B()               Right-hand side vector of Eq. SC(5) going into
c                     SOLVE0,1;  returns as solution vector
c                     vector  L, the constants of integration
c
c   BDR(IQ/2,0:IQ/2)  Bottom-boundary bidirectional reflectivity for a
c                     given azimuthal component.  First index always
c                     refers to a computational angle.  Second index:
c                     if zero, refers to incident beam angle UMU0;
c                     if non-zero, refers to a computational angle.
c
c   BEM(IQ/2)         Bottom-boundary directional emissivity at compu-
c                     tational angles.
c
c   BPLANK            Intensity emitted from bottom boundary
c
c   CBAND()           Matrix of left-hand side of the linear system
c                     Eq. SC(5), scaled by Eq. SC(12);  in banded
c                     form required by LINPACK solution routines
c
c   CC(IQ,IQ)         C-sub-IJ in Eq. SS(5)
c
c   CMU(IQ)           Computational polar angles (Gaussian)
c
c   CWT(IQ)           Quadrature weights corresponding to CMU
c
c   DELM0             Kronecker delta, delta-sub-M0, where M = MAZIM
c                     is the number of the Fourier component in the
c                     azimuth cosine expansion
c
c   DITHER            Small quantity subtracted from single-scattering
c                     albedos of unity, in order to avoid using special
c                     case formulas;  prevents an eigenvalue of exactly
c                     zero from occurring, which would cause an
c                     immediate overflow
c
c   DTAUCP(LC)        Computational-layer optical depths (delta-M-scaled
c                     if DELTAM = TRUE, otherwise equal to DTAUC)
c
c   EMU(IU)           Bottom-boundary directional emissivity at user
c                     angles.
c
c   EVAL(IQ)          Temporary storage for eigenvalues of Eq. SS(12)
c
c   EVECC(IQ,IQ)      Complete eigenvectors of SS(7) on return from
c                     SOLEIG; stored permanently in  GC
c
c   EXPBEA(LC)        Transmission of direct beam in delta-M optical
c                     depth coordinates
c
c   FLYR(LC)          Truncated fraction in delta-M method
c
c   GL(K,LC)          Phase function Legendre polynomial expansion
c                     coefficients, calculated from PMOM by
c                     including single-scattering albedo, factor
c                     2K+1, and (if DELTAM=TRUE) the delta-M
c                     scaling
c
c   GC(IQ,IQ,LC)      Eigenvectors at polar quadrature angles,
c                     g  in Eq. SC(1)
c
c   GU(IU,IQ,LC)      Eigenvectors interpolated to user polar angles
c                     ( g  in Eqs. SC(3) and S1(8-9), i.e.
c                       G without the L factor )
c
c   HLPR()            Legendre coefficients of bottom bidirectional
c                     reflectivity (after inclusion of 2K+1 factor)
c
c   IPVT(LC*IQ)       Integer vector of pivot indices for LINPACK
c                     routines
c
c   KK(IQ,LC)         Eigenvalues of coeff. matrix in Eq. SS(7)
c
c   KCONV             Counter in azimuth convergence test
c
c   LAYRU(LU)         Computational layer in which user output level
c                     UTAU(LU) is located
c
c   LL(IQ,LC)         Constants of integration L in Eq. SC(1),
c                     obtained by solving scaled version of Eq. SC(5)
c
c   LYRCUT            TRUE, radiation is assumed zero below layer
c                     NCUT because of almost complete absorption
c
c   NAZ               Number of azimuthal components considered
c
c   NCUT              Computational layer number in which absorption
c                     optical depth first exceeds ABSCUT
c
c   OPRIM(LC)         Single scattering albedo after delta-M scaling
c
c   PASS1             TRUE on first entry, FALSE thereafter
c
c   PKAG(0:LC)        Integrated Planck function for internal emission
c
c   PSI(IQ)           Sum just after square bracket in  Eq. SD(9)
c
c   RMU(IU,0:IQ)      Bottom-boundary bidirectional reflectivity for a
c                     given azimuthal component.  First index always
c                     refers to a user angle.  Second index:
c                     if zero, refers to incident beam angle UMU0;
c                     if non-zero, refers to a computational angle.
c
c   SQT(k)            Square root of k (used only in LEPOLY for
c                     computing associated Legendre polynomials)
c
c   TAUC(0:LC)        Cumulative optical depth (un-delta-M-scaled)
c
c   TAUCPR(0:LC)      Cumulative optical depth (delta-M-scaled if
c                     DELTAM = TRUE, otherwise equal to TAUC)
c
c   TPLANK            Intensity emitted from top boundary
c
c   UUM(IU,LU)        Expansion coefficients when the intensity
c                     (u-super-M) is expanded in Fourier cosine series
c                     in azimuth angle
c
c   U0C(IQ,LU)        Azimuthally-averaged intensity
c
c   UTAUPR(LU)        Optical depths of user output levels in delta-M
c                     coordinates;  equal to  UTAU(LU) if no delta-M
c
c   WK()              scratch array
c
c   XR0(LC)           X-sub-zero in expansion of thermal source func-
c                     tion preceding Eq. SS(14) (has no mu-dependence)
c
c   XR1(LC)           X-sub-one in expansion of thermal source func-
c                     tion;  see  Eqs. SS(14-16)
c
c   YLM0(L)           Normalized associated Legendre polynomial
c                     of subscript L at the beam angle (not saved
c                     as function of superscipt M)
c
c   YLMC(L,IQ)        Normalized associated Legendre polynomial
c                     of subscript L at the computational angles
c                     (not saved as function of superscipt M)
c
c   YLMU(L,IU)        Normalized associated Legendre polynomial
c                     of subscript L at the user angles
c                     (not saved as function of superscipt M)
c
c   Z()               scratch array used in SOLVE0, ALBTRN to solve
c                     a linear system for the constants of integration
c
c   Z0(IQ)            Solution vectors Z-sub-zero of Eq. SS(16)
c
c   Z0U(IU,LC)        Z-sub-zero in Eq. SS(16) interpolated to user
c                     angles from an equation derived from SS(16)
c
c   Z1(IQ)            Solution vectors Z-sub-one  of Eq. SS(16)
c
c   Z1U(IU,LC)        Z-sub-one in Eq. SS(16) interpolated to user
c                     angles from an equation derived from SS(16)
c
c   ZBEAM(IU,LC)      Particular solution for beam source
c
c   ZJ(IQ)            Right-hand side vector  X-sub-zero in
c                     Eq. SS(19), also the solution vector
c                     Z-sub-zero after solving that system
c
c   ZZ(IQ,LC)         Permanent storage for the beam source vectors ZJ
c
c   ZPLK0(IQ,LC)      Permanent storage for the thermal source
c                     vectors  Z0  obtained by solving  Eq. SS(16)
c
c   ZPLK1(IQ,LC)      Permanent storage for the thermal source
c                     vectors  Z1  obtained by solving  Eq. SS(16)
c
c +-------------------------------------------------------------------+
c
c  LOCAL SYMBOLIC DIMENSIONS (have big effect on storage requirements):
c
c       MXCLY  = Max no. of computational layers
c       MXULV  = Max no. of output levels
c       MXCMU  = Max no. of computation polar angles
c       MXUMU  = Max no. of output polar angles
c       MXPHI  = Max no. of output azimuthal angles
c       MXSQT  = Max no. of square roots of integers (for LEPOLY)
c +-------------------------------------------------------------------+

c     .. Parameters ..


      INTEGER   MXCLY, MXULV, MXCMU, MXUMU, MXPHI, MI, MI9M2, NNLYRI,
     &          MXSQT
      parameter (nstrms=40) ! pjr
      parameter (mxly=50)   ! pjr
      PARAMETER ( MXCLY=mxly, MXULV = mxly+1,
     &            MXCMU = nstrms, MXUMU = nstrms, MXPHI = nstrms,
     &            MI = MXCMU / 2, MI9M2 = 9*MI - 2,
     &            NNLYRI = MXCMU*MXCLY, MXSQT = 1000 )
c     ..
c     .. Scalar Arguments ..

      CHARACTER HEADER*127
      LOGICAL   DELTAM, LAMBER, ONLYFL, PLANK, USRANG, USRTAU
      INTEGER   IBCND, MAXCLY, MAXCMU, MAXPHI, MAXULV, MAXUMU, NLYR,
     &          NPHI, NSTR, NTAU, NUMU
      REAL      ACCUR, ALBEDO, BTEMP, FBEAM, FISOT, PHI0, TEMIS, TTEMP,
     &          UMU0, WVNMHI, WVNMLO
c     ..
c     .. Array Arguments ..

      LOGICAL   PRNT( 7 )
      REAL      ALBMED( MAXUMU ), DFDT( MAXULV ), DTAUC( MAXCLY ),
     &          FLUP( MAXULV ), HL( 0:MAXCMU ), PHI( MAXPHI ),
     &          PMOM( 0:MAXCMU, MAXCLY ), RFLDIR( MAXULV ),
     &          RFLDN( MAXULV ), SSALB( MAXCLY ), TEMPER( 0:MAXCLY ),
     &          TRNMED( MAXUMU ), U0U( MAXUMU, MAXULV ), UAVG( MAXULV ),
     &          UMU( MAXUMU ), UTAU( MAXULV ),
     &          UU( MAXUMU, MAXULV, MAXPHI )
c     ..
c     .. Local Scalars ..

      LOGICAL   COMPAR, LYRCUT, PASS1
      INTEGER   IQ, IU, J, KCONV, L, LC, LEV, LU, MAZIM, NAZ, NCOL,
     &          NCOS, NCUT, NN, NS
      REAL      ANGCOS, AZERR, AZTERM, BPLANK, COSPHI, DELM0, DITHER,
     &          DUM, PI, RPD, SGN, TPLANK
c     ..
c     .. Local Arrays ..

      INTEGER   IPVT( NNLYRI ), LAYRU( MXULV )

      REAL      AMB( MI, MI ), APB( MI, MI ), ARRAY( MXCMU, MXCMU ),
     &          B( NNLYRI ), BDR( MI, 0:MI ), BEM( MI ),
     &          CBAND( MI9M2, NNLYRI ), CC( MXCMU, MXCMU ),
     &          CMU( MXCMU ), CWT( MXCMU ), DTAUCP( MXCLY ),
     &          EMU( MXUMU ), EVAL( MI ), EVECC( MXCMU, MXCMU ),
     &          EXPBEA( 0:MXCLY ), FLDIR( MXULV ), FLDN( MXULV ),
     &          FLYR( MXCLY ), GC( MXCMU, MXCMU, MXCLY ),
     &          GL( 0:MXCMU, MXCLY ), GU( MXUMU, MXCMU, MXCLY ),
     &          HLPR( 0:MXCMU ), KK( MXCMU, MXCLY ), LL( MXCMU, MXCLY ),
     &          OPRIM( MXCLY ), PHIRAD( MXPHI ), PKAG( 0:MXCLY ),
     &          PSI( MXCMU ), RMU( MXUMU, 0:MI ), TAUC( 0:MXCLY ),
     &          TAUCPR( 0:MXCLY ), U0C( MXCMU, MXULV ), UTAUPR( MXULV ),
     &          UUM( MXUMU, MXULV ), WK( MXCMU ), XR0( MXCLY ),
     &          XR1( MXCLY ), YLM0( 0:MXCMU ), YLMC( 0:MXCMU, MXCMU ),
     &          YLMU( 0:MXCMU, MXUMU ), Z( NNLYRI ), Z0( MXCMU ),
     &          Z0U( MXUMU, MXCLY ), Z1( MXCMU ), Z1U( MXUMU, MXCLY ),
     &          ZBEAM( MXUMU, MXCLY ), ZJ( MXCMU ),
     &          ZPLK0( MXCMU, MXCLY ), ZPLK1( MXCMU, MXCLY ),
     &          ZZ( MXCMU, MXCLY ), SQT( MXSQT )

      DOUBLE PRECISION AAD( MI, MI ), EVALD( MI ), EVECCD( MI, MI ),
     &                 WKD( MXCMU )
c     ..
c     .. External Functions ..

      REAL      PLKAVG, R1MACH, RATIO
      EXTERNAL  PLKAVG, R1MACH, RATIO
c     ..
c     .. External Subroutines ..

      EXTERNAL  ALBTRN, CHEKIN, CMPINT, FLUXES, LEPOLY, PRAVIN, PRTINP,
     &          PRTINT, SETDIS, SETMTX, SLFTST, SOLEIG, SOLVE0, SURFAC,
     &          TERPEV, TERPSO, UPBEAM, UPISOT, USRINT, ZEROAL, ZEROIT
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ABS, ASIN, COS, LEN, MAX
c     ..
      SAVE      PASS1, PI, DITHER, RPD, SQT
      DATA      PASS1 / .TRUE. /


      IF( PASS1 ) THEN

         PI     = 2.*ASIN( 1.0 )
         DITHER = 10.*R1MACH( 4 )

c                            ** Must dither more on Cray (14-digit prec)

         IF( DITHER.LT.1.E-10 ) DITHER = 10.*DITHER

         RPD  = PI / 180.0

         DO 5 NS = 1, MXSQT
            SQT( NS ) = SQRT( FLOAT( NS ) )
    5    CONTINUE
c                            ** Set input values for self-test.
c                            ** Be sure SLFTST sets all print flags off.
         COMPAR = .FALSE.

c         CALL SLFTST( ACCUR, ALBEDO, BTEMP, DELTAM, DTAUC( 1 ), FBEAM,
c     &                FISOT, IBCND, LAMBER, NLYR, PLANK, NPHI, NUMU,
c     &                NSTR, NTAU, ONLYFL, PHI( 1 ), PHI0, PMOM( 0,1 ),
c     &                PRNT, SSALB( 1 ), TEMIS, TEMPER( 0 ), TTEMP,
c     &                UMU( 1 ), USRANG, USRTAU, UTAU( 1 ), UMU0, WVNMHI,
c     &                WVNMLO, COMPAR, DUM, DUM, DUM, DUM )
      END IF


   10 CONTINUE

c      IF( .NOT.PASS1 .AND. LEN(HEADER).NE.0 ) 
c     &    WRITE( *,'(//,1X,100(''*''),/,A,/,1X,100(''*''))' )
c     &    ' DISORT: '//HEADER

c                                  ** Calculate cumulative optical depth
c                                  ** and dither single-scatter albedo
c                                  ** to improve numerical behavior of
c                                  ** eigenvalue/vector computation
      CALL ZEROIT( TAUC, MXCLY + 1 )

      DO 20 LC = 1, NLYR

         IF( SSALB( LC ).EQ.1.0 ) SSALB( LC ) = 1.0 - DITHER
         TAUC( LC ) = TAUC( LC - 1 ) + DTAUC( LC )

   20 CONTINUE
c                                ** Check input dimensions and variables

      CALL CHEKIN( NLYR, DTAUC, SSALB, PMOM, TEMPER, WVNMLO, WVNMHI,
     &             USRTAU, NTAU, UTAU, NSTR, USRANG, NUMU, UMU, NPHI,
     &             PHI, IBCND, FBEAM, UMU0, PHI0, FISOT, LAMBER, ALBEDO,
     &             HL, BTEMP, TTEMP, TEMIS, PLANK, ONLYFL, ACCUR, TAUC,
     &             MAXCLY, MAXULV, MAXUMU, MAXCMU, MAXPHI, MXCLY, MXULV,
     &             MXUMU, MXCMU, MXPHI, MXSQT )

c                                 ** Zero internal and output arrays

      CALL  ZEROAL( MXCLY, EXPBEA(1), FLYR, OPRIM, TAUCPR(1), XR0, XR1,
     &              MXCMU, CMU, CWT, PSI, WK, Z0, Z1, ZJ,
     &              MXCMU+1, HLPR, YLM0,
     &              MXCMU**2, ARRAY, CC, EVECC,
     &              (MXCMU+1)*MXCLY, GL,
     &              (MXCMU+1)*MXCMU, YLMC,
     &              (MXCMU+1)*MXUMU, YLMU,
     &              MXCMU*MXCLY, KK, LL, ZZ, ZPLK0, ZPLK1,
     &              MXCMU**2*MXCLY, GC,
     &              MXULV, LAYRU, UTAUPR,
     &              MXUMU*MXCMU*MXCLY, GU,
     &              MXUMU*MXCLY, Z0U, Z1U, ZBEAM,
     &              MI, EVAL,
     &              MI**2, AMB, APB,
     &              NNLYRI, IPVT, Z,
     &              MAXULV, RFLDIR, RFLDN, FLUP, UAVG, DFDT,
     &              MAXUMU, ALBMED, TRNMED,
     &              MAXUMU*MAXULV, U0U,
     &              MAXUMU*MAXULV*MAXPHI, UU )

c                                 ** Perform various setup operations

      CALL SETDIS( CMU, CWT, DELTAM, DTAUC, DTAUCP, EXPBEA, FBEAM, FLYR,
     &             GL, HL, HLPR, IBCND, LAMBER, LAYRU, LYRCUT, MAXUMU,
     &             MAXCMU, MXCMU, NCUT, NLYR, NTAU, NN, NSTR, PLANK,
     &             NUMU, ONLYFL, OPRIM, PMOM, SSALB, TAUC, TAUCPR, UTAU,
     &             UTAUPR, UMU, UMU0, USRTAU, USRANG )

c                                 ** Print input information
      IF ( PRNT(1) )
     &     CALL PRTINP( NLYR, DTAUC, DTAUCP, SSALB, PMOM, TEMPER,
     &                  WVNMLO, WVNMHI, NTAU, UTAU, NSTR, NUMU, UMU,
     &                  NPHI, PHI, IBCND, FBEAM, UMU0, PHI0, FISOT,
     &                  LAMBER, ALBEDO, HL, BTEMP, TTEMP, TEMIS,
     &                  DELTAM, PLANK, ONLYFL, ACCUR, FLYR, LYRCUT,
     &                  OPRIM, TAUC, TAUCPR, MAXCMU, PRNT(7) )

c                              ** Handle special case for getting albedo
c                              ** and transmissivity of medium for many
c                              ** beam angles at once
      IF( IBCND.EQ.1 ) THEN

         CALL ALBTRN( ALBEDO, AMB, APB, ARRAY, B, BDR, CBAND, CC, 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 )
         RETURN

      END IF
c                                   ** Calculate Planck functions
      IF( .NOT.PLANK ) THEN

         BPLANK = 0.0
         TPLANK = 0.0
         CALL ZEROIT( PKAG, MXCLY + 1 )

      ELSE

         TPLANK = TEMIS*PLKAVG( WVNMLO, WVNMHI, TTEMP )
         BPLANK =       PLKAVG( WVNMLO, WVNMHI, BTEMP )

         DO 30 LEV = 0, NLYR
            PKAG( LEV ) = PLKAVG( WVNMLO, WVNMHI, TEMPER( LEV ) )
   30    CONTINUE

      END IF


c ========  BEGIN LOOP TO SUM AZIMUTHAL COMPONENTS OF INTENSITY  =======
c           (EQ STWJ 5)

      KCONV  = 0
      NAZ  = NSTR - 1
c                                    ** Azimuth-independent case

      IF( FBEAM.EQ.0.0 .OR. ( 1.- UMU0 ).LT.1.E-5 .OR. ONLYFL .OR.
     &      ( NUMU.EQ.1 .AND. ( 1.- UMU(1) ).LT.1.E-5 ) )
     &   NAZ = 0


      DO 160 MAZIM = 0, NAZ

         IF( MAZIM.EQ.0 ) DELM0  = 1.0
         IF( MAZIM.GT.0 ) DELM0  = 0.0

c                             ** Get normalized associated Legendre
c                             ** polynomials for
c                             ** (a) incident beam angle cosine
c                             ** (b) computational and user polar angle
c                             **     cosines
         IF( FBEAM.GT.0.0 ) THEN

            NCOS   = 1
            ANGCOS = -UMU0

            CALL LEPOLY( NCOS, MAZIM, MXCMU, NSTR-1, ANGCOS, SQT, YLM0 )

         END IF


         IF( .NOT.ONLYFL .AND. USRANG )
     &       CALL LEPOLY( NUMU, MAZIM, MXCMU, NSTR-1, UMU, SQT, YLMU )

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

c                       ** Get normalized associated Legendre polys.
c                       ** with negative arguments from those with
c                       ** positive arguments; Dave/Armstrong Eq. (15)
         SGN  = - 1.0

         DO 50 L = MAZIM, NSTR - 1

            SGN  = - SGN

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

   50    CONTINUE
c                                 ** Specify users bottom reflectivity
c                                 ** and emissivity properties
      IF ( .NOT.LYRCUT )
     &   CALL  SURFAC( ALBEDO, DELM0, FBEAM, HLPR, LAMBER,
     &                 MI, MAZIM, MXCMU, MXUMU, NN, NUMU, NSTR, ONLYFL,
     &                 UMU, USRANG, YLM0, YLMC, YLMU, BDR, EMU, BEM,
     &                 RMU, SQT )


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

         DO 60 LC = 1, NCUT

c                        ** Solve eigenfunction problem in Eq. STWJ(8B);
c                        ** return eigenvalues and eigenvectors

            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                                  ** Calculate particular solutions of
c                                  ** Eq.SS(18) for incident beam source
         IF ( FBEAM.GT.0.0 )
     &        CALL  UPBEAM( ARRAY, CC, CMU, DELM0, FBEAM, GL(0,LC),
     &                      IPVT, MAZIM, MXCMU, NN, NSTR, PI, UMU0, WK,
     &                      YLM0, YLMC, ZJ, ZZ(1,LC) )

c                              ** Calculate particular solutions of
c                              ** Eq. SS(15) for thermal emission source

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

               XR1( LC ) = 0.0

               IF( DTAUCP( LC ).GT.0.0 ) XR1( LC ) =
     &             ( PKAG( LC ) - PKAG( LC-1 ) ) / DTAUCP( LC )

               XR0( LC ) = PKAG( LC-1 ) - XR1( LC )*TAUCPR( LC-1 )

               CALL UPISOT( ARRAY, CC, CMU, IPVT, MXCMU, NN, NSTR,
     &                      OPRIM( LC ), WK, XR0( LC ), XR1( LC ), Z0,
     &                      Z1, ZPLK0( 1,LC ), ZPLK1( 1,LC ) )
            END IF


            IF( .NOT.ONLYFL .AND. USRANG ) THEN

c                                            ** Interpolate eigenvectors
c                                            ** to user angles

               CALL TERPEV( CWT, EVECC, GL( 0,LC ), GU( 1,1,LC ), MAZIM,
     &                      MXCMU, MXUMU, NN, NSTR, NUMU, WK, YLMC,
     &                      YLMU )
c                                            ** Interpolate source terms
c                                            ** to user angles

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

   60    CONTINUE

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


c                      ** Set coefficient matrix of equations combining
c                      ** boundary and layer interface conditions

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

c                      ** Solve for constants of integration in homo-
c                      ** geneous solution (general boundary conditions)

         CALL SOLVE0( B, BDR, BEM, BPLANK, CBAND, CMU, CWT, EXPBEA,
     &                FBEAM, FISOT, IPVT, LAMBER, LL, LYRCUT, MAZIM, MI,
     &                MI9M2, MXCMU, NCOL, NCUT, NN, NSTR, NNLYRI, PI,
     &                TPLANK, TAUCPR, UMU0, Z, ZZ, ZPLK0, ZPLK1 )

c                                  ** Compute upward and downward fluxes

      IF ( MAZIM.EQ.0 )
     &     CALL FLUXES( CMU, CWT, FBEAM, GC, KK, LAYRU, LL, LYRCUT,
     &                  MAXULV, MXCMU, MXULV, NCUT, NN, NSTR, NTAU,
     &                  PI, PRNT, SSALB, TAUCPR, UMU0, UTAU, UTAUPR,
     &                  XR0, XR1, ZZ, ZPLK0, ZPLK1, DFDT, FLUP,
     &                  FLDN, FLDIR, RFLDIR, RFLDN, UAVG, U0C )

         IF( ONLYFL ) THEN

            IF( MAXUMU.GE.NSTR ) THEN
c                                     ** Save azimuthal-avg intensities
c                                     ** at quadrature angles
               DO 80 LU = 1, NTAU

                  DO 70 IQ = 1, NSTR
                     U0U( IQ, LU ) = U0C( IQ, LU )
   70             CONTINUE

   80          CONTINUE

            END IF

            GO TO  170

         END IF


         CALL ZEROIT( UUM, MXUMU*MXULV )

         IF( USRANG ) THEN
c                                     ** Compute azimuthal intensity
c                                     ** components at user angles

            CALL USRINT( BPLANK, CMU, CWT, DELM0, DTAUCP, EMU, EXPBEA,
     &                   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 )

         ELSE
c                                     ** Compute azimuthal intensity
c                                     ** components at quadrature angles

            CALL CMPINT( FBEAM, GC, KK, LAYRU, LL, LYRCUT, MAZIM, MXCMU,
     &                   MXULV, MXUMU, NCUT, NN, NSTR, PLANK, NTAU,
     &                   TAUCPR, UMU0, UTAUPR, ZZ, ZPLK0, ZPLK1, UUM )
         END IF


         IF( MAZIM.EQ.0 ) THEN
c                               ** Save azimuthally averaged intensities

            DO 110 LU = 1, NTAU

               DO 100 IU = 1, NUMU
                  U0U( IU, LU ) = UUM( IU, LU )

                  DO 90 J = 1, NPHI
                     UU( IU, LU, J ) = UUM( IU, LU )
   90             CONTINUE

  100          CONTINUE

  110       CONTINUE
c                              ** Print azimuthally averaged intensities
c                              ** at user angles

            IF( PRNT( 4 ) ) CALL PRAVIN( UMU, NUMU, MAXUMU, UTAU, NTAU,
     &                                   U0U )
            IF( NAZ.GT.0 ) THEN

               CALL ZEROIT( PHIRAD, MXPHI )
               DO 120 J = 1, NPHI
                  PHIRAD( J ) = RPD*( PHI( J ) - PHI0 )
  120          CONTINUE

            END IF


         ELSE
c                                ** Increment intensity by current
c                                ** azimuthal component (Fourier
c                                ** cosine series);  Eq SD(2)
            AZERR  = 0.0

            DO 150 J = 1, NPHI

               COSPHI = COS( MAZIM*PHIRAD( J ) )

               DO 140 LU = 1, NTAU

                  DO 130 IU = 1, NUMU
                     AZTERM = UUM( IU, LU )*COSPHI
                     UU( IU, LU, J ) = UU( IU, LU, J ) + AZTERM
                     AZERR = MAX( AZERR,
     &                       RATIO( ABS(AZTERM), ABS(UU(IU,LU,J)) ) )
  130             CONTINUE

  140          CONTINUE

  150       CONTINUE

            IF( AZERR.LE.ACCUR ) KCONV  = KCONV + 1

            IF( KCONV.GE.2 ) GO TO  170

         END IF

  160 CONTINUE

c ===================  END LOOP ON AZIMUTHAL COMPONENTS  ===============


c                                          ** Print intensities
  170 CONTINUE
      IF( PRNT( 5 ) .AND. .NOT.ONLYFL ) CALL PRTINT( UU, UTAU, NTAU,
     &    UMU, NUMU, PHI, NPHI, MAXULV, MAXUMU )


      IF( PASS1 ) THEN
c                                    ** Compare test case results with
c                                    ** correct answers and abort if bad
         COMPAR = .TRUE.

c         CALL SLFTST( ACCUR, ALBEDO, BTEMP, DELTAM, DTAUC( 1 ), FBEAM,
c     &                FISOT, IBCND, LAMBER, NLYR, PLANK, NPHI, NUMU,
c     &                NSTR, NTAU, ONLYFL, PHI( 1 ), PHI0, PMOM( 0,1 ),
c     &                PRNT, SSALB( 1 ), TEMIS, TEMPER( 0 ), TTEMP,
c     &                UMU( 1 ), USRANG, USRTAU, UTAU( 1 ), UMU0, WVNMHI,
c     &                WVNMLO, COMPAR, FLUP( 1 ), RFLDIR( 1 ),
c     &                RFLDN( 1 ), UU( 1,1,1 ) )

         PASS1  = .FALSE.
         GO TO  10

      END IF

      RETURN
      END