      SUBROUTINE SWTRAN2(REFL, TRAN, CUMDTR, TRAN0, TAUA, TAUR, TAUG,
     1                   TAUOMA, TAUOMGA, F1, F2, TAUCS, TAUOMC, 
     2                   TAUOMGC, CLD, CLDM, A1, RMU, C1, C2, ALBSUR,
     3                   CSALB, NBLK, NCT, 
     4                   CUT, LEV1, IL1, IL2, ILG, LAY, LEV)
C
C     * NOV 22/2006 - M.LAZARE. NEW VERSION BASED ON SWTRAN:
C     *                         - PROMOTION OF SELECTED ARRAYS/VARIABLES
C     *                           TO REAL*8 FOR SUFFICIENT ACCURACY
C     *                           WHEN RUNNING IN 32-BIT.
C     *                         - BOUND OF LOWBND USED FOR DENOMINATOR
C     *                           FOR DM (USING "SIGN") IN 2 PLACES.
C     *                         - WORK ARAYS FOR THIS ROUTINE NOW
C     *                           ARE AUTOMATIC ARRAYS AND THEIR
C     *                           SPACE IS NOT PASSED IN.   
C     * MAY 14,2003 - J.LI.     PREVIOUS VERSION SWTRAN UP TO GCM15E.
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     REFL:    REFLECTIVITY (1) CLEAR SKY; (2) ALL SKY                 C
C     TRAN:    TRANSMITIVITY                                           C
C     CUMDTR:  DIRECT TRANSMISSION FOR MULT-LAYERS                     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     TAUCS:   CLOUD OPTICAL DEPTH                                     C
C     TAUOMC:  CLOUD OPTICAL DEPTH TIMES CLOUD SINGLE SCATTERING ALBEDOC
C     TAUOMGC: TAUOMC TIMES CLOUD ASYMMETRY FACTOR                     C
C     CLD:     CLOUD FRACTION                                          C
C     CLDM:    MAXIMUM PORTION IN EACH CLOUD BLOCK, IN WHICH THE EXACT C
C              SOLUTION FOR SUBGRID VARIABILITY IS APPLIED             C
C     A1:      VARIOUS RELATION FOR CLOUD OVERLAP                      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     ALBSUR:  SURFACE ALBEDO                                          C
C     CSALB:   CLEAR-SKY SURFACE ALBEDO                                C
C     NBLK:    NUMBER OF CLOUD BLOCKS ACCOUNTED FROM SURFACE           C
C     NCT:     THE HIGHEST CLOUD TOP LEVEL FOR THE LONGITUDE AND       C
C              LATITUDE LOOP (ILG)                                     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     RMDF:    BLOCK DIFFUSE REFLECTION FROM MODEL TOP LEVEL           C
C     TMDR:    BLOCK DIRECT TRANSMISSION FROM MODEL TOP LEVEL          C
C     RMUR:    BLOCK DIRECT REFLECTION FROM MODEL BOTTOM LEVEL         C
C     RMUF:    BLOCK DIFFUSE REFLECTION FROM MODEL BOTTOM LEVEL        C
C----------------------------------------------------------------------C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL   REFL(ILG,2,LEV), TRAN(ILG,2,LEV), CUMDTR(ILG,4,LEV), 
     1       TRAN0(ILG)
      REAL   TAUA(ILG,LAY),TAUR(ILG,LAY),TAUG(ILG,LAY),TAUOMA(ILG,LAY),
     1       TAUOMGA(ILG,LAY),F1(ILG,LAY),F2(ILG,LAY),TAUCS(ILG,LAY),
     2       TAUOMC(ILG,LAY),TAUOMGC(ILG,LAY),CLD(ILG,LAY), 
     3       CLDM(ILG,LAY),A1(ILG,11),RMU(ILG),C1(ILG),C2(ILG),
     4       ALBSUR(ILG),CSALB(ILG)
      REAL*8 RDF (ILG,4,LAY), TDF (ILG,4,LAY), RDR (ILG,4,LAY),
     1       TDR (ILG,4,LAY), DTR (ILG,4,LAY), RMDF(ILG,4,LEV),
     2       TMDR(ILG,4,LEV), RMUR(ILG,4,LEV), RMUF(ILG,4,LEV)
      INTEGER NBLK(ILG, LAY), NCT(ILG)
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
c
      data lowbnd /1.e-15/
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
      DO 200 K = LEV1, LAY
      DO 200 I = IL1, IL2
        EXTOPT                    =  TAUA(I,K) + TAUR(I,K) + TAUG(I,K)
        OMARS                     =  TAUOMA(I,K) + TAUR(I,K)
        SSALB                     =  OMARS / (EXTOPT + 1.E-20)
        SF1                       =  F1(I,K) / OMARS
        SF                        =  SSALB * SF1
        TAU1                      =  EXTOPT * (1.0 - SF)
        OM1                       = (SSALB - SF) / (1.0 - SF)
        COW                       =  1.0 - OM1 + 1.E-10
        SSGAS1                    = (TAUOMGA(I,K) / OMARS - SF1) / 
     1                              (1.0 - SF1)
        COWG                      =  1.0 - OM1 * SSGAS1
C
        DTR(I,1,K)                =  EXP( - TAU1 / RMU(I))
        ALAMD                     =  SQRT(3.0 * COW * COWG)
        U                         =  1.50 * COWG / ALAMD
        U2                        =  U + U
        UU                        =  U * U
        EFUN                      =  EXP(- ALAMD * TAU1)
        EFUN2                     =  EFUN * EFUN
        X1                        = (UU - U2 + 1.0) * EFUN2
        RN                        =  1.0 / (UU + U2 + 1.0 - X1)
        X2                        =  ALAMD * RMU(I) 
        Y                         =  1.0 - X2 * X2
        YX                        =  SIGN( MAX( ABS(Y), LOWBND), Y ) 
        DM                        =  OM1 / YX
        GSCW                      =  SSGAS1 * COW
        APPGM                     = (C1(I) + 0.50 + 
     1                               GSCW * (C1(I) + C2(I))) * DM
        APMGM                     = (C1(I) - 0.50 + 
     1                               GSCW * (C1(I) - C2(I))) * DM
        RDF(I,1,K)                = (UU - 1.0) * (1.0  - EFUN2) * RN 
        TDF(I,1,K)                = (U2 + U2) * EFUN * RN
        RDR(I,1,K)                =  APPGM * RDF(I,1,K) + APMGM *
     1                              (TDF(I,1,K) * DTR(I,1,K) - 1.0)
        TDR(I,1,K)                =  APPGM * TDF(I,1,K) +
     1                              (APMGM * RDF(I,1,K) - APPGM + 1.0) *
     2                               DTR(I,1,K)
C
        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)
          RDF(I,3,K)              =  RDF(I,1,K)
          TDF(I,3,K)              =  TDF(I,1,K)
          RDR(I,3,K)              =  RDR(I,1,K)
          TDR(I,3,K)              =  TDR(I,1,K)
          DTR(I,3,K)              =  DTR(I,1,K)
          RDF(I,4,K)              =  RDF(I,1,K)
          TDF(I,4,K)              =  TDF(I,1,K)
          RDR(I,4,K)              =  RDR(I,1,K)
          TDR(I,4,K)              =  TDR(I,1,K)
          DTR(I,4,K)              =  DTR(I,1,K)
        ELSE
          EXTOPT                  =  TAUCS(I,K) + EXTOPT
          OMARCS                  =  TAUOMC(I,K) + TAUR(I,K)
          SSALB                   =  OMARCS / EXTOPT
          SF2                     =  F2(I,K) / OMARCS
          SF                      =  SSALB * SF2
          TAU2                    =  EXTOPT * (1.0 - SF)  
          OM2                     = (SSALB - SF) / (1.0 - SF)
          COW                     =  1.0 - OM2
          SSGAS2                  = (TAUOMGC(I,K) / OMARCS - SF2) / 
     1                              (1.0 - SF2)
          COWG                    =  1.0 - OM2 * SSGAS2
          ALAMD                   =  SQRT(3.0 * COW * COWG)
          U                       =  1.50 * COWG / ALAMD
          U2                      =  U + U
          UU                      =  U * U
          SDTR                    =  EXP(- TAU2 / RMU(I))
          EFUN                    =  EXP(- ALAMD * TAU2)
          EFUN2                   =  EFUN * EFUN
          X3                      = (UU - U2 + 1.0) * EFUN2
          RN                      =  1.0 / (UU + U2 + 1.0 - X3)
          X4                      =  ALAMD * RMU(I)
          Y                       =  1.0 - X4 * X4
          YX                      =  SIGN( MAX( ABS(Y), LOWBND), Y ) 
          DM                      =  OM2 / YX
          GSCW                    =  SSGAS2 * COW
          APPGM                   = (C1(I) + 0.50 +
     1                               GSCW * (C1(I) + C2(I))) * DM
          APMGM                   = (C1(I) - 0.50 +
     1                               GSCW * (C1(I) - C2(I))) * DM
          SRDF                    = (UU - 1.0) * (1.0 - EFUN2) * RN
          STDF                    = (U2 + U2) * EFUN * RN 
          SRDR                    =  APPGM * SRDF + APMGM *
     1                              (STDF * SDTR - 1.0) 
          STDR                    =  APPGM * STDF + (APMGM * SRDF - 
     1                               APPGM + 1.0) * SDTR
          IF (NBLK(I,K) .EQ. 3)                                     THEN
            X5                    =  A1(I,9) * CLD(I,K)
            RDF(I,2,K)            =  RDF(I,1,K) + X5 * 
     1                              (SRDF - RDF(I,1,K))
            TDF(I,2,K)            =  TDF(I,1,K) + X5 * 
     1                              (STDF - TDF(I,1,K))
            RDR(I,2,K)            =  RDR(I,1,K) + X5 * 
     1                              (SRDR - RDR(I,1,K))
            TDR(I,2,K)            =  TDR(I,1,K) + X5 * 
     1                              (STDR - TDR(I,1,K))
            DTR(I,2,K)            =  DTR(I,1,K) + X5 * 
     1                              (SDTR - DTR(I,1,K))
            X6                    =  A1(I,10) * CLD(I,K)
            RDF(I,3,K)            =  RDF(I,1,K) + X6 *
     1                              (SRDF - RDF(I,1,K))
            TDF(I,3,K)            =  TDF(I,1,K) + X6 *
     1                              (STDF - TDF(I,1,K))
            RDR(I,3,K)            =  RDR(I,1,K) + X6 *
     1                              (SRDR - RDR(I,1,K))
            TDR(I,3,K)            =  TDR(I,1,K) + X6 *
     1                              (STDR - TDR(I,1,K))
            DTR(I,3,K)            =  DTR(I,1,K) + X6 *
     1                              (SDTR - DTR(I,1,K))
            X7                    =  A1(I,11) * CLD(I,K)
            RDF(I,4,K)            =  RDF(I,1,K) + X7 *
     1                              (SRDF - RDF(I,1,K))
            TDF(I,4,K)            =  TDF(I,1,K) + X7 *
     1                              (STDF - TDF(I,1,K))
            RDR(I,4,K)            =  RDR(I,1,K) + X7 *
     1                              (SRDR - RDR(I,1,K))
            TDR(I,4,K)            =  TDR(I,1,K) + X7 *
     1                              (STDR - TDR(I,1,K))
            DTR(I,4,K)            =  SDTR
          ELSE IF (NBLK(I,K) .EQ. 1)                                THEN
            RDF(I,4,K)            =  RDF(I,1,K)
            TDF(I,4,K)            =  TDF(I,1,K)
            RDR(I,4,K)            =  RDR(I,1,K)
            TDR(I,4,K)            =  TDR(I,1,K)
            DTR(I,4,K)            =  DTR(I,1,K)
            X                     =  CLD(I,K) / CLDM(I,K)
            RDF(I,2,K)            =  RDF(I,1,K) + X * 
     1                              (SRDF - RDF(I,1,K))
            TDF(I,2,K)            =  TDF(I,1,K) + X *
     1                              (STDF - TDF(I,1,K))
            RDR(I,2,K)            =  RDR(I,1,K) + X *
     1                              (SRDR - RDR(I,1,K))
            TDR(I,2,K)            =  TDR(I,1,K) + X *
     1                              (STDR - TDR(I,1,K))
            DTR(I,2,K)            =  SDTR 
            IF (A1(I,2) .GE. CUT)                                   THEN
              Y                   =  X * A1(I,8)
              RDF(I,3,K)          =  RDF(I,1,K) + Y *
     1                              (SRDF - RDF(I,1,K))
              TDF(I,3,K)          =  TDF(I,1,K) + Y *
     1                              (STDF - TDF(I,1,K))
              RDR(I,3,K)          =  RDR(I,1,K) + Y *
     1                              (SRDR - RDR(I,1,K))
              TDR(I,3,K)          =  TDR(I,1,K) + Y *
     1                              (STDR - TDR(I,1,K))
              DTR(I,3,K)          =  DTR(I,1,K) + Y *
     1                              (SDTR - DTR(I,1,K))
            ELSE
              RDF(I,3,K)          =  RDF(I,1,K)
              TDF(I,3,K)          =  TDF(I,1,K)
              RDR(I,3,K)          =  RDR(I,1,K)
              TDR(I,3,K)          =  TDR(I,1,K)
              DTR(I,3,K)          =  DTR(I,1,K)
            ENDIF
          ELSE IF (NBLK(I,K) .EQ. 2)                                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)
            RDF(I,4,K)            =  RDF(I,1,K)
            TDF(I,4,K)            =  TDF(I,1,K)
            RDR(I,4,K)            =  RDR(I,1,K)
            TDR(I,4,K)            =  TDR(I,1,K)
            DTR(I,4,K)            =  DTR(I,1,K)
            X8                    =  CLD(I,K) / CLDM(I,K)
            RDF(I,3,K)            =  RDF(I,1,K) + X8 *
     1                              (SRDF - RDF(I,1,K))
            TDF(I,3,K)            =  TDF(I,1,K) + X8 *
     1                              (STDF - TDF(I,1,K))
            RDR(I,3,K)            =  RDR(I,1,K) + X8 *
     1                              (SRDR - RDR(I,1,K))
            TDR(I,3,K)            =  TDR(I,1,K) + X8 *
     1                              (STDR - TDR(I,1,K))
            DTR(I,3,K)            =  SDTR
          ENDIF
        ENDIF
  200 CONTINUE
C
      DO 300 I = IL1, IL2
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE FIRST LEVEL (LEV1).                       C
C----------------------------------------------------------------------C
C
        ATRAN0                    =  1.0 - TRAN0(I)
        TMDR(I,1,LEV1)            =  TRAN0(I)
        RMDF(I,1,LEV1)            =  ATRAN0
        CUMDTR(I,1,LEV1)          =  TRAN0(I)
        TMDR(I,2,LEV1)            =  TRAN0(I)
        RMDF(I,2,LEV1)            =  ATRAN0
        CUMDTR(I,2,LEV1)          =  TRAN0(I)
        TMDR(I,3,LEV1)            =  TRAN0(I)
        RMDF(I,3,LEV1)            =  ATRAN0
        CUMDTR(I,3,LEV1)          =  TRAN0(I)
        TMDR(I,4,LEV1)            =  TRAN0(I)
        RMDF(I,4,LEV1)            =  ATRAN0
        CUMDTR(I,4,LEV1)          =  TRAN0(I)
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR THE GROUND LAYER.                             C
C----------------------------------------------------------------------C
C
        RMUR(I,1,LEV)             =  CSALB(I)
        RMUF(I,1,LEV)             =  CSALB(I)
        RMUR(I,2,LEV)             =  ALBSUR(I)
        RMUF(I,2,LEV)             =  ALBSUR(I)
        RMUR(I,3,LEV)             =  ALBSUR(I)
        RMUF(I,3,LEV)             =  ALBSUR(I)
        RMUR(I,4,LEV)             =  ALBSUR(I)
        RMUF(I,4,LEV)             =  ALBSUR(I)
C
        CUMDTR(I,2,LEV)           =  0.0
        CUMDTR(I,3,LEV)           =  0.0
        CUMDTR(I,4,LEV)           =  0.0
  300 CONTINUE
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS DOWNWARD FROM THE SECOND LAYER TO THE SURFACE.    C
C----------------------------------------------------------------------C
C
      DO 450 K = LEV1 + 1, LEV  
        KM1 = K - 1
        L = LEV - K + LEV1
        LP1 = L + 1
        DO 400 I = IL1, IL2
          DMM                     =  TDF(I,1,KM1) / 
     1                              (1.0 - RDF(I,1,KM1) * RMDF(I,1,KM1))
          FMM                     =  RMDF(I,1,KM1) * DMM
          TMDR(I,1,K)             =  CUMDTR(I,1,KM1) * (TDR(I,1,KM1) +
     1                               RDR(I,1,KM1) * FMM) + 
     2                              (TMDR(I,1,KM1) - CUMDTR(I,1,KM1)) * 
     3                               DMM
          RMDF(I,1,K)             =  RDF(I,1,KM1) + TDF(I,1,KM1) * FMM
          CUMDTR(I,1,K)           =  CUMDTR(I,1,KM1) * DTR(I,1,KM1)
C
          IF (A1(I,1) .GE. CUT)                                     THEN
            IF (K .LE. NCT(I))                                      THEN
              TMDR(I,2,K)         =  TMDR(I,1,K)
              RMDF(I,2,K)         =  RMDF(I,1,K)
              CUMDTR(I,2,K)       =  CUMDTR(I,1,K)
            ELSE
              DPP                 =  TDF(I,2,KM1) /
     1                              (1.0 - RMDF(I,2,KM1) * RDF(I,2,KM1))
              FPP                 =  RMDF(I,2,KM1) * DPP
              TMDR(I,2,K)         =  CUMDTR(I,2,KM1) * (TDR(I,2,KM1) +
     1                               RDR(I,2,KM1) * FPP) +
     2                              (TMDR(I,2,KM1) - CUMDTR(I,2,KM1)) *
     3                               DPP
              RMDF(I,2,K)         =  RDF(I,2,KM1) + TDF(I,2,KM1) * FPP
              CUMDTR(I,2,K)       =  CUMDTR(I,2,KM1) * DTR(I,2,KM1)
            ENDIF
          ELSE
            TMDR(I,2,K)           =  1.0
            RMDF(I,2,K)           =  0.0
            CUMDTR(I,2,K)         =  0.0
          ENDIF
C
          IF (A1(I,2) .GE. CUT)                                     THEN
            IF (K .LE. NCT(I))                                      THEN
              TMDR(I,3,K)         =  TMDR(I,1,K)
              RMDF(I,3,K)         =  RMDF(I,1,K)
              CUMDTR(I,3,K)       =  CUMDTR(I,1,K)
            ELSE
              DPP                 =  TDF(I,3,KM1) /
     1                              (1.0 - RMDF(I,3,KM1) * RDF(I,3,KM1))
              FPP                 =  RMDF(I,3,KM1) * DPP
              TMDR(I,3,K)         =  CUMDTR(I,3,KM1) * (TDR(I,3,KM1) +
     1                               RDR(I,3,KM1) * FPP) +
     2                              (TMDR(I,3,KM1) - CUMDTR(I,3,KM1)) *
     3                               DPP
              RMDF(I,3,K)         =  RDF(I,3,KM1) + TDF(I,3,KM1) * FPP
              CUMDTR(I,3,K)       =  CUMDTR(I,3,KM1) * DTR(I,3,KM1)
            ENDIF
C
            IF (A1(I,3) .GE. CUT)                                   THEN
              IF (K .LE. NCT(I))                                    THEN
                TMDR(I,4,K)       =  TMDR(I,1,K)
                RMDF(I,4,K)       =  RMDF(I,1,K)
                CUMDTR(I,4,K)     =  CUMDTR(I,1,K)
              ELSE
                DPP               =  TDF(I,4,KM1) /
     1                              (1.0 - RMDF(I,4,KM1) * RDF(I,4,KM1))
                FPP               =  RMDF(I,4,KM1) * DPP
                TMDR(I,4,K)       =  CUMDTR(I,4,KM1) * (TDR(I,4,KM1) +
     1                               RDR(I,4,KM1) * FPP) +
     2                              (TMDR(I,4,KM1) - CUMDTR(I,4,KM1)) *
     3                               DPP
                RMDF(I,4,K)       =  RDF(I,4,KM1) + TDF(I,4,KM1) * FPP
                CUMDTR(I,4,K)     =  CUMDTR(I,4,KM1) * DTR(I,4,KM1)
              ENDIF
            ELSE
              TMDR(I,4,K)         =  1.0
              RMDF(I,4,K)         =  0.0
              CUMDTR(I,4,K)       =  0.0
            ENDIF
          ELSE
            TMDR(I,3,K)           =  1.0
            RMDF(I,3,K)           =  0.0
            CUMDTR(I,3,K)         =  0.0
            TMDR(I,4,K)           =  1.0
            RMDF(I,4,K)           =  0.0
            CUMDTR(I,4,K)         =  0.0
          ENDIF
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD FROM ONE LAYER ABOVE SURFACE TO THE LEV1.  C
C----------------------------------------------------------------------C
C
          UMM                     =  TDF(I,1,L) /
     1                              (1.0 - RDF(I,1,L) * RMUF(I,1,LP1))
          FMM                     =  RMUF(I,1,LP1) * UMM
          RMUR(I,1,L)             =  RDR(I,1,L) + DTR(I,1,L) *
     1                               RMUR(I,1,LP1) * UMM + (TDR(I,1,L) -
     2                               DTR(I,1,L)) * FMM
          RMUF(I,1,L)             =  RDF(I,1,L) + TDF(I,1,L) * FMM
C
          IF (A1(I,1) .GE. CUT)                                     THEN
            UPP                   =  TDF(I,2,L) /
     1                              (1.0 - RMUF(I,2,LP1) * RDF(I,2,L))
            FPP                   =  RMUF(I,2,LP1) * UPP
            RMUR(I,2,L)           =  RDR(I,2,L) + DTR(I,2,L) *
     1                               RMUR(I,2,LP1) * UPP + (TDR(I,2,L) -
     2                               DTR(I,2,L)) * FPP
            RMUF(I,2,L)           =  RDF(I,2,L) + TDF(I,2,L) * FPP
          ELSE
            RMUR(I,2,L)           =  0.0
            RMUF(I,2,L)           =  0.0
          ENDIF
C
          IF (A1(I,2) .GE. CUT)                                     THEN
            UPP                   =  TDF(I,3,L) /
     1                              (1.0 - RMUF(I,3,LP1) * RDF(I,3,L))
            FPP                   =  RMUF(I,3,LP1) * UPP
            RMUR(I,3,L)           =  RDR(I,3,L) + DTR(I,3,L) *
     1                               RMUR(I,3,LP1) * UPP + (TDR(I,3,L) -
     2                               DTR(I,3,L)) * FPP
            RMUF(I,3,L)           =  RDF(I,3,L) + TDF(I,3,L) * FPP
C
            IF (A1(I,3) .GE. CUT)                                   THEN
              UPP                 =  TDF(I,4,L) /
     1                              (1.0 - RMUF(I,4,LP1) * RDF(I,4,L))
              FPP                 =  RMUF(I,4,LP1) * UPP
              RMUR(I,4,L)         =  RDR(I,4,L) + DTR(I,4,L) *
     1                               RMUR(I,4,LP1) * UPP + (TDR(I,4,L) -
     2                               DTR(I,4,L)) * FPP
              RMUF(I,4,L)         =  RDF(I,4,L) + TDF(I,4,L) * FPP
            ELSE
              RMUR(I,4,L)         =  0.0
              RMUF(I,4,L)         =  0.0
            ENDIF
          ELSE
            RMUR(I,3,L)           =  0.0
            RMUF(I,3,L)           =  0.0
            RMUR(I,4,L)           =  0.0
            RMUF(I,4,L)           =  0.0
          ENDIF
  400   CONTINUE
  450 CONTINUE
C
C----------------------------------------------------------------------C
C     ADD DOWNWARD TO CALCULATE THE RESULTANT REFLECTANCE AND          C
C     TRANSMITTANCE AT FLUX LEVELS.                                    C
C----------------------------------------------------------------------C
C
      DO 550 K = LEV1, LEV
        KM1 = K - 1
        DO 500 I = IL1, IL2
          DMM                     =  1.0 /
     1                              (1.0 - RMUF(I,1,K) * RMDF(I,1,K))
          X                       =  CUMDTR(I,1,K) * RMUR(I,1,K)
          Y                       =  TMDR(I,1,K) - CUMDTR(I,1,K)
          TRAN(I,1,K)             =  CUMDTR(I,1,K) +
     1                              (X * RMDF(I,1,K) + Y) * DMM
          REFL(I,1,K)             = (X + Y * RMUF(I,1,K)) * DMM
C
          IF (A1(I,1) .GE. CUT)                                     THEN
            DPP                   =  1.0 /
     1                              (1.0 - RMUF(I,2,K) * RMDF(I,2,K))
            X                     =  CUMDTR(I,2,K) * RMUR(I,2,K)
            Y                     =  TMDR(I,2,K) - CUMDTR(I,2,K)
            TRAN(I,2,K)           =  A1(I,1) * (CUMDTR(I,2,K) +
     1                              (X * RMDF(I,2,K) + Y) * DPP) +
     2                               A1(I,7) * TRAN(I,1,K)
            REFL(I,2,K)           =  A1(I,1) * (X + Y * RMUF(I,2,K)) *
     1                               DPP + A1(I,7) * REFL(I,1,K)
          ELSE
            TRAN(I,2,K)           =  A1(I,7) * TRAN(I,1,K)
            REFL(I,2,K)           =  A1(I,7) * REFL(I,1,K)
          ENDIF
C
          IF (A1(I,2) .GE. CUT)                                     THEN
            DPP                   =  1.0 /
     1                              (1.0 - RMUF(I,3,K) * RMDF(I,3,K))
            X                     =  CUMDTR(I,3,K) * RMUR(I,3,K)
            Y                     =  TMDR(I,3,K) - CUMDTR(I,3,K)
            TRANPP                =  CUMDTR(I,3,K) +
     1                              (X * RMDF(I,3,K) + Y) * DPP
            REFLPP                = (X + Y * RMUF(I,3,K)) * DPP
            TRAN(I,2,K)           =  A1(I,2) * TRANPP + TRAN(I,2,K)
            REFL(I,2,K)           =  A1(I,2) * REFLPP + REFL(I,2,K)
C
            IF (A1(I,3) .GE. CUT)                                   THEN
              DPP                 =  1.0 /
     1                              (1.0 - RMUF(I,4,K) * RMDF(I,4,K))
              X                   =  CUMDTR(I,4,K) * RMUR(I,4,K)
              Y                   =  TMDR(I,4,K) - CUMDTR(I,4,K)
              TRANPP              =  CUMDTR(I,4,K) +
     1                              (X * RMDF(I,4,K) + Y) * DPP
              REFLPP              = (X + Y * RMUF(I,4,K)) * DPP
              TRAN(I,2,K)         =  A1(I,3) * TRANPP + TRAN(I,2,K)
              REFL(I,2,K)         =  A1(I,3) * REFLPP + REFL(I,2,K)
            ENDIF
          ENDIF
  500   CONTINUE
  550 CONTINUE
C
      RETURN
      END
