      SUBROUTINE SWTRAN_MCICA (REFL, TRAN, ITILE,
     1                         UCUMDTR, TRAN0, TAUA, TAUR, TAUG,
     2                         TAUOMA, TAUOMGA, F1, F2, TAUC,
     3                         TAUOMC, TAUOMGC, CLD, RMU, C1, C2,
     4                         ALBSURT, CSALBT, NCT, CUT, LEV1,
     5                         IL1, IL2, ILG, LAY, LEV,
     6                         NTILE)
C
C     * Sept 3/2018 - J. COLE  Promoted more variables to real*8
C     *                        Hopefully this will make code more robust.
C
C     * Feb 10/2014 - J.LI.    Revised version for gcm18:
C     *                        - Compute the direct beam using the
C     *                          unscaled cloud optical properties.
C     * FEB 20/2007 - J. COLE: UPDATED CODE TO FOLLOW CHANGES MADE BETWEEN
C                              SWTRAN AND SWTRAN2
C
C     * APR 11/2005 - J. COLE: USED PETRI'S HARD WORK AND MODIFIED TO SUIT MY VERSION.
C                              NO LONGER HAS SUBCOLUMNS, WORKS ON ONLY ONE COLUMN PER
C                              GCM COLUMN.
C     * MAY 31/2004 - P.RAISANEN:
C
C     THIS CODE VERSION HAS BEEN UPDATED FOR McICA RADIATION CALCULATIONS
C
C       - Clear-sky reflectances and transmittances are calculated
C         using mean profiles for the GCM column, but cloudy-sky
C         reflectances and transmittances are calculated separately
C         for NSUB subcolumns and then averaged. It is assumed that
C         there are equally many subcolumns for each GCM column I1...I2.
C
C       - Treatment of cloud overlap and horizontal variability eliminated.
C         Cloud fraction for subcolumns (CLDSUB) is still included as
C         input, but the code assumes it is either 0 or 1.
C
C       - Security checks added to avoid values of X2 and X4 equal to 1.
C
C     * APR 25,2003 - J.LI.
C----------------------------------------------------------------------C
C     DELTA-EDDINGTON APPROXIMATION AND ADDING PROCESS FOR CLEAR AND   C
C     ALL SKY, THE ADDING METHOD BY COAKLEY ET AL (1983). THIS CODE    C
C     CAN DEAL WITH SOLAR RADIATIVE TRANSFER THROUGH ATMOSPHERE WITH   C
C     PROPER TREATMENT OF CLOUD OVERLAP (RANDOM + MAXIMUM OR RANDOM +  C
C     SLANTWISE) AND CLOUD SUB-GRID VARIABILITY. THE THEORY FOR ADDING,C
C     CLOUD OVERLAP LI AND DOBBIE (2003). CLOUD SUB-GRID VARIABILITY   C
C     SIMILAR TO WITH ADJUSTMENT OF CLOUD OPTICAL DEPTH                C
C                                                                      C
C     REFLT:     TILED REFLECTIVITY (1) CLEAR SKY; (2) ALL SKY         C
C     TRANT:     TILED TRANSMITIVITY                                   C
C     CUMDTR:    DIRECT TRANSMISSION FOR MULT-LAYERS                   C
C    UCUMDTR:    DIRECT TRANSMISSION FOR MULT-LAYERS USING UNSCALE     C
C                CLOUD PROPERTIES                                      C
C     TAUA:      AEROSOL OPTICAL DEPTH                                 C
C     TAUR:      RAYLEIGH OPTICAL DEPTH                                C
C     TAUG:      GASEOUS OPTICAL DEPTH                                 C
C     TAUOMA:    AEROSOL OPTICAL DEPTH TIMES AEROSOL SINGLE SCATTERING C
C                ALBEDO                                                C
C     TAUOMGA:   TAUOMA TIMES AEROSOL ASYMMETRY FACTOR                 C
C     F1:        SQUARE OF AEROSOL ASYMMETRY FACTOR                    C
C     F2:        SQUARE OF CLOUD ASYMMETRY FACTOR                      C
C     TAUC:      CLOUD OPTICAL DEPTH                                   C
C     TAUOMC:    CLOUD OPTICAL DEPTH TIMES CLOUD SINGLE                C
C                SCATTERING ALBEDO                                     C
C     TAUOMG:    TAUOMC TIMES CLOUD ASYMMETRY FACTOR                   C
C     CLD:       CLOUD FRACTION (ASSUMED TO BE 0 OR 1)                 C
C     RMU:       COS OF SOLAR ZENITH ANGLE                             C
C     C1 @ C2:   TWO FACTORS NOT DEPENDENT ON IB AND IG CALCULATED     C
C                OUTSIDE FOR EFFICIENCY                                C
C     ALBSURT:   TILED ALL-SKY SURFACE ALBEDO                          C
C     CSALBT:    TILED CLEAR-SKY SURFACE ALBEDO                        C
C     ITILE:     INDICATOR IF CALCULATIONS IS TO BE PERFORMED ON TILE  C
C     NCT:       THE HIGHEST CLOUD TOP LEVEL                           C
C     RDF:       LAYER DIFFUSE REFLECTION                              C
C     TDF:       LAYER DIFFUSE TRANSMISSION                            C
C     RDR:       LAYER DIRECT REFLECTION                               C
C     TDR:       LAYER DIRECT TRANSMISSION                             C
C     DTR:       DIRECT TRANSMISSION                                   C
C    UDTR:       DIRECT TRANSMISSION USING UNSCALE CLOUD PROPERTIES    C
C     RMDF:      DIFFUSE REFLECTION FROM MODEL TOP LEVEL               C
C     TMDR:      DIRECT TRANSMISSION FROM MODEL TOP LEVEL              C
C     RMUR:      DIRECT REFLECTION FROM MODEL BOTTOM LEVEL             C
C     RMUF:      DIFFUSE REFLECTION FROM MODEL BOTTOM LEVEL            C
C----------------------------------------------------------------------C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

