SUBROUTINE FLUXES( CMU, CWT, FBEAM, GC, KK, LAYRU, LL, LYRCUT, 1,3
     &                   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 )

c       Calculates the radiative fluxes, mean intensity, and flux
c       derivative with respect to optical depth from the m=0 intensity
c       components (the azimuthally-averaged intensity)
c
c
c    I N P U T     V A R I A B L E S:
c
c       CMU      :  Abscissae for Gauss quadrature over angle cosine
c
c       CWT      :  Weights for Gauss quadrature over angle cosine
c
c       GC       :  Eigenvectors at polar quadrature angles, SC(1)
c
c       KK       :  Eigenvalues of coeff. matrix in Eq. SS(7)
c
c       LAYRU    :  Layer number of user level UTAU
c
c       LL       :  Constants of integration in Eq. SC(1), obtained
c                   by solving scaled version of Eq. SC(5);
c                   exponential term of Eq. SC(12) not included
c
c       LYRCUT   :  Logical flag for truncation of comput. layer
c
c       NN       :  Order of double-Gauss quadrature (NSTR/2)
c
c       NCUT     :  Number of computational layer where absorption
c                   optical depth exceeds ABSCUT
c
c       TAUCPR   :  Cumulative optical depth (delta-M-scaled)
c
c       UTAUPR   :  Optical depths of user output levels in delta-M
c                   coordinates;  equal to UTAU if no delta-M
c
c       XR0      :  Expansion of thermal source function in Eq. SS(14)
c
c       XR1      :  Expansion of thermal source function Eqs. SS(16)
c
c       ZZ       :  Beam source vectors in Eq. SS(19)
c
c       ZPLK0    :  Thermal source vectors Z0, by solving Eq. SS(16)
c
c       ZPLK1    :  Thermal source vectors Z1, by solving Eq. SS(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       U0C      :  Azimuthally averaged intensities
c                   ( at polar quadrature angles )
c
c       (RFLDIR, RFLDN, FLUP, DFDT, UAVG are DISORT output variables)
c
c
c    I N T E R N A L       V A R I A B L E S:
c
c       DIRINT   :  Direct intensity attenuated
c       FDNTOT   :  Total downward flux (direct + diffuse)
c       FLDIR    :  Direct-beam flux (delta-M scaled)
c       FLDN     :  Diffuse down-flux (delta-M scaled)
c       FNET     :  Net flux (total-down - diffuse-up)
c       FACT     :  EXP( - UTAUPR / UMU0 )
c       PLSORC   :  Planck source function (thermal)
c       ZINT     :  Intensity of m = 0 case, in Eq. SC(1)
c
c   Called by- DISORT
c   Calls- ZEROIT
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

      LOGICAL   LYRCUT
      INTEGER   MAXULV, MXCMU, MXULV, NCUT, NN, NSTR, NTAU
      REAL      FBEAM, PI, UMU0
c     ..
c     .. Array Arguments ..

      LOGICAL   PRNT( * )
      INTEGER   LAYRU( MXULV )
      REAL      CMU( MXCMU ), CWT( MXCMU ), DFDT( MAXULV ),
     &          FLDIR( MXULV ), FLDN( MXULV ), FLUP( MAXULV ),
     &          GC( MXCMU, MXCMU, * ), KK( MXCMU, * ), LL( MXCMU, * ),
     &          RFLDIR( MAXULV ), RFLDN( MAXULV ), SSALB( * ),
     &          TAUCPR( 0:* ), U0C( MXCMU, MXULV ), UAVG( MAXULV ),
     &          UTAU( MAXULV ), UTAUPR( MXULV ), XR0( * ), XR1( * ),
     &          ZPLK0( MXCMU, * ), ZPLK1( MXCMU, * ), ZZ( MXCMU, * )
c     ..
c     .. Local Scalars ..

      INTEGER   IQ, JQ, LU, LYU
      REAL      ANG1, ANG2, DIRINT, FACT, FDNTOT, FNET, PLSORC, ZINT
c     ..
c     .. External Subroutines ..

      EXTERNAL  ZEROIT
c     ..
c     .. Intrinsic Functions ..

      INTRINSIC ACOS, EXP
c     ..


      IF( PRNT(2) )  WRITE(*,'(//,21X,A,/,2A,/,2A,/)')
     & '<----------------------- FLUXES ----------------------->',
     & '   Optical  Compu    Downward    Downward    Downward     ',
     & ' Upward                    Mean      Planck   d(Net Flux)',
     & '     Depth  Layer      Direct     Diffuse       Total     ',
     & 'Diffuse         Net   Intensity      Source   / d(Op Dep)'

c                                        ** Zero DISORT output arrays
      CALL ZEROIT( U0C, MXULV*MXCMU )
      CALL ZEROIT( FLDIR, MXULV )
      CALL ZEROIT( FLDN, MXULV )

c                                        ** Loop over user levels
      DO 80 LU = 1, NTAU

         LYU  = LAYRU( LU )

         IF( LYRCUT .AND. LYU.GT.NCUT ) THEN
c                                                ** No radiation reaches
c                                                ** this level
            FDNTOT = 0.0
            FNET   = 0.0
            PLSORC = 0.0
            GO TO  70

         END IF


         IF( FBEAM.GT.0.0 ) THEN

            FACT         = EXP( -UTAUPR( LU ) / UMU0 )
            DIRINT       = FBEAM*FACT
            FLDIR( LU )  = UMU0*( FBEAM*FACT )
            RFLDIR( LU ) = UMU0*FBEAM * EXP( -UTAU( LU ) / UMU0 )

         ELSE

            DIRINT       = 0.0
            FLDIR( LU )  = 0.0
            RFLDIR( LU ) = 0.0

         END IF


         DO 30 IQ = 1, NN

            ZINT   = 0.0

            DO 10 JQ = 1, NN
               ZINT   = ZINT + GC( IQ, JQ, LYU )*LL( JQ, LYU )*
     &                  EXP( -KK( JQ,LYU )*( UTAUPR( LU ) -
     &                  TAUCPR( LYU ) ) )
   10       CONTINUE

            DO 20 JQ = NN + 1, NSTR
               ZINT   = ZINT + GC( IQ, JQ, LYU )*LL( JQ, LYU )*
     &                  EXP( -KK( JQ,LYU )*( UTAUPR( LU ) -
     &                  TAUCPR( LYU - 1 ) ) )
   20       CONTINUE

            U0C( IQ, LU ) = ZINT

            IF( FBEAM.GT.0.0 ) U0C( IQ, LU ) = ZINT + ZZ( IQ, LYU )*FACT

            U0C( IQ, LU ) = U0C( IQ, LU ) + ZPLK0( IQ, LYU ) +
     &                      ZPLK1( IQ, LYU )*UTAUPR( LU )
            UAVG( LU ) = UAVG( LU ) + CWT( NN + 1 - IQ )*U0C( IQ, LU )
            FLDN( LU ) = FLDN( LU ) + CWT( NN + 1 - IQ )*
     &                   CMU( NN + 1 - IQ )*U0C( IQ, LU )
   30    CONTINUE


         DO 60 IQ = NN + 1, NSTR

            ZINT   = 0.0

            DO 40 JQ = 1, NN
               ZINT   = ZINT + GC( IQ, JQ, LYU )*LL( JQ, LYU )*
     &                  EXP( -KK( JQ,LYU )*( UTAUPR( LU ) -
     &                  TAUCPR( LYU ) ) )
   40       CONTINUE

            DO 50 JQ = NN + 1, NSTR
               ZINT   = ZINT + GC( IQ, JQ, LYU )*LL( JQ, LYU )*
     &                  EXP( -KK( JQ,LYU )*( UTAUPR( LU ) -
     &                  TAUCPR( LYU - 1 ) ) )
   50       CONTINUE

            U0C( IQ, LU ) = ZINT

            IF( FBEAM.GT.0.0 ) U0C( IQ, LU ) = ZINT + ZZ( IQ, LYU )*FACT

            U0C( IQ, LU ) = U0C( IQ, LU ) + ZPLK0( IQ, LYU ) +
     &                      ZPLK1( IQ, LYU )*UTAUPR( LU )
            UAVG( LU ) = UAVG( LU ) + CWT( IQ - NN )*U0C( IQ, LU )
            FLUP( LU ) = FLUP( LU ) + CWT( IQ - NN )*CMU( IQ - NN )*
     &                   U0C( IQ, LU )
   60    CONTINUE


         FLUP( LU )  = 2.*PI*FLUP( LU )
         FLDN( LU )  = 2.*PI*FLDN( LU )
         FDNTOT      = FLDN( LU ) + FLDIR( LU )
         FNET        = FDNTOT - FLUP( LU )
         RFLDN( LU ) = FDNTOT - RFLDIR( LU )
         UAVG( LU )  = ( 2.*PI*UAVG( LU ) + DIRINT ) / ( 4.*PI )
         PLSORC      = XR0( LYU ) + XR1( LYU )*UTAUPR( LU )
         DFDT( LU )  = ( 1.- SSALB( LYU ) ) * 4.*PI *
     &                 ( UAVG( LU ) - PLSORC )

   70    CONTINUE
         IF( PRNT(2) ) WRITE(*,'(F10.4,I7,1P,7E12.3,E14.3)') UTAU( LU ),
     &       LYU, RFLDIR( LU ), RFLDN( LU ), FDNTOT, FLUP( LU ), FNET,
     &       UAVG( LU ), PLSORC, DFDT( LU )

   80 CONTINUE


      IF( PRNT(3) ) THEN

         WRITE(*,'(//,2A)') ' ******** AZIMUTHALLY AVERAGED ',
     &      'INTENSITIES ( at polar quadrature angles ) *******'

         DO 100 LU = 1, NTAU

            WRITE( *, '(/,A,F10.4,//,2A)' ) 
     &         ' Optical depth =', UTAU( LU ),
     &         '     Angle (deg)   cos(Angle)     Intensity',
     &         '     Angle (deg)   cos(Angle)     Intensity'

            DO 90 IQ = 1, NN
               ANG1 = (180./ PI) * ACOS( CMU( 2*NN - IQ + 1 ) )
               ANG2 = (180./ PI) * ACOS( CMU( IQ ) )
               WRITE(*,'(     &             ANG1, CMU(2*NN-IQ+1), U0C(IQ,LU),
     &             ANG2, CMU(IQ),        U0C(IQ+NN,LU)
   90       CONTINUE

  100    CONTINUE

      END IF

      RETURN
      END