      SUBROUTINE LWTRAN4(FU, FD, SLWF, TAUCI, OMCI, GCI, FL, TAUAL,
     1                   TAUG, BF, BS, URBF, DBF, EM0, CLD, CLDM, ANU,
     2                   NCT, NCD, NCU, NCUM, NCDM, LEV1, 
     3                   CUT, MAXC, IL1, IL2, ILG, LAY, LEV)
C
C     * FEB 11,2009 - J.COLE.   NEW VERSION FOR GCM15H:
C     *                         - CORRECT BUG IN SPECIFICATION OF ABSE0.
C     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION LWTRAN3 FOR GCM15G:
C     *                         - S,EMISW,SCATBK,SCATFW,TERM1,S1,T
C     *                           NOW INTERNAL WORK ARRAYS.
C     *                         - SUPPORT ADDED FOR EMISSIVITY<1.  
C     * NOV 22/2006 - M.LAZARE. PREVIOUS VERSION LWTRAN2 FOR GCM15F:
C     *                         - BOUND OF 1.E-10 USED FOR 1-CLD
C     *                           INSTEAD OF SUBTRACTING AN EPSILON.
C     *                         - BOUND OF 1.E-10 USED FOR DENOMINATOR
C     *                           FOR SSALB. 
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 LWTRAN UP TO GCM15E.
C----------------------------------------------------------------------C
C     CALCULATION OF LONGWAVE RADIATIVE TRANSFER USING ABSORPTION      C
C     APPROXIMATION. THE FINITE CLOUD EFFECT IS PROPERLY CONSIDERED    C
C     WITH RANDOM AND FULL OVERLAP ASSUMPTION. CLOUD SUBGRID           C
C     VARIABILITY IS INCLUDED (BASED ON LI, 2002 JAS P3302; LI AND     C
C     BARKER JAS P3321).                                               C
C                                                                      C
C     FU:     UPWARD INFRARED FLUX                                     C
C     FD:     DOWNWARD INFRARED FLUX                                   C
C     SLWF:   INPUT SOLAR FLUX AT MODEL TOP LEVEL FOR EACH BAND        C
C     TAUCI:  CLOUD OPTICAL DEPTH FOR THE INFRARED                     C
C     OMCI:   CLOUD SINGLE SCATTERING ALBEDO TIMES OPTICAL DEPTH       C
C     GCI:    CLOUD ASYMMETRY FACTOR TIMES OMCI                        C
C     FL:     SQUARE OF CLOUD ASYMMETRY FACTOR                         C
C     TAUAL:  AEROSOL OPTICAL DEPTH FOR THE INFRARED                   C
C     TAUG:   GASEOUS OPTICAL DEPTH FOR THE INFRARED                   C
C     BF:     BLACKBODY INTENSITY INTEGRATED OVER EACH BAND AT EACH    C
C             LEVEL IN UNITS W / M^2 / SR. THEREFOR A PI FACTOR NEEDED C
C             FOR FLUX                                                 C
C     BS:     THE BLACKBODY INTENSITY AT THE SURFACE.                  C
C     URBF:   U TIMES THE DIFFERENCE OF LOG(BF) FOR TWO NEIGHBOR       C
C             LEVELS USED FOR EXPONENTIAL SOURCE FUNCTION              C
C     DBF:    DIFFERENCE OF BF FOR TWO NEIGHBOR LEVELS USED FOR        C
C             LINEAR SOURCE FUNCTION                                   C
C     EM0:    SURFACE EMISSION                                         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 LI AND       C
C             BARKER JAS P3321).                                       C
C     ANU:    NU FACTOR FOR CLOUD SUBGRID VARIABILITY                  C
C     S:      TOTAL SCATTERING                                         C
C     SCATBK: BACKWARD SCATTERING                                      C
C     SCATFW: FORWARD SCATTERING                                       C
C     SCATSM: INTERNAL SCATTERING                                      C
C     TAUM:   TAUM(1) A FACTOR RELATED TO TAU IN ZETA FACTOR FOR       C
C             LINEAR SOURCE TERM; TAUM(2) THE CUMULATED TAUM(1) FOR    C
C             SUBGRID VARIABILITY CALCULATION                          C
C     XU:     THE EMISSION PART IN THE UPWARD FLUX TRANSMISSION        C
C             (LI, 2002 JAS P3302)                                     C
C     XD:     THE EMISSION PART IN THE DOWNWARD FLUX TRANSMISSION      C
C     DTR:    DIRECT TRANSMISSION                                      C
C     FY:     UPWARD FLUX FOR PURE CLEAR PORTION (1) AND PURE CLOUD    C
C             PORTION (2)                                              C
C     FX:     THE SAME AS FY BUT FOR THE DOWNWARD FLUX                 C
C     FW:     A TERM FOR TRANSFER WITHIN CLDM                          C
C     NCT:    THE HIGHEST CLOUD TOP LEVEL FOR THE LONGITUDE AND        C
C             LATITUDE LOOP (ILG)                                      C
C     NCD:    LAYER INSIDE A CLOUD BLOCK ACCOUNTED FROM CLOUD TOP      C
C     NCU:    LAYER INSIDE A CLOUD BLOCK ACCOUNTED FROM CLOUD BOTTOM   C
C     NCUM:   MAXIMUM LOOP NUMBER CLOUD VERTICAL CORRELATION ACCOUNTED C
C             FROM LOWER LEVEL TO HIGHER LEVEL                         C
C     NCDM:   MAXIMUM LOOP NUMBER CLOUD VERTICAL CORRELATION ACCOUNTED C
C             FROM HIGHER LEVEL TO LOWER LEVEL                         C
C----------------------------------------------------------------------C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL   FU(ILG,2,LEV), FD(ILG,2,LEV)
      REAL   SLWF(ILG),TAUCI(ILG,LAY),OMCI(ILG,LAY),GCI(ILG,LAY),
     1       FL(ILG,LAY),TAUAL(ILG,LAY),TAUG(ILG,LAY),BF(ILG,LEV), 
     2       BS(ILG),URBF(ILG,LAY),DBF(ILG,LAY),EM0(ILG),CLD(ILG,LAY),
     3       CLDM(ILG,LAY),ANU(ILG,LAY)
      REAL   SCATSM(ILG,4,LAY), TAUM(ILG,4,LAY),
     2       XD(ILG,4,LAY), XU(ILG,4,LAY), DTR(ILG,4,LAY), 
     3       FX(ILG,4,LEV), FY(ILG,4,LEV), FW(ILG,4,LEV)
C
C     * INTERNAL WORK ARRAYS.
C
      REAL  S(ILG,LAY), EMISW(ILG,LAY), SCATBK(ILG,LAY), SCATFW(ILG,LAY)
      REAL  TERM1(ILG), S1(ILG), T(ILG)
      INTEGER NCT(ILG), NCD(ILG,LAY), NCU(ILG,LAY), NCUM(LAY), NCDM(LAY)
C
      DATA  RU / 1.6487213 /
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR FIRST LAYER. CALCULATE THE DOWNWARD FLUX IN   C
C     THE SECOND LAYER                                                 C
C     COMBINE THE OPTICAL PROPERTIES FOR THE INFRARED,                 C
C     1, AEROSOL + GAS; 2, CLOUD + AEROSOL + GAS.                      C
C     FD (FU) IS DOWN (UPWARD) FLUX, FX (FY) IS THE INCIDENT FLUX      C
C     ABOVE (BELOW) THE CONSIDERED LAYER.                              C
C     GAUSSIAN INTEGRATION AND DIFFUSIVITY FACTOR, RU (LI JAS 2000)    C
C     ABOVE MAXC, EXPONENTIAL SOURCE FUNCTION IS USED                  C
C     BELOW MAXC, LINEAR SOURCE FUNCTION IS USED                       C
C----------------------------------------------------------------------C
C
      L1 =  LEV1
      L2 =  LEV1 + 1
      DO 100 I = IL1, IL2
        FD(I,1,LEV1)              =  SLWF(I)
        FD(I,2,LEV1)              =  SLWF(I)
        FX(I,1,LEV1)              =  FD(I,1,LEV1)
        FX(I,2,LEV1)              =  FD(I,2,LEV1)
  100 CONTINUE
C
      DO 150 K = L2, MAXC
        KM1 = K - 1
        DO 125 I = IL1, IL2
          TAUL1                   =  TAUAL(I,KM1) + TAUG(I,KM1)
          RTAUL1                  =  TAUL1 * RU
          DTR(I,1,KM1)            =  EXP (- RTAUL1)
          UBETA                   =  URBF(I,KM1) / (TAUL1 + 1.E-20)
          EPSD                    =  UBETA + 1.0
          EPSU                    =  UBETA - 1.0
C
          IF (ABS(EPSD) .GT. 0.001)                                 THEN
            XD(I,1,KM1)           = (BF(I,K) - BF(I,KM1) * 
     1                               DTR(I,1,KM1)) / EPSD
          ELSE
            XD(I,1,KM1)           =  RTAUL1 * BF(I,KM1) * DTR(I,1,KM1)
          ENDIF
C
          IF (ABS(EPSU) .GT. 0.001)                                 THEN
            XU(I,1,KM1)           = (BF(I,K) * DTR(I,1,KM1) - 
     1                               BF(I,KM1)) / EPSU
          ELSE
            XU(I,1,KM1)           =  RTAUL1 * BF(I,K) * DTR(I,1,KM1)
          ENDIF
C
          FD(I,1,K)               =  FD(I,1,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
          FD(I,2,K)               =  FD(I,2,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
          FX(I,1,K)               =  FD(I,2,K)
          FX(I,2,K)               =  FD(I,2,K)
  125   CONTINUE
  150 CONTINUE 
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS DOWNWARD FROM THE SECOND LAYER TO THE SURFACE.    C
C     DETERMINE THE XU FOR THE UPWARD PATH.                            C
C     USING EXPONENTIAL SOURCE FUNCTION FOR CLR FLUX CALCULATION AND   C
C     ALSO FOR ALL SKY FLUX IN CLOUD FREE LAYERS.                      C
C----------------------------------------------------------------------C
C
      IF (MAXC .LT. LEV)                                            THEN
        DO 250 K = MAXC + 1, LEV
          KM1 = K - 1
          KM2 = KM1 - 1
          DO 225 I = IL1, IL2
            TAUL1                 =  TAUAL(I,KM1) + TAUG(I,KM1)
            RTAUL1                =  TAUL1 * RU
            DTR(I,1,KM1)          =  EXP (- RTAUL1)
C
            UBETA                 =  URBF(I,KM1) / (TAUL1 + 1.E-20)
            EPSD                  =  UBETA + 1.0
            EPSU                  =  UBETA - 1.0
C
            IF (ABS(EPSD) .GT. 0.001)                               THEN
              XD(I,1,KM1)         = (BF(I,K) - BF(I,KM1) * 
     1                               DTR(I,1,KM1)) / EPSD
            ELSE
              XD(I,1,KM1)         =  RTAUL1 * BF(I,KM1) * DTR(I,1,KM1)
            ENDIF
C
            IF (ABS(EPSU) .GT. 0.001)                               THEN
              XU(I,1,KM1)         = (BF(I,K) * DTR(I,1,KM1) - 
     1                               BF(I,KM1)) / EPSU
            ELSE
              XU(I,1,KM1)         =  RTAUL1 * BF(I,K) * DTR(I,1,KM1)
            ENDIF
C
            FD(I,1,K)             =  FD(I,1,KM1) * DTR(I,1,KM1) +
     1                               XD(I,1,KM1)
            IF (CLD(I,KM1) .LT. CUT)                                THEN
              FD(I,2,K)           =  FD(I,2,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
              FX(I,1,K)           =  FD(I,2,K)
              FX(I,2,K)           =  FD(I,2,K)
            ELSE
              TAUL2               =  TAUCI(I,KM1) + TAUG(I,KM1)
              SSALB               =  OMCI(I,KM1) / (TAUL2 + 1.E-10)
              SF                  =  SSALB * FL(I,KM1)
              W                   = (SSALB - SF) / (1.0 - SF)
              COW                 =  1.0 - W
              TAUM(I,1,KM1)       = (COW * TAUL2 * (1.0 - SF) + 
     1                               TAUAL(I,KM1)) * RU
              ZETA                =  DBF(I,KM1) / TAUM(I,1,KM1)
              TAU2                =  TAUM(I,1,KM1) + TAUM(I,1,KM1)
C
              SANU                =  ANU(I,KM1)
              X                   =  SANU / (SANU + TAUM(I,1,KM1))
              Y                   =  SANU / (SANU + TAU2)
C
              IF (SANU .LE. 0.50)                                   THEN
                DTR(I,2,KM1)      =  SQRT(X)
                DTR2              =  SQRT(Y)
                EMISW(I,KM1)      =  ZETA * (SQRT(1.0 + TAU2) - 1.0)
                EMBK              =  1.0 - SQRT(1.0 + TAU2 + TAU2)
              ELSE IF (SANU .GT. 0.50 .AND. SANU .LE. 1.0)          THEN
                WT                =  2.0 * SANU - 1.0
                SX                =  SQRT(X)
                SY                =  SQRT(Y)
                DTR(I,2,KM1)      =  SX + (X - SX) * WT
                DTR2              =  SY + (Y - SY) * WT
                P1                =  SQRT(1.0 + TAU2) - 1.0
                EMISW(I,KM1)      =  ZETA * (P1 + (LOG(1.0 + 
     1                               TAUM(I,1,KM1)) - P1) * WT)
                P2                =  1.0 - SQRT(1.0 + TAU2 + TAU2)
                EMBK              =  P2 - (LOG(1.0 + TAU2) + P2) * WT
              ELSE IF (SANU .GT. 1.0 .AND. SANU .LE. 2.0)           THEN
                WT                =  SANU - 1.0
                DTR(I,2,KM1)      =  X + (X * X - X) * WT
                DTR2              =  Y + (Y * Y - Y) * WT
                Z                 =  SANU / (SANU - 1.0)
                P1                =  LOG(1.0 + TAUM(I,1,KM1))
                EMISW(I,KM1)      =  ZETA * (P1 + (Z * (1.0 - X) - P1) *
     1                               WT)
                P2                =  - LOG(1.0 + TAU2)
                EMBK              =  P2 + (Z * (Y - 1.0) - P2) * WT
              ELSE IF (SANU .GT. 2.0 .AND. SANU .LE. 3.0)           THEN
                X2                =  X * X
                Y2                =  Y * Y
                WT                =  SANU - 2.0
                DTR(I,2,KM1)      =  X2 + (X * X2 - X2) * WT
                DTR2              =  Y2 + (Y * Y2 - Y2) * WT
                Z                 =  SANU / (SANU - 1.0)
                EMISW(I,KM1)      =  Z * ZETA *
     1                              (1.0 - X + (X - X2) * WT)
                EMBK              =  Z * (Y - 1.0 + (Y2 - Y) * WT)
              ELSE IF (SANU .GT. 3.0 .AND. SANU .LE. 4.0)           THEN
                X2                =  X * X
                Y2                =  Y * Y
                X3                =  X2 * X
                Y3                =  Y2 * Y
                WT                =  SANU - 3.0
                DTR(I,2,KM1)      =  X3 + (X2 * X2 - X3) * WT
                DTR2              =  Y3 + (Y2 * Y2 - Y3) * WT
                Z                 =  SANU / (SANU - 1.0)
                EMISW(I,KM1)      =  Z * ZETA *
     1                              (1.0 - X2 + (X2 - X3) * WT)
                EMBK              =  Z * (Y2 - 1.0 + (Y3 - Y2) * WT)
C
C----------------------------------------------------------------------C
C     FOR ANU > 4, THE INHOMOENEITY EFFECT IS VERY WEAK, FOR SAVING    C
C     THE INTEGER ANU IS ASSUMED. FOR ANU > 20, HOMOGENOUS IS ASSUMED  C 
C----------------------------------------------------------------------C
C
              ELSE IF (SANU .GT. 4.0 .AND. SANU .LE. 20.0)          THEN
                NANU              =  NINT(SANU)
                DTR(I,2,KM1)      =  X ** NANU
                DTR2              =  Y ** NANU
                Z                 =  SANU / (SANU - 1.0)
                EMISW(I,KM1)      =  Z * ZETA * (1.0 - DTR(I,2,KM1) / X)
                EMBK              =  Z * (DTR2 / Y - 1.0)
              ELSE
                DTR(I,2,KM1)      =  EXP(- TAUM(I,1,KM1))
                DTR2              =  EXP(- TAU2)
                Z                 =  SANU / (SANU - 1.0)
                EMISW(I,KM1)      =  Z * ZETA * (1.0 - DTR(I,2,KM1) / X)
                EMBK              =  Z * (DTR2 / Y - 1.0) 
              ENDIF
C
              XD(I,2,KM1)         =  BF(I,K) - BF(I,KM1) * 
     1                               DTR(I,2,KM1) - EMISW(I,KM1) 
              XU(I,2,KM1)         =  BF(I,KM1) - BF(I,K) * 
     1                               DTR(I,2,KM1) + EMISW(I,KM1)
C
              WGRCOW              =  W * GCI(I,KM1) / COW
              TAUDTR              =  TAUM(I,1,KM1) * DTR(I,2,KM1)
C
              SCATFW(I,KM1)       =  WGRCOW * X * TAUDTR
              SCATBK(I,KM1)       =  0.5 * WGRCOW * (DTR2 - 1.0)
C
              X                   =  WGRCOW * (2.0 * EMISW(I,KM1) +
     1                              (0.5 * EMBK - TAUDTR) * ZETA)
              SCATSM(I,1,KM1)     =  - SCATBK(I,KM1) * BF(I,K) -
     1                               SCATFW(I,KM1) * BF(I,KM1) - X
              SCATSM(I,2,KM1)     =  - SCATBK(I,KM1) * BF(I,KM1) -
     1                               SCATFW(I,KM1) * BF(I,K) + X
C
              IF (K .EQ. L2)                                        THEN
                FX(I,1,K)         =  FX(I,1,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
                FX(I,2,K)         =  FX(I,2,KM1) * DTR(I,2,KM1) + 
     1                               XD(I,2,KM1)
              ELSE IF (CLD(I,KM1) .LE. CLD(I,KM2))                  THEN
                FX(I,1,K)         = ( FX(I,2,KM1) + (1.0 - CLD(I,KM2)) /
     1                              (MAX(1.0 - CLD(I,KM1),1.E-10)) * 
     2                              (FX(I,1,KM1) - FX(I,2,KM1)) ) * 
     3                               DTR(I,1,KM1) + XD(I,1,KM1)
                FX(I,2,K)         =  FX(I,2,KM1) * DTR(I,2,KM1) + 
     1                               XD(I,2,KM1)
              ELSE IF (CLD(I,KM1) .GT. CLD(I,KM2))                  THEN
                FX(I,1,K)         =  FX(I,1,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
                FX(I,2,K)         = (FX(I,1,KM1) + 
     1                               CLD(I,KM2) / CLD(I,KM1) * 
     2                              (FX(I,2,KM1) - FX(I,1,KM1))) * 
     3                               DTR(I,2,KM1) + XD(I,2,KM1)
              ENDIF
C
              FD(I,2,K)           =  FX(I,1,K) + CLD(I,KM1) * 
     1                              (FX(I,2,K) - FX(I,1,K))
            ENDIF   
  225     CONTINUE
  250   CONTINUE
      ENDIF
C
C----------------------------------------------------------------------C
C     INITIALIZATION FOR SURFACE                                       C
C----------------------------------------------------------------------C
C
      K = LEV - 1
      DO 300 I = IL1, IL2
        EMBS                      =  EM0(I) * BS(I)
        ABSE0                     =  1.0 - EM0(I)
        FU(I,1,LEV)               =  EMBS + ABSE0 * FD(I,1,LEV)
        FY(I,1,LEV)               =  EMBS + ABSE0 * FX(I,1,LEV)
        FY(I,2,LEV)               =  EMBS + ABSE0 * FX(I,2,LEV)
        FU(I,2,LEV)               =  FY(I,1,LEV) + CLD(I,K) *
     1                              (FY(I,2,LEV) - FY(I,1,LEV))
        FW(I,2,LEV)               =  FY(I,2,LEV)
C
C----------------------------------------------------------------------C
C     DETERMINING THE UPWARD FLUX FOR THE FIRST LAY ABOVE SURFACE      C
C----------------------------------------------------------------------C
C
        FU(I,1,K)                 =  FU(I,1,LEV) * DTR(I,1,K) +
     1                               XU(I,1,K) 
        FY(I,1,K)                 =  FY(I,1,LEV) * DTR(I,1,K) +
     1                               XU(I,1,K)
C
        IF (CLD(I,K) .LT. CUT)                                      THEN
          S(I,K)                  =  0.0
          FY(I,2,K)               =  FY(I,2,LEV) * DTR(I,1,K) +
     1                               XU(I,1,K)
          FU(I,2,K)               =  FY(I,1,K)
          FW(I,2,K)               =  FY(I,2,K)
          TAUM(I,2,K)             =  0.0
        ELSE
          S(I,K)                  =  SCATBK(I,K) * FX(I,2,K) +
     1                               SCATFW(I,K) * FY(I,2,LEV) +
     2                               SCATSM(I,2,K)
C
          FY(I,2,K)               =  FY(I,2,LEV) * DTR(I,2,K) +
     1                               XU(I,2,K) + S(I,K)                
          FU(I,2,K)               =  FY(I,1,K) + CLD(I,K) * 
     1                              (FY(I,2,K) - FY(I,1,K))
          FW(I,2,K)               =  FY(I,2,K)
          TAUM(I,2,K)             =  TAUM(I,1,K)
        ENDIF
  300 CONTINUE
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD FROM THE SECOND LAYER TO MAXC              C
C     SCATTERING EFFECT FOR UPWARD PATH IS INCLUDED                    C
C----------------------------------------------------------------------C
C
      DO 450 K = LEV - 2, MAXC, - 1
        KP1 = K + 1
        KP2 = K + 2
        DO 400 I = IL1, IL2
          IF (K .GE. NCT(I))                                        THEN
            FU(I,1,K)             =  FU(I,1,KP1) * DTR(I,1,K) + 
     1                               XU(I,1,K)
C
            IF (CLD(I,K) .LT. CUT)                                  THEN
              FU(I,2,K)           =  FU(I,2,KP1) * DTR(I,1,K) + 
     1                               XU(I,1,K)
              FY(I,1,K)           =  FU(I,2,K)
              FY(I,2,K)           =  FU(I,2,K)
              FW(I,2,K)           =  FU(I,2,K)
              TAUM(I,2,K)         =  0.0
            ELSE
C
C----------------------------------------------------------------------C
C     FY(I,2,K) CONTAINS UNPERTURBED + BACKWARD SCATTERING EFFECT +    C
C     FORWARD SCATTERING EFFECT + INTERNAL SCATTERING EFFECT           C
C    (LI AND FU, JAS 2000)                                             C
C----------------------------------------------------------------------C
C
              IF (CLD(I,K) .LE. CLD(I,KP1) .OR.
     1                          CLD(I,K) - CLDM(I,K) .LT. CUT)      THEN
C
                FY(I,1,K)         = ( FY(I,2,KP1) + (1.0 - CLD(I,KP1)) /
     1                              (MAX(1.0 - CLD(I,K),1.E-10)) *
     2                              (FY(I,1,KP1) - FY(I,2,KP1)) ) * 
     3                               DTR(I,1,K) + XU(I,1,K)
                T(I)              =  FY(I,2,KP1)
              ELSE
                FY(I,1,K)         =  FY(I,1,KP1) * DTR(I,1,K) + 
     1                               XU(I,1,K)
                T(I)              =  FY(I,1,KP1) + 
     1                              (CLD(I,KP1) - CLDM(I,KP1)) / 
     2                              (CLD(I,K) - CLDM(I,K)) *
     3                              (FY(I,2,KP1) - FY(I,1,KP1))
              ENDIF
C
              BKINS               =  SCATBK(I,K) * FX(I,2,K) + 
     1                               SCATSM(I,2,K)
              FY(I,2,K)           =  T(I) * (DTR(I,2,K) + SCATFW(I,K)) +
     1                               XU(I,2,K) + BKINS 
              TAUM(I,2,K)         =  TAUM(I,2,KP1) + TAUM(I,1,K)
              S1(I)               =  0.0
              S(I,K)              =  BKINS + SCATFW(I,K) * FY(I,2,KP1) 
              TERM1(I)            =  0.0 
C
              IF (NCU(I,K) .GT. 1)                                  THEN
                KX = K + NCU(I,K)
                KXM = KX - 1
C
                SANU              =  ANU(I,KXM)
                ANUTAU            =  SANU / (SANU + TAUM(I,2,K))
                IF (SANU .LE. 0.50)                                 THEN
                  DTRGW           =  SQRT(ANUTAU)
                ELSE IF (SANU .GT. 0.50 .AND. SANU .LE. 1.0)        THEN
                  X               =  SQRT(ANUTAU)
                  DTRGW           =  X + 2.0 * (SANU - 0.50) *
     1                              (ANUTAU - X)
                ELSE IF (SANU .GT. 1.0 .AND. SANU .LE. 2.0)         THEN
                  DTRGW           =  ANUTAU + (SANU - 1.0) * ANUTAU *
     1                              (ANUTAU - 1.0)
                ELSE IF (SANU .GT. 2.0 .AND. SANU .LE. 3.0)         THEN
                  X               =  ANUTAU * ANUTAU
                  DTRGW           =  X + (SANU - 2.0) * X *
     1                              (ANUTAU - 1.0)
                ELSE IF (SANU .GT. 3.0 .AND. SANU .LE. 4.0)         THEN
                  X               =  ANUTAU * ANUTAU * ANUTAU
                  DTRGW           =  X + (SANU - 3.0) * X *
     1                              (ANUTAU - 1.0) 
                ELSE IF (SANU .GT. 4.0 .AND. SANU .LE. 20.0)        THEN
                  DTRGW           =  ANUTAU ** (NINT(SANU))
                ELSE
                  DTRGW           =  EXP(- TAUM(I,2,K))
                ENDIF
C
                TERM1(I)          = (FW(I,2,KX) - BF(I,KX)) * DTRGW
                S1(I)             = (EMISW(I,KP1) + S(I,KP1)) * 
     1                               DTR(I,2,K)
              ENDIF
            ENDIF
          ENDIF 
  400   CONTINUE
C
C----------------------------------------------------------------------C
C     DETERMINING THE TERMS GOING INTO THE CORRELATION CALCULATIONS    C
C     FOR SUBGRID VARIABILITY FOR CLDM PORTION.                        C
C----------------------------------------------------------------------C
C
        IF (NCUM(K) .GT. 2)                                         THEN
          DO 420 KK = KP2, K + NCUM(K) - 1
          DO 420 I = IL1, IL2
            IF (K .GE. NCT(I) .AND. CLD(I,K) .GE. CUT .AND.
     1          NCU(I,K) .GT. 2 .AND. KK .LE. K + NCU(I,K) - 1)     THEN
C
              SANU                =  ANU(I,KK)
              ANUTAU              =  SANU / (SANU + 
     1                               TAUM(I,2,K) - TAUM(I,2,KK))
              IF (SANU .LE. 0.50)                                   THEN
                DTRGW             =  SQRT(ANUTAU)
              ELSE IF (SANU .GT. 0.50 .AND. SANU .LE. 1.0)          THEN
                X                 =  SQRT(ANUTAU)
                DTRGW             =  X + 2.0 * (SANU - 0.50) *
     1                              (ANUTAU - X)
              ELSE IF (SANU .GT. 1.0 .AND. SANU .LE. 2.0)           THEN
                DTRGW             =  ANUTAU + (SANU - 1.0) * ANUTAU *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 2.0 .AND. SANU .LE. 3.0)           THEN
                X                 =  ANUTAU * ANUTAU
                DTRGW             =  X + (SANU - 2.0) * X *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 3.0 .AND. SANU .LE. 4.0)           THEN
                X                 =  ANUTAU * ANUTAU * ANUTAU
                DTRGW             =  X + (SANU - 3.0) * X *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 4.0 .AND. SANU .LE. 20.0)          THEN
                DTRGW             =  ANUTAU ** (NINT(SANU))
              ELSE
                DTRGW             =  EXP(-(TAUM(I,2,K) - TAUM(I,2,KK)))
              ENDIF
C
              S1(I)               =  S1(I) + 
     1                              (EMISW(I,KK) + S(I,KK)) * DTRGW
            ENDIF
  420     CONTINUE
        ENDIF
C
C----------------------------------------------------------------------C
C     IN CLDM REGION CONSIDER THE CORRELATION BETWEEN DIFFERENT LAYERS C
C----------------------------------------------------------------------C
C
        DO 430 I = IL1, IL2
          IF (K .GE. NCT(I))                                        THEN
            IF (CLD(I,K) .GE. CUT)                                  THEN
              IF (NCU(I,K) .EQ. 1)                                  THEN
                FW(I,2,K)         =  FY(I,2,K) 
                FU(I,2,K)         =  FY(I,1,K) + CLD(I,K) * (FY(I,2,K) -
     1                               FY(I,1,K)) 
              ELSE
                FW(I,2,K)         =  TERM1(I) + S1(I) + BF(I,K) + 
     1                               EMISW(I,K) + S(I,K)
                FU(I,2,K)         =  CLDM(I,K) * (FW(I,2,K) - 
     1                               FY(I,2,K)) + FY(I,1,K) + 
     2                               CLD(I,K) * (FY(I,2,K) - FY(I,1,K))
              ENDIF
            ENDIF
          ENDIF
  430   CONTINUE        
  450 CONTINUE
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD ABOVE THE HIGHEST CLOUD  TO THE TOA, NO    C
C     SCATTERING                                                       C
C----------------------------------------------------------------------C
C
      DO 550 K = LEV - 1, L1, - 1
        KP1 = K + 1
C
        DO 500 I = IL1, IL2
          IF (KP1 .LE. NCT(I))                                      THEN
            FU(I,1,K)             =  FU(I,1,KP1) * DTR(I,1,K) + 
     1                               XU(I,1,K)
            FU(I,2,K)             =  FU(I,2,KP1) * DTR(I,1,K) + 
     1                               XU(I,1,K)
          ENDIF
C
C----------------------------------------------------------------------C
C     SCATTERING EFFECT FOR DOWNWARD PATH AT THE TOP LAYER OF THE      C
C     HIGHEST CLOUD                                                    C
C----------------------------------------------------------------------C
C
          IF (K .EQ. NCT(I))                                        THEN
            FW(I,1,K)             =  FX(I,1,K)
            FWINS                 =  SCATSM(I,1,K) + 
     1                               SCATFW(I,K) * FX(I,2,K)
            FMBK                  =  FX(I,2,K) * DTR(I,2,K) + 
     1                               XD(I,2,K) + FWINS
            FX(I,2,KP1)           =  FMBK + SCATBK(I,K) * FY(I,2,KP1)
            TAUM(I,2,K)           =  TAUM(I,1,K)
            S(I,K)                =  SCATBK(I,K) * FW(I,2,KP1) + FWINS
C
            FW(I,1,KP1)           =  FMBK + SCATBK(I,K) * FW(I,2,KP1)
            FD(I,2,KP1)           =  FX(I,1,KP1) + CLD(I,K) *
     1                              (FX(I,2,KP1) - FX(I,1,KP1))
          ENDIF
  500   CONTINUE
  550 CONTINUE
C
C----------------------------------------------------------------------C
C     SCATTERING EFFECT FOR DOWNWARD PATH IN FROM MAXC TO THE SURFACE  C
C----------------------------------------------------------------------C
C
      DO 750 K = MAXC + 2, LEV
        KM1 = K - 1
        KM2 = K - 2
        KM3 = K - 3  
        DO 700 I = IL1, IL2
          IF (KM2 .GE. NCT(I))                                      THEN
            IF (CLD(I,KM1) .LT. CUT)                                THEN
              FD(I,2,K)           =  FD(I,2,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
              FX(I,1,K)           =  FD(I,2,K)
              FX(I,2,K)           =  FD(I,2,K)
              FW(I,1,K)           =  FD(I,2,K)
              TAUM(I,2,KM1)       =  0.0
            ELSE 
              IF (CLD(I,KM1) .LE. CLD(I,KM2) .OR. 
     1                        CLD(I,KM1) - CLDM(I,KM1) .LT. CUT)    THEN
C
                FX(I,1,K)         = (FX(I,2,KM1) + (1.0 - CLD(I,KM2)) /
     1                              (MAX(1.0 - CLD(I,KM1),1.E-10)) * 
     2                              (FX(I,1,KM1) - FX(I,2,KM1))) * 
     3                               DTR(I,1,KM1) + XD(I,1,KM1)
                T(I)              =  FX(I,2,KM1)
              ELSE
                FX(I,1,K)         =  FX(I,1,KM1) * DTR(I,1,KM1) + 
     1                               XD(I,1,KM1)
                T(I)              =  FX(I,1,KM1) + 
     1                              (CLD(I,KM2) - CLDM(I,KM2)) / 
     2                              (CLD(I,KM1) - CLDM(I,KM1)) * 
     3                              (FX(I,2,KM1) -  FX(I,1,KM1))
              ENDIF
C
              FX(I,2,K)           =  T(I) * DTR(I,2,KM1) + XD(I,2,KM1) +
     1                               SCATBK(I,KM1) * FY(I,2,K) + 
     2                               SCATFW(I,KM1) * T(I) + 
     3                               SCATSM(I,1,KM1)
C
              TAUM(I,2,KM1)       =  TAUM(I,2,KM2) + TAUM(I,1,KM1)
              S1(I)               =  0.0
              S(I,KM1)            =  SCATBK(I,KM1) * FW(I,2,K) + 
     1                               SCATFW(I,KM1) * FW(I,1,KM1) + 
     2                               SCATSM(I,1,KM1)
              TERM1(I)            =  0.0 
C
              IF (NCD(I,KM1) .GT. 1)                                THEN
                KX = K - NCD(I,KM1)
                SANU              =  ANU(I,KX)
                ANUTAU            =  SANU / (SANU + TAUM(I,2,KM1))
                IF (SANU .LE. 0.50)                                 THEN
                  DTRGW           =  SQRT(ANUTAU)
                ELSE IF (SANU .GT. 0.50 .AND. SANU .LE. 1.0)        THEN
                  X               =  SQRT(ANUTAU)
                  DTRGW           =  X + 2.0 * (SANU - 0.50) *
     1                              (ANUTAU - X)
                ELSE IF (SANU .GT. 1.0 .AND. SANU .LE. 2.0)         THEN
                  DTRGW           =  ANUTAU + (SANU - 1.0) * ANUTAU *
     1                              (ANUTAU - 1.0)
                ELSE IF (SANU .GT. 2.0 .AND. SANU .LE. 3.0)         THEN
                  X               =  ANUTAU * ANUTAU
                  DTRGW           =  X + (SANU - 2.0) * X *
     1                              (ANUTAU - 1.0)
                ELSE IF (SANU .GT. 3.0 .AND. SANU .LE. 4.0)         THEN
                  X               =  ANUTAU * ANUTAU * ANUTAU
                  DTRGW           =  X + (SANU - 3.0) * X *
     1                              (ANUTAU - 1.0)
                ELSE IF (SANU .GT. 4.0 .AND. SANU .LE. 20.0)        THEN
                  DTRGW           =  ANUTAU ** (NINT(SANU))
                ELSE
                  DTRGW           =  EXP(- TAUM(I,2,KM1))
                ENDIF
C
                TERM1(I)          = (FW(I,1,KX) - BF(I,KX)) * DTRGW
                S1(I)             = (S(I,KM2) - EMISW(I,KM2)) * 
     1                               DTR(I,2,KM1)
              ENDIF
            ENDIF
          ENDIF
  700   CONTINUE          
C
C----------------------------------------------------------------------C
C     DETERMINING THE TERMS GOING INTO THE CORRELATION CALCULATIONS    C
C     FOR CLDM PORTION.                                                C
C----------------------------------------------------------------------C
C
        IF (NCDM(KM1) .GT. 2)                                       THEN
C
C----------------------------------------------------------------------C
C     NOTE THAT IN THE FOLLOWING LOOP, "KM1" IS ACTUALLY THE           C
C     REPRESENTATIVE VARIABLE, SO THAT K-NCD(I,KM1) IS ACTUALLY        C
C     KM1-NCD(I,KM1)+1. THE SIMPLER FORM IS USED ONLY FOR              C
C     COMPUTATIONAL EFFICIENCY.                                        C
C----------------------------------------------------------------------C
C
          DO 720 KK = KM3, K - NCDM(KM1), - 1
          DO 720 I = IL1, IL2
            IF (KM2 .GE. NCT(I) .AND. CLD(I,KM1) .GE. CUT .AND.
     1          NCD(I,KM1) .GT. 2 .AND. KK .GE. K - NCD(I,KM1))     THEN
C              
              SANU                =  ANU(I,KK)
              ANUTAU              =  SANU / (SANU + 
     1                               TAUM(I,2,KM1) - TAUM(I,2,KK))
              IF (SANU .LE. 0.50)                                   THEN
                DTRGW             =  SQRT(ANUTAU)
              ELSE IF (SANU .GT. 0.50 .AND. SANU .LE. 1.0)          THEN
                X                 =  SQRT(ANUTAU)
                DTRGW             =  X + 2.0 * (SANU - 0.50) *
     1                              (ANUTAU - X)
              ELSE IF (SANU .GT. 1.0 .AND. SANU .LE. 2.0)           THEN
                DTRGW             =  ANUTAU + (SANU - 1.0) * ANUTAU *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 2.0 .AND. SANU .LE. 3.0)           THEN
                X                 =  ANUTAU * ANUTAU
                DTRGW             =  X + (SANU - 2.0) * X *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 3.0 .AND. SANU .LE. 4.0)           THEN
                X                 =  ANUTAU * ANUTAU * ANUTAU
                DTRGW             =  X + (SANU - 3.0) * X *
     1                              (ANUTAU - 1.0)
              ELSE IF (SANU .GT. 4.0 .AND. SANU .LE. 20.0)          THEN
                DTRGW             =  ANUTAU ** (NINT(SANU))
              ELSE
                DTRGW             =  EXP(-(TAUM(I,2,KM1)-TAUM(I,2,KK)))
              ENDIF
C
              S1(I)               =  S1(I) - 
     1                              (EMISW(I,KK) - S(I,KK)) * DTRGW
            ENDIF
  720     CONTINUE
        ENDIF
C
        DO 730 I = IL1, IL2
          IF (KM2 .GE. NCT(I))                                      THEN
            IF (CLD(I,KM1) .GE. CUT)                                THEN
              IF (NCD(I,KM1) .EQ. 1)                                THEN
                FW(I,1,K)         =  FX(I,2,K)
                FD(I,2,K)         =  FX(I,1,K) + CLD(I,KM1) * 
     1                              (FX(I,2,K) - FX(I,1,K)) 
              ELSE
                FW(I,1,K)         =  TERM1(I) + S1(I) + BF(I,K) - 
     1                               EMISW(I,KM1) + S(I,KM1)
                FD(I,2,K)         =  CLDM(I,KM1) * 
     1                              (FW(I,1,K) - FX(I,2,K)) + 
     2                               FX(I,1,K) + CLD(I,KM1) * 
     3                              (FX(I,2,K) - FX(I,1,K))
              ENDIF
            ENDIF
          ENDIF
  730   CONTINUE
  750 CONTINUE
C
      RETURN
      END