C     
      REAL CUMDTR(ILG,2,LEV), UCUMDTR(ILG,2,LEV), TRAN0(ILG)
      REAL REFL(ILG,NTILE,2,LEV), TRAN(ILG,NTILE,2,LEV),
     1     ALBSURT(ILG,NTILE), CSALBT(ILG,NTILE)
      REAL TAUA(ILG,LAY), TAUR(ILG,LAY), TAUG(ILG,LAY), TAUOMA(ILG,LAY),
     1     TAUOMGA(ILG,LAY), F1(ILG,LAY), F2(ILG,LAY),
     2     TAUC(ILG,LAY), TAUOMC(ILG,LAY),
     3     TAUOMGC(ILG,LAY), CLD(ILG,LAY),
     3     RMU(ILG), C1(ILG), C2(ILG)
      REAL*8 RDF(ILG,2,LAY), TDF(ILG,2,LAY), RDR(ILG,2,LAY),
     1     TDR(ILG,2,LAY), DTR(ILG,2,LAY), UDTR(ILG,2,LAY)

      REAL*8 RMDF(ILG,NTILE,2,LEV), TMDR(ILG,NTILE,2,LEV),
     1       RMUR(ILG,NTILE,2,LEV), RMUF(ILG,NTILE,2,LEV)
      INTEGER NCT(ILG)
      INTEGER ITILE(ILG,NTILE)

c
c     * scalars promoted to 64-bit in loop 200 to avoid pole singularity
c     * from two-stream calculation in clear-sky.
c
      real*8 extopt,omars,ssalb,sf1,sf,tau1,om1,cow
      real*8 ssgas1,cowg,alamd,u,u2,uu,efun,efun2,x1,rn,x2
      real*8 y,yx,dm,gscw,appgm,apmgm,omarcs,sf2
      real*8 tau2,om2,ssgas2,sdtr,x3,x4,srdf,stdf,srdr,stdr
      real*8 x5,x6,x7,x,lowbnd

      real*8 rdf8,tdf8,rdr8,tdr8,dtr8,udtr8
      real*8 arg_alamd, arg_dtr, arg_udtr, arg_efun
      real*8 rmu8(ilg),c18(ilg),c28(ilg)
      real*8 taua8(ilg,lay), taur8(ilg,lay), taug8(ilg,lay), 
     +       tauoma8(ilg,lay),tauomga8(ilg,lay), f18(ilg,lay), 
     +       f28(ilg,lay), tauc8(ilg,lay),tauomc8(ilg,lay),
     +       tauomgc8(ilg,lay)

      real*8, parameter :: r_0_5=0.5_8, 
     +                     r_one=1.0_8, 
     +                     r_1_5=1.5_8, 
     +                     r_three=3.0_8,
     +                     r_zero=0.0_8,
     +                     r_mone=-1.0_8, 
     +                     eps_10=1.0D-10,
     +                     eps_20=1.0D-20

      data lowbnd /1.e-15/

C-------------------
C INITIALIZATIONS
C-------------------

      DO 10 K = LEV1, LEV
        DO 20 I = IL1,IL2
          CUMDTR(I,1:2,K) = r_zero
         UCUMDTR(I,1:2,K) = r_zero
 20     CONTINUE
 10   CONTINUE

      DO K = LEV1, LEV
        DO M = 1, NTILE
          DO I = IL1,IL2
            REFL(I,M,1:2,K) = r_zero
            TRAN(I,M,1:2,K) = r_zero
          END DO ! I
        END DO ! M
      END DO ! K

      DO I = IL1, IL2
         rmu8(i) = REAL(RMU(I),8)
         c18(i)  = REAL(C1(I),8)
         c28(i)  = REAL(C2(I),8)
      END DO

      DO K = 1, LAY
         DO I = IL1, IL2
            taua8(i,k)    = REAL(TAUA(I,K),8)
            taur8(i,k)    = REAL(TAUR(I,K),8)
            taug8(i,k)    = REAL(TAUG(I,K),8)
            tauoma8(i,k)  = REAL(TAUOMA(I,K),8)
            tauomga8(i,k) = REAL(TAUOMGA(I,K),8)
            tauc8(i,k)    = REAL(TAUC(I,K),8)
            tauomc8(i,k)  = REAL(TAUOMC(I,K),8)
            tauomgc8(i,k) = REAL(TAUOMGC(I,K),8)
            f18(i,k)      = REAL(F1(I,K),8)
            f28(i,k)      = REAL(F2(I,K),8)
         END DO ! I
      END DO ! K
C
C----------------------------------------------------------------------C
C     COMBINE THE OPTICAL PROPERTIES FOR SOLAR,                        C
C     1, AEROSOL + RAYLEIGH + GAS; 2, CLOUD + AEROSOL + RAYLEIGH + GAS C
C     CALCULATE THE DIRECT AND DIFFUSE REFLECTION AND TRANSMISSION IN  C
C     THE SCATTERING LAYERS USING THE DELTA-EDDINGTON METHOD.          C
C----------------------------------------------------------------------C
C
C----------------------------------------------------------------------C
C FIRST, COMPUTATION OF CLEAR-SKY FLUXES                   1           C
C----------------------------------------------------------------------C

      DO 200 K = LEV1, LAY
        DO 200 I = IL1, IL2
          EXTOPT          =  TAUA8(I,K) +  TAUR8(I,K) +  TAUG8(I,K)
          OMARS           =  TAUOMA8(I,K) + TAUR8(I,K)
          SSALB           =  OMARS / (EXTOPT + EPS_20)
          SF1             =  F18(I,K) / OMARS
          SF              =  SSALB * SF1
          TAU1            =  EXTOPT * (r_one - SF)
          OM1             = (SSALB - SF) / (r_one - SF)
          COW             =  r_one - OM1 + eps_10
          SSGAS1          = (TAUOMGA8(I,K) / OMARS - SF1) /
     1                      (r_one - SF1)
          COWG            =  r_one - OM1 * SSGAS1
          arg_alamd       =  r_three * COW * COWG
          ALAMD           =  SQRT(arg_alamd)
C
          arg_dtr         =  r_mone * TAU1 / RMU8(I)
          arg_udtr        =  r_mone * EXTOPT / RMU8(I)
          DTR8            =  EXP(arg_dtr)
         UDTR8            =  EXP(arg_udtr)          
          U               =  r_1_5 * COWG / ALAMD
          U2              =  U + U
          UU              =  U * U
          arg_efun        =  r_mone * ALAMD * TAU1
          EFUN            =  EXP(arg_efun)
          EFUN2           =  EFUN * EFUN
          X1              = (UU - U2 + r_one) * EFUN2
          RN              =  r_one / (UU + U2 + r_one - X1)
          X2              =  ALAMD * RMU8(I)
          Y               =  r_one - X2 * X2
          YX              =  SIGN( MAX( ABS(Y), LOWBND), Y )
          DM              =  OM1 / YX
          GSCW            =  SSGAS1 * COW
          APPGM           = (C18(I) + r_0_5 +
     1                       GSCW * (C18(I) + C28(I)))*DM
          APMGM           = (C18(I) - r_0_5 +
     1                       GSCW * (C18(I) - C28(I)))*DM
          RDF8            = (UU - r_one) * (r_one  - EFUN2) * RN
          TDF8            = (U2 + U2) * EFUN * RN
          RDR8            =  APPGM * RDF8 
     1                    +  APMGM*(TDF8 * DTR8 - r_one)
          TDR8            =  APPGM * TDF8 
     1                    + (APMGM * RDF8 - APPGM + r_one) * DTR8

          RDF(I,1,K)      = REAL(rdf8,4)
          TDF(I,1,K)      = REAL(tdf8,4)
          RDR(I,1,K)      = REAL(rdr8,4)
          TDR(I,1,K)      = REAL(tdr8,4)
          DTR(I,1,K)      = REAL(dtr8,4)
          UDTR(I,1,K)     = REAL(udtr8,4)

 200  CONTINUE

C
      DO 300 I = IL1, IL2
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE FIRST LEVEL (LEV1).                       C
C----------------------------------------------------------------------C
C
        CUMDTR(I,1,LEV1)        =  TRAN0(I)
       UCUMDTR(I,1,LEV1)        =  TRAN0(I)
 300  CONTINUE

      DO M = 1, NTILE
        DO I = IL1, IL2
          IF (ITILE(I,M) .GT. 0) THEN
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE FIRST LEVEL (LEV1).                       C
C----------------------------------------------------------------------C
C
            TMDR(I,M,1,LEV1)            =  TRAN0(I)
            RMDF(I,M,1,LEV1)            =  1.0 - TRAN0(I)
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE GROUND LAYER.                             C
C----------------------------------------------------------------------C
C
            RMUR(I,M,1,LEV)             =  CSALBT(I,M)
            RMUF(I,M,1,LEV)             =  CSALBT(I,M)
          END IF
        END DO ! I
      END DO ! M
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS DOWNWARD FROM THE SECOND LAYER TO THE SURFACE.    C
C----------------------------------------------------------------------C
C
      DO K = LEV1 + 1, LEV
        KM1 = K - 1
        DO I = IL1, IL2
           CUMDTR(I,1,K) =  CUMDTR(I,1,KM1) *  DTR(I,1,KM1)
          UCUMDTR(I,1,K) = UCUMDTR(I,1,KM1) * UDTR(I,1,KM1)
        END DO ! I
      END DO ! K

      DO 451 K = LEV1 + 1, LEV
        KM1 = K - 1
        L = LEV - K + LEV1
        LP1 = L + 1
        DO 450 M = 1, NTILE
          DO 400 I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
              DMM           = TDF(I,1,KM1) /
     1                        (1.0 - RDF(I,1,KM1)*RMDF(I,M,1,KM1))
              FMM           = RMDF(I,M,1,KM1) * DMM
              TMDR(I,M,1,K) = CUMDTR(I,1,KM1) * (TDR(I,1,KM1) +
     1                        RDR(I,1,KM1) * FMM) +
     2                        (TMDR(I,M,1,KM1) - CUMDTR(I,1,KM1))* DMM
              RMDF(I,M,1,K) =  RDF(I,1,KM1) + TDF(I,1,KM1) * FMM

C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD FROM ONE LAYER ABOVE SURFACE TO THE LEV1.  C
C----------------------------------------------------------------------C

              UMM           = TDF(I,1,L) /
     1                        (1.0 - RDF(I,1,L) * RMUF(I,M,1,LP1))
              FMM           = RMUF(I,M,1,LP1) * UMM
              RMUR(I,M,1,L) = RDR(I,1,L) + DTR(I,1,L) *
     1                        RMUR(I,M,1,LP1) * UMM +
     2                        (TDR(I,1,L) - DTR(I,1,L)) * FMM
              RMUF(I,M,1,L) = RDF(I,1,L) + TDF(I,1,L) * FMM
            END IF
 400      CONTINUE
 450    CONTINUE
 451  CONTINUE
C
C----------------------------------------------------------------------C
C     ADD DOWNWARD TO CALCULATE THE RESULTANT REFLECTANCES AND         C
C     TRANSMITTANCE AT FLUX LEVELS.                                    C
C----------------------------------------------------------------------C
C
      DO 550 K = LEV1, LEV
        DO 551 M = 1, NTILE
          DO 500 I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
              DMM           = 1.0 /
     1                       (1.0 - RMUF(I,M,1,K)*RMDF(I,M,1,K))
              X             = CUMDTR(I,1,K) * RMUR(I,M,1,K)
              Y             = TMDR(I,M,1,K) - CUMDTR(I,1,K)
              TRAN(I,M,1,K) = CUMDTR(I,1,K) +
     1                        (X * RMDF(I,M,1,K) + Y) * DMM
              REFL(I,M,1,K) = (X + Y * RMUF(I,M,1,K)) * DMM
            END IF
 500      CONTINUE
 551    CONTINUE
 550  CONTINUE

C----------------------------------------------------------------------C
C SECOND, COMPUTATION OF ALL-SKY FLUXES                                C
C----------------------------------------------------------------------C

C----------------------------------------------------------------------C
C ALL-SKY FLUXES.                                                      C
C----------------------------------------------------------------------C

      DO 600 K = LEV1, LAY
        DO 600 I = IL1, IL2
        IF (CLD(I,K) .LT. CUT)                               THEN
          RDF(I,2,K)      =  RDF(I,1,K)
          TDF(I,2,K)      =  TDF(I,1,K)
          RDR(I,2,K)      =  RDR(I,1,K)
          TDR(I,2,K)      =  TDR(I,1,K)
          DTR(I,2,K)      =  DTR(I,1,K)
         UDTR(I,2,K)      = UDTR(I,1,K)
        ELSE
          EXTOPT          =  TAUC8(I,K) + TAUA8(I,K) +
     1                       TAUR8(I,K) + TAUG8(I,K)
          OMARCS          =  TAUOMC8(I,K) + TAUR8(I,K)
          SSALB           =  OMARCS / EXTOPT
          SF2             =  F28(I,K) / OMARCS
          SF              =  SSALB * SF2
          TAU2            =  EXTOPT * (r_one - SF)
          OM2             = (SSALB - SF) / (r_one - SF)
          COW             =  r_one - OM2
          SSGAS2          = (TAUOMGC8(I,K)/ OMARCS - SF2) /
     1                      (r_one - SF2)
          COWG            =  r_one - OM2 * SSGAS2
          arg_alamd       =  r_three * COW * COWG
          ALAMD           =  SQRT(arg_alamd)
C
          arg_dtr         =  r_mone * TAU2 / RMU8(I)
          arg_udtr        =  r_mone * EXTOPT / RMU8(I)
          DTR8            =  EXP(arg_dtr)
         UDTR8            =  EXP(arg_udtr)
          U               =  r_1_5 * COWG / ALAMD
          U2              =  U + U
          UU              =  U * U
          arg_efun        =  r_mone * ALAMD * TAU2
          EFUN            =  EXP(arg_efun)
          EFUN2           =  EFUN * EFUN
          X3              = (UU - U2 + r_one) * EFUN2
          RN              =  r_one / (UU + U2 + r_one - X3)
          X4              =  ALAMD * RMU8(I)
          Y               =  r_one - X4 * X4
          YX              =  SIGN( MAX( ABS(Y), LOWBND), Y )
          DM              =  OM2 / YX
          GSCW            =  SSGAS2 * COW
          APPGM           = (C18(I) + r_0_5 +
     1                       GSCW * (C18(I) + C28(I))) * DM
          APMGM           = (C1(I) - r_0_5 +
     1                       GSCW * (C18(I) - C28(I))) * DM

          RDF8            = (UU - r_one) * (r_one - EFUN2) * RN
          TDF8            = (U2 + U2) * EFUN * RN
          RDR8            =  APPGM * RDF8 + APMGM *
     1                      (TDF8 * DTR8 - r_one)
          TDR8            =  APPGM * TDF8 +
     1                      (APMGM * RDF8 - APPGM + r_one) *
     2                       DTR8

          RDF(I,2,K)      = REAL(rdf8,4)
          TDF(I,2,K)      = REAL(tdf8,4)
          RDR(I,2,K)      = REAL(rdr8,4)
          TDR(I,2,K)      = REAL(tdr8,4)
          DTR(I,2,K)      = REAL(dtr8,4)
          UDTR(I,2,K)     = REAL(udtr8,4)

        ENDIF
 600  CONTINUE
C
      DO 700 I = IL1, IL2
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE FIRST LEVEL (LEV1).                       C
C----------------------------------------------------------------------C
C
        CUMDTR(I,2,LEV1)        =  TRAN0(I)
       UCUMDTR(I,2,LEV1)        =  TRAN0(I)

 700  CONTINUE

      DO M = 1, NTILE
        DO I = IL1, IL2
          IF (ITILE(I,M) .GT. 0) THEN
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE FIRST LEVEL (LEV1).                       C
C----------------------------------------------------------------------C
C

            TMDR(I,M,2,LEV1)     =  TRAN0(I)
            RMDF(I,M,2,LEV1)     =  1.0 - TRAN0(I)
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE GROUND LAYER.                             C
C----------------------------------------------------------------------C
C
            RMUR(I,M,2,LEV)      =  ALBSURT(I,M)
            RMUF(I,M,2,LEV)      =  ALBSURT(I,M)
          END IF
        END DO ! I
      END DO ! M
