      SUBROUTINE LWTRAGH4(FUT, FDT, 
     1                    SLWF, TAUCI, OMCI, TAUAL, TAUG, BF,
     2                    URBF, CLD, EM0T, BST, ITILE, CUT,
     3                    IL1, IL2, ILG, LAY, LEV,
     4                    NTILE)
C
C     * JUN 02,2015 - M.LAZARE/ NEW VERSION FOR GCM19:
C     *               J.COLE:   - ADD TILED RADIATION CALCULATIONS
C     *                           (IE "FUT","FDT")
C     * FEB 11,2009 - J.COLE.  PREVIOUS VERSION LWTRAGH4 FOR GCM15H
C     *                        THROUGH GCM18:
C     *                         - CORRECT BUG IN SPECIFICATION OF ABSE0.
C     *                         - INITIALIZE CLEAR-SKY FX IN TOP LAYER.
C     *                         - CHANGE ANY WORK ARRAYS HAVING A MIDDLE 
C     *                           DIMENSION OF "4" TO "2", AS THIS IS
C     *                           WHAT IS USED.   
C     * DEC 05,2007 - M.LAZARE. PREVIOUS VERSION LWTRAGH3 FOR GCM15G:
C     *                         - SUPPORT ADDED FOR EMISSIVITY<1.  
C     * NOV 22/2006 - M.LAZARE. PREVIOUS VERSION LWTRAGH2 FOR GCM15F:
C     *                         - BOUND OF 1.E-10 USED FOR 1-CLD
C     *                           INSTEAD OF SUBTRACTING AN EPSILON.
C     *                         - WORK ARAYS FOR THIS ROUTINE NOW
C     *                           ARE AUTOMATIC ARRAYS AND THEIR
C     *                           SPACE IS NOT PASSED IN.   
C     * AUG 29,2003 - J.LI.       PREVIOUS VERSION LWTRAGH UP TO GCM15E.
C----------------------------------------------------------------------C
C     IN THE G SPACE WITH INTERVAL CLOSE 1 (VERY LARGE OPTICAL DEPTH)  C
C     OR IN THE CASE WITH CLOUD ABSORPTION IS VERY SMALL OR THE WEIGHT C
C     OF FLUX AND COOLING RATE ARE VERY SMALL. THE CLOUD RADIATIVE     C
C     PROCESS CAN BE HIGHLY SIMPLIFIED. THE ABSORPTION APPROXIMATION   C
C     METHOD IS USED AND CLOUD RANDOM AND MAXIMUM OVERLAP IS           C
C     CONSIDERED, BUT CLOUD SCATTERING AND INHOMOGENEITY ARE IGNORED.  C
C     THE EXPONENTIAL SOURCE PLANCK FUNCTION IS USED WHICH IS MORE     C
C     ACCURATE IN THE REGION ABOVE 200 MB IN COMPARISON WITH LINEAR    C
C     SOURCE FUNCTION                                                  C
C                                                                      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     TAUCI: CLOUD OPTICAL DEPTH FOR THE INFRARED                      C
C     OMCI:  CLOUD SINGLE SCATTERING ALBEDO TIMES OPTICAL DEPTH        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.                              C
C     BS:    THE BLACKBODY INTENSITY AT THE SURFACE.                   C
C     BST:    THE BLACKBODY INTENSITY FOR EACH TILED SURFACE.          C
C            USED FOR EXPONENTIAL SOURCE FUNCTION (LI, 2002 JAS P3302) C
C     CLD:   CLOUD FRACTION                                            C
C     EM0T:   SURFACE EMISSION FOR EACH TILED SURFACE                  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     FYT:   UPWARD FLUX FOR PURE CLEAR PORTION (1) AND PURE CLOUD     C
C            PORTION (2) FOR EACH TILE                                 C
C     FX:    THE SAME AS FY BUT FOR THE DOWNWARD FLUX                  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),
     2       TAUAL(ILG,LAY), TAUG(ILG,LAY), BF(ILG,LEV), URBF(ILG,LAY),
     3       CLD(ILG,LAY)
      REAL   XU(ILG,2,LAY),XD(ILG,2,LAY),DTR(ILG,2,LAY),
     1       FX(ILG,2,LEV)
      REAL FUT(ILG,NTILE,2,LEV),FDT(ILG,NTILE,2,LEV),
     1     FYT(ILG,NTILE,2,LEV)
      REAL EM0T(ILG,NTILE), BST(ILG,NTILE)
      INTEGER ITILE(ILG,NTILE)
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                                    C
C     THE OVERLAP BETWEEN SOLAR AND INFRARED IN 4 - 10 UM IS           C
C     CONSIDERED, SLWF IS THE INCOMING SOLAR FLUX                      C
C     SINGULARITY FOR XD AND XU HAS BEEN CONSIDERED AS LI JAS 2002     C
C----------------------------------------------------------------------C
C
      DO 100 I = IL1, IL2
        FD(I,1,1)         =  SLWF(I)
        FD(I,2,1)         =  SLWF(I)
        FX(I,1,1)         =  SLWF(I)
        FX(I,2,1)         =  SLWF(I)
C
        TAUL1             =  TAUAL(I,1) + TAUG(I,1)
        RTAUL1            =  TAUL1 * RU
        DTR(I,1,1)        =  EXP (- RTAUL1)
        UBETA             =  URBF(I,1) / (TAUL1 + 1.E-20)
        EPSD              =  UBETA + 1.0
        EPSU              =  UBETA - 1.0
C
        IF (ABS(EPSD) .GT. 0.001)                                   THEN
          XD(I,1,1)       = (BF(I,2) - BF(I,1) * DTR(I,1,1)) / EPSD
        ELSE
          XD(I,1,1)       =  RTAUL1 * BF(I,1) * DTR(I,1,1)
        ENDIF
        IF (ABS(EPSU) .GT. 0.001)                                   THEN
          XU(I,1,1)       = (BF(I,2) * DTR(I,1,1) - BF(I,1)) / EPSU
        ELSE
          XU(I,1,1)       =  RTAUL1 * BF(I,2) * DTR(I,1,1)
        ENDIF
C
        FD(I,1,2)         =  FD(I,1,1) * DTR(I,1,1) + XD(I,1,1)
C
        IF (CLD(I,1) .LT. CUT)                                      THEN
          FX(I,1,2)       =  FD(I,1,2)
          FX(I,2,2)       =  FD(I,1,2)
          FD(I,2,2)       =  FD(I,1,2)
        ELSE
          TAUL2           =  TAUCI(I,1) + TAUL1
          COW             =  1.0 - OMCI(I,1) / TAUL2
          CTAUL2          =  COW * TAUL2
          CRTAUL2         =  CTAUL2 * RU
          DTR(I,2,1)      =  EXP (- CRTAUL2)
          UBETA           =  URBF(I,1) / (CTAUL2)
          EPSD            =  UBETA + 1.0
          EPSU            =  UBETA - 1.0
C
          IF (ABS(EPSD) .GT. 0.001)                                 THEN
            XD(I,2,1)     = (BF(I,2) - BF(I,1) * DTR(I,2,1)) / EPSD
          ELSE
            XD(I,2,1)     =  CRTAUL2 * BF(I,1) * DTR(I,2,1)
          ENDIF
          IF (ABS(EPSU) .GT. 0.001)                                 THEN
            XU(I,2,1)     = (BF(I,2) * DTR(I,2,1) - BF(I,1)) / EPSU
          ELSE
            XU(I,2,1)     =  CRTAUL2 * BF(I,2) * DTR(I,2,1)
          ENDIF
C
          FX(I,1,2)       =  FX(I,1,1) * DTR(I,1,1) + XD(I,1,1)
          FX(I,2,2)       =  FX(I,2,1) * DTR(I,2,1) + XD(I,2,1)
          FD(I,2,2)       =  FX(I,1,2) +
     1                       CLD(I,1) * (FX(I,2,2) - FX(I,1,2))
        ENDIF
  100 CONTINUE
C
      DO 250 K = 3, LEV
        KM1 = K - 1
        KM2 = KM1 - 1
        DO 200 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) * DTR(I,1,KM1)) / EPSD
          ELSE
            XD(I,1,KM1)   =  RTAUL1 * BF(I,KM1) * DTR(I,1,KM1) 
          ENDIF
          IF (ABS(EPSU) .GT. 0.001)                                 THEN
            XU(I,1,KM1)   = (BF(I,K) * DTR(I,1,KM1) - 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) + XD(I,1,KM1)
C
          IF (CLD(I,KM1) .LT. CUT)                                  THEN
            FD(I,2,K)     =  FD(I,2,KM1) * DTR(I,1,KM1) + 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) + TAUL1
            COW           =  1.0 - OMCI(I,KM1) / TAUL2
            CTAUL2        =  COW * TAUL2
            CRTAUL2       =  CTAUL2 * RU
            DTR(I,2,KM1)  =  EXP (- CRTAUL2)
            UBETA         =  URBF(I,KM1) / (CTAUL2)
            EPSD          =  UBETA + 1.0
            EPSU          =  UBETA - 1.0
C
            IF (ABS(EPSD) .GT. 0.001)                               THEN
              XD(I,2,KM1) = (BF(I,K) - BF(I,KM1) * DTR(I,2,KM1)) / EPSD
            ELSE
              XD(I,2,KM1) =  CRTAUL2 * BF(I,KM1) * DTR(I,2,KM1)
            ENDIF
            IF (ABS(EPSU) .GT. 0.001)                               THEN
              XU(I,2,KM1) = (BF(I,K) * DTR(I,2,KM1) - BF(I,KM1)) / EPSU
            ELSE
              XU(I,2,KM1) =  CRTAUL2 * BF(I,K) * DTR(I,2,KM1)
            ENDIF
C
            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) + 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) + XD(I,1,KM1)
              FX(I,2,K)   = (FX(I,1,KM1) + CLD(I,KM2) / CLD(I,KM1) *
     1                      (FX(I,2,KM1) - FX(I,1,KM1))) *
     2                       DTR(I,2,KM1) + XD(I,2,KM1)
            ENDIF
C
            FD(I,2,K)     =  FX(I,1,K) + CLD(I,KM1) * (FX(I,2,K) -
     1                       FX(I,1,K))
          ENDIF
  200   CONTINUE
  250 CONTINUE
C
      DO K = 1, 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)
                  FDT(I,M,2,K) = FD(I,2,K)
               END IF
            END DO ! I
         END DO ! M
      END DO ! K

      DO M = 1, NTILE
         DO I = IL1, IL2
            IF (ITILE(I,M) .GT. 0) THEN
               EMBS             =  EM0T(I,M) * BST(I,M)
               ABSE0            =  1.0 - EM0T(I,M)
               FUT(I,M,1,LEV)   =  EMBS + ABSE0 * FDT(I,M,1,LEV)
               FYT(I,M,1,LEV)   =  EMBS + ABSE0 * FX(I,1,LEV)
               FYT(I,M,2,LEV)   =  EMBS + ABSE0 * FX(I,2,LEV)
               FUT(I,M,2,LEV)   =  FYT(I,M,1,LEV) +
     1                             CLD(I,LAY) * 
     2                             (FYT(I,M,2,LEV) - FYT(I,M,1,LEV))
C
               FUT(I,M,1,LAY)   =  FUT(I,M,1,LEV) * DTR(I,1,LAY) + 
     1                             XU(I,1,LAY)
               FYT(I,M,1,LAY)   =  FYT(I,M,1,LEV) * DTR(I,1,LAY) + 
     1                             XU(I,1,LAY)
C
               IF (CLD(I,LAY) .LT. CUT)                            THEN
                  FYT(I,M,2,LAY) =  FYT(I,M,2,LEV) * DTR(I,1,LAY) + 
     1                              XU(I,1,LAY)
                  FUT(I,M,2,LAY) =  FYT(I,M,1,LAY)
               ELSE
                  FYT(I,M,2,LAY) =  FYT(I,M,2,LEV) * DTR(I,2,LAY) + 
     1                              XU(I,2,LAY)
                  FUT(I,M,2,LAY) =  FYT(I,M,1,LAY) +
     1                              CLD(I,LAY) * 
     2                              (FYT(I,M,2,LAY) - FYT(I,M,1,LAY))
               ENDIF
            ENDIF
         END DO ! I
      END DO ! M
C
      DO K = LEV - 2, 1, - 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) + XU(I,1,K)
C
                 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)
                    FYT(I,M,1,K) = FUT(I,M,2,K)
                    FYT(I,M,2,K) = FUT(I,M,2,K)
                 ELSE
                    IF (CLD(I,K) .LT. CLD(I,KP1))                   THEN
                     FYT(I,M,1,K) = (FYT(I,M,2,KP1) + (1.0 - CLD(I,KP1))
     1                            / (1.0 - CLD(I,K)) * (FYT(I,M,1,KP1)
     2                            -  FYT(I,M,2,KP1)) ) * DTR(I,1,K)
     3                            +  XU(I,1,K)
                     FYT(I,M,2,K) = FYT(I,M,2,KP1) * DTR(I,2,K) + 
     1                              XU(I,2,K)
                    ELSE
                     FYT(I,M,1,K) = FYT(I,M,1,KP1) * DTR(I,1,K) + 
     1                              XU(I,1,K)
                     FYT(I,M,2,K) = (FYT(I,M,1,KP1) + CLD(I,KP1) 
     1                            / CLD(I,K)
     2                            * (FYT(I,M,2,KP1) - FYT(I,M,1,KP1)))
     3                            * DTR(I,2,K) + XU(I,2,K)
                   ENDIF
C
                  FUT(I,M,2,K) = FYT(I,M,1,K) +
     1                           CLD(I,K)*(FYT(I,M,2,K) - FYT(I,M,1,K))
               ENDIF
              ENDIF
            END DO ! I
         END DO ! M
      END DO ! K
C
      RETURN
      END
