      SUBROUTINE LWTRAN_MCICA2(FUT, FDT,
     1                         SLWF, TAUC, OMC, GC, FL, TAUAL,
     2                         TAUG, BF, BST, URBF, DBF, EM0T,
     3                         CLD, NCT, LEV1, CUT, MAXC,
     4                         ITILE,
     5                         IL1, IL2, ILG, LAY, LEV,
     6                         NTILE)
C
C     * JUN 02,2015 - M.LAZARE/ NEW VERSION FOR GCM19:
C     *               J.COLE:   - ADD TILED RADIATION CALCULATIONS
C     *                           (IE "FUT","FDT"), UNDER CONTROL OF
C     *                           "ITILRAD".
C     * FEB 11,2009 - J.COLE.  PREVIOUS VERSION LWTRAN_MCICA2 FOR GCM15H
C     *                        THROUGH GCM18:
C     * FEB 20/2007 - J. COLE. PREVIOUS VERSION LWTRAN_MCICA FOR GCM15G:
C     *                        - UPDATE TO MATCH LWTRAN2, AUTOMATIC ARRAYS.
C     *                        - USED PETRI'S HARD WORK AND MODIFIED TO 
C     *                          SUIT MY VERSION. NO LONGER HAS SUBCOLUMNS,
C     *                          WORKS ON ONLY ONE SUBCOLUMN PER GCM 
C     *                          COLUMN.
C     * MAY 31/2004 - P.RAISANEN:
C 
C     THIS CODE VERSION HAS BEEN UPDATED FOR McICA RADIATION CALCULATIONS
C
C       - Clear-sky fluxes are calculated using mean profiles for the GCM 
C         column, but cloudy-sky fluxes are calculated separately
C         for NSUB subcolumns and then averaged. It is assumed that
C         there are equally many subcolumns for each GCM column IL1...IL2.
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
C     * FEB 04,2004 - J.LI, M.LAZARE.
C----------------------------------------------------------------------C
C     CALCULATION OF LONGWAVE RADIATIVE TRANSFER USING ABSORPTION      C
C     APPROXIMATION, BUT INCLUDING CORRECTIONS FOR SCATTERING          C   
C     (UNPERTURBED + BACKWARD SCATTERING EFFECT +  FORWARD SCATTERING  C
C      EFFECT + INTERNAL SCATTERING EFFECT  (LI AND FU, JAS 2000)      C
C----------------------------------------------------------------------C
C     FU:      UPWARD INFRARED FLUX                                    C
C     FD:      DOWNWARD INFRARED FLUX                                  C
C    FUT:      TILED UPWARD   INFRARED FLUX                            C
C    FDT:      TILED DOWNWARD INFRARED FLUX                            C
C     SLWF:    INPUT SOLAR FLUX AT MODEL TOP LEVEL FOR EACH BAND       C
C     TAUCSUB: CLOUD SUBCOLUMN OPTICAL DEPTH FOR THE INFRARED          C
C     OMCSUB:  CLOUD SUBCOLUMN SINGLE SCATTERING ALBEDO TIMES OPTICAL  C
C              DEPTH                                                   C
C     GCSUB:   CLOUD SUBCOLUMN ASYMMETRY FACTOR TIMES OMCI             C
C     FLSUB:   SQUARE OF CLOUD SUBCOLUMN 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 NEEDEDC
C              FOR FLUX                                                C
C     BST:     THE BLACKBODY INTENSITY FOR EACH TILED 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     EM0T:    SURFACE EMISSION FOR EACH TILED SURFACE                 C
C     CLDSUB:  SUBCOLUMN CLOUD FRACTION (ASSUMED TO BE 0 OR 1)         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     NCTSUB:  THE HIGHEST CLOUD TOP LEVEL FOR SUBCOLUMNS              C
C     MAXC:    MINIMUM VALUE OF NCT FOR ALL COLUMNS AND SUBCOLUMNS     C      
C     NSUB:    NUMBER OF SUBCOLUMNS PER GCM COLUMN                     C
C     NSUBMX:  MAXIMUM NUMBER OF SUBCOLUMNS PER GCM COLUMN             C  
C----------------------------------------------------------------------C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
      REAL FUT(ILG,NTILE,2,LEV), FDT(ILG,NTILE,2,LEV)
      REAL SLWF(ILG), TAUC(ILG,LAY), OMC(ILG,LAY), 
     1     GC(ILG,LAY), FL(ILG,LAY), 
     2     TAUAL(ILG,LAY), TAUG(ILG,LAY), BF(ILG,LEV), 
     3     URBF(ILG,LAY), DBF(ILG,LAY),
     3     CLD(ILG,LAY)
      REAL BST(ILG,NTILE), EM0T(ILG,NTILE)
      INTEGER ITILE(ILG,NTILE)
      REAL SCATSM(ILG,2,LAY), TAUM(ILG,2,LAY),
     2     XD(ILG,2,LAY), XU(ILG,2,LAY), DTR(ILG,2,LAY)
      INTEGER NCT(ILG)
C
C     * INTERNAL WORK ARRAYS.
C
      REAL EMISW(ILG,LAY), SCATBK(ILG,LAY), SCATFW(ILG,LAY)
      REAL EMBST(ILG,NTILE), ABSE0T(ILG,NTILE)
      REAL FD(ILG,2,LEV), FU(ILG,2,LEV)
C
      DATA  RU / 1.6487213 /
C=======================================================================
C
C-------------------
C INITIALIZATIONS
C-------------------
C
      L1 =  LEV1
      L2 =  LEV1 + 1
      
      DO 10 K = L1, LEV
         DO 20 I = IL1,IL2
            FU(I,2,K) = 0.
            FD(I,2,K) = 0.
 20      CONTINUE
 10   CONTINUE
      DO K = LEV1,LEV
         DO M = 1, NTILE
            DO I = IL1, IL2
               FDT(I,M,1,K) = 0.0
               FUT(I,M,1,K) = 0.0               
               FDT(I,M,2,K) = 0.0
               FUT(I,M,2,K) = 0.0
            END DO ! I
         END DO ! M
      END DO ! K

      DO M = 1, NTILE
         DO I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
               EMBST(I,M)  = EM0T(I,M) * BST(I,M)
               ABSE0T(I,M) = 1.0 - EM0T(I,M)
            END IF
         END DO ! I
      END DO ! M
C
C----------------------------------------------------------------------C
C     CALCULATE THE DOWNWARD FLUXES FIRST WITHOUT SCATTERING           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.                                   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
C----------------------------------------------------------------------C
C COMPUTATION OF CLEAR-SKY DOWNWARD FLUXES USING EXPONENTIAL SOURCE    C
C FUNCTION                                                             C
C----------------------------------------------------------------------C

      L1 =  LEV1
      L2 =  LEV1 + 1

      DO 100 I = IL1, IL2
         FD(I,1,LEV1) =  SLWF(I)
 100  CONTINUE
C
      DO 110 K = L2, LEV
         KM1 = K - 1
         DO 120 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)
 120     CONTINUE
 110  CONTINUE 
      DO K = LEV1,LEV
         DO M = 1, NTILE
            DO I = IL1, IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  FDT(I,M,1,K) = FD(I,1,K)
               END IF 
            END DO ! I
         END DO ! M
      END DO ! K

C----------------------------------------------------------------------C
C COMPUTATION OF CLEAR-SKY UPWARD FLUXES                               C
C----------------------------------------------------------------------C

      DO M = 1, NTILE
         DO I = IL1,IL2
            IF (ITILE(I,M) .GT. 0) THEN
               FUT(I,M,1,LEV) = EMBST(I,M) + ABSE0T(I,M)*FDT(I,M,1,LEV)
            END IF
         END DO ! I
      END DO ! M

      DO K = LEV-1, L1, -1
         KP1 = K+1
         DO M = 1, NTILE
            DO I = IL1,IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  FUT(I,M,1,K)     = FUT(I,M,1,KP1) * DTR(I,1,K) +
     1                                XU(I,1,K)
               END IF
            END DO ! I
         END DO ! M
      END DO ! K

C----------------------------------------------------------------------C
C ALL-SKY DOWNWARD FLUXES ABOVE THE HIGHEST CLOUD TOP                  C
C----------------------------------------------------------------------C

      DO 400 I = IL1,IL2
         FD(I,2,L1)               = SLWF(I) 
 400  CONTINUE

      DO 410 K = L2, LEV 
         DO 420 I = IL1,IL2
            IF (K.LE.NCT(I)) THEN
               FD(I,2,K)            = FD(I,1,K)
            END IF
 420     CONTINUE
 410  CONTINUE  

C----------------------------------------------------------------------C
C     ADD THE LAYERS DOWNWARD FROM THE HIGHEST CLOUD LAYER TO THE      C
C     SURFACE. DETERMINE THE XU FOR THE UPWARD PATH.                   C
C     USING EXPONENTIAL SOURCE FUNCTION ALSO FOR ALL SKY FLUX IN CLOUD C
C     FREE LAYERS.                                                     C
C----------------------------------------------------------------------C
C
      DO 430 K = MAXC+1, LEV
         KM1 = K - 1
         DO 440 I = IL1, IL2
            IF (K.GT.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)
               ELSE
                  TAUL2               =  TAUC(I,KM1) + TAUG(I,KM1)
                  SSALB               =  OMC(I,KM1) / TAUL2
                  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 
                  DTR(I,2,KM1)        =  EXP(-TAUM(I,1,KM1))
                  DTR2                =  EXP(-TAU2)
                  EMISW(I,KM1)        =  ZETA * (1.0 - DTR(I,2,KM1))
                  EMBK                =  DTR2 - 1.0 
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 * GC(I,KM1) / COW
                  TAUDTR              =  TAUM(I,1,KM1) * DTR(I,2,KM1)
                  SCATFW(I,KM1)       =  WGRCOW * 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
                  FD(I,2,K)           =  FD(I,2,KM1) * DTR(I,2,KM1) + 
     1                                  XD(I,2,KM1)
               ENDIF
            ENDIF
 440     CONTINUE
 430  CONTINUE
      DO K = L1,LEV
         DO M = 1, NTILE
            DO I = IL1, IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  FDT(I,M,2,K) = FD(I,2,K)
               END IF 
            END DO ! I
         END DO ! M
      END DO ! K

C----------------------------------------------------------------------C
C COMPUTATION OF ALL-SKY UPWARD FLUXES                                 C
C----------------------------------------------------------------------C
C INITIALIZATION FOR SURFACE                                           C
C----------------------------------------------------------------------C
      K = LEV - 1
      DO M = 1, NTILE
         DO I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
               FUT(I,M,2,LEV) = EMBST(I,M) + ABSE0T(I,M)*FDT(I,M,2,LEV)
            END IF
         END DO ! I
      END DO ! M

C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD FROM THE FIRST LAYER TO MAXC               C
C     SCATTERING EFFECT FOR UPWARD PATH IS INCLUDED                    C
C----------------------------------------------------------------------C
C
      DO K = LEV - 1, MAXC, - 1
         KP1 = K + 1
         DO M = 1, NTILE
            DO I = IL1, IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  IF (K .GE. NCT(I)) THEN
                     IF (CLD(I,K) .LT. CUT) THEN
                        FUT(I,M,2,K) =  FUT(I,M,2,KP1) * DTR(I,1,K) +
     1                                  XU(I,1,K)
                     ELSE
                        FUT(I,M,2,K) =  FUT(I,M,2,KP1) *
     1                                 (DTR(I,2,K) + SCATFW(I,K)) +
     2                                  FDT(I,M,2,K) * SCATBK(I,K) +
     3                                  SCATSM(I,2,K) + XU(I,2,K)
                     ENDIF
                  END IF
               END IF
            END DO ! I
         END DO ! M
      END DO ! K
C
C----------------------------------------------------------------------C
C     ADD THE LAYERS UPWARD ABOVE THE HIGHEST CLOUD TO THE TOA, NO     C
C     SCATTERING                                                       C
C----------------------------------------------------------------------C
C
      DO K = LEV - 1, L1, - 1
         KP1 = K + 1
         DO M = 1, NTILE
            DO I = IL1, IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  IF (KP1 .LE. NCT(I)) THEN
                     FUT(I,M,2,K) = FUT(I,M,2,KP1) * DTR(I,1,K) +
     1                              XU(I,1,K)
                  ENDIF
               END IF
            END DO ! I
         END DO ! M 
      END DO ! K
C
C----------------------------------------------------------------------C
C     SCATTERING EFFECT FOR DOWNWARD PATH IN FROM MAXC TO THE SURFACE  C
C----------------------------------------------------------------------C
C
      DO K = MAXC + 1, LEV
         KM1 = K - 1
         DO M = 1, NTILE
            DO I = IL1, IL2
               IF (ITILE(I,M) .GT. 0) THEN
                  IF (KM1 .GE. NCT(I)) THEN
                     IF (CLD(I,KM1) .LT. CUT) THEN
                        FDT(I,M,2,K) = FDT(I,M,2,KM1) * DTR(I,1,KM1) +
     1                                 XD(I,1,KM1)
                     ELSE
                        FDT(I,M,2,K) = FDT(I,M,2,KM1) *
     1                                (DTR(I,2,KM1) + SCATFW(I,KM1)) +
     2                                 FUT(I,M,2,K) * SCATBK(I,KM1) +
     3                                 SCATSM(I,1,KM1) + XD(I,2,KM1)
                     ENDIF
                  ENDIF  
               END IF
            END DO ! I
         END DO ! M
      END DO ! K
C

      RETURN
      END