C----------------------------------------------------------------------C
C     ADD THE LAYERS DOWNWARD FROM THE SECOND LAYER TO THE SURFACE.    C
C----------------------------------------------------------------------C
C
      DO K = LEV1 + 1, LEV
        KM1 = K - 1
        DO I = IL1, IL2
          IF (NCT(I) .LE. LAY) THEN
            IF (K .LE. NCT(I)) THEN
              CUMDTR(I,2,K) =  CUMDTR(I,1,K)
             UCUMDTR(I,2,K) = UCUMDTR(I,1,K)
            ELSE
              CUMDTR(I,2,K) =  CUMDTR(I,2,KM1)* DTR(I,2,KM1)
             UCUMDTR(I,2,K) = UCUMDTR(I,2,KM1)*UDTR(I,2,KM1)
            END IF
          ELSE
            CUMDTR(I,2,K) = 0.0
           UCUMDTR(I,2,K) = 0.0
          END IF
        END DO ! I
      END DO ! K

      DO 851 K = LEV1 + 1, LEV
        KM1 = K - 1
        L = LEV - K + LEV1
        LP1 = L + 1
        DO 850 M = 1, NTILE
          DO 800 I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
              IF (NCT(I) .LE. LAY) THEN
                IF (K .LE. NCT(I)) THEN
                  TMDR(I,M,2,K) =  TMDR(I,M,1,K)
                  RMDF(I,M,2,K) =  RMDF(I,M,1,K)
                ELSE
                  DPP           = TDF(I,2,KM1) /
     1                            (1.0 - RMDF(I,M,2,KM1)*
     2                            RDF(I,2,KM1))
                  FPP           = RMDF(I,M,2,KM1) * DPP
                  TMDR(I,M,2,K) = CUMDTR(I,2,KM1) *
     1                            (TDR(I,2,KM1) +
     1                            RDR(I,2,KM1) * FPP) +
     2                            (TMDR(I,M,2,KM1) -
     3                            CUMDTR(I,2,KM1))*DPP
                  RMDF(I,M,2,K) = RDF(I,2,KM1) +
     1                            TDF(I,2,KM1) * FPP
                ENDIF
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD FROM ONE LAYER ABOVE SURFACE TO THE LEV1.  C
C----------------------------------------------------------------------C
C
                UPP           = TDF(I,2,L) /
     1                          (1.0 - RMUF(I,M,2,LP1) *
     2                          RDF(I,2,L))
                FPP           = RMUF(I,M,2,LP1) * UPP
                RMUR(I,M,2,L) = RDR(I,2,L) + DTR(I,2,L) *
     1                          RMUR(I,M,2,LP1)*UPP +
     2                          (TDR(I,2,L) - DTR(I,2,L))*FPP
                RMUF(I,M,2,L) = RDF(I,2,L) + TDF(I,2,L)*FPP
              ELSE
                TMDR(I,M,2,K) = 1.0
                RMDF(I,M,2,K) = 0.0
                RMUR(I,M,2,L) = 0.0
                RMUF(I,M,2,L) = 0.0
              ENDIF
            END IF
 800      CONTINUE
 850    CONTINUE
 851  CONTINUE

C
C----------------------------------------------------------------------C
C     ADD DOWNWARD TO CALCULATE THE RESULTANT REFLECTANCES AND         C
C     TRANSMITTANCE AT FLUX LEVELS.                                    C
C----------------------------------------------------------------------C
C

      DO 950 K = LEV1, LEV
        DO 951 M = 1, NTILE
          DO 900 I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
              IF (NCT(I) .LE. LAY)                          THEN
                DPP           = 1.0 /
     1                         (1.0-RMUF(I,M,2,K)*RMDF(I,M,2,K))
                X             = CUMDTR(I,2,K) * RMUR(I,M,2,K)
                Y             = TMDR(I,M,2,K) - CUMDTR(I,2,K)
                TRAN(I,M,2,K) = CUMDTR(I,2,K)
     1                         + (X * RMDF(I,M,2,K) + Y) * DPP
                REFL(I,M,2,K) = (X + Y * RMUF(I,M,2,K)) * DPP
              ELSE
                TRAN(I,M,2,K) = TRAN(I,M,1,K)
                REFL(I,M,2,K) = REFL(I,M,1,K)
              ENDIF
            END IF
 900      CONTINUE
 951    CONTINUE
 950  CONTINUE

      RETURN
      END
