      SUBROUTINE RADDRIV10(HRS, HRL, HRSC, HRLC,
     1                    HRSH2O, HRSO3, HRSCO2, HRSO2, HRSCH4,
     2                    HRLH2O, HRLO3, HRLCO2, HRLCH4, HRLN2O,
     +                    FSAU,FSAD,FLAU,FLAD,
     +                    FSCU,FSCD,FLCU,FLCD,
     4                    FSG, FSD, FSF, FSV, FSI, FDL,
     5                    FLG, FDLC, CSB, CLB, FSR, FSRC, OLR, OLRC,
     6                    PAR, CSD, CSF, FSLO, FSAMOON, FLAMOON, FSO,
     7                    FSDB, FSFB, CSDB, CSFB, FSSB, FSSCB,
     8                    WRKA, WRKB,
     +                    FSGT, FSDT, FSFT, FSVT, FSIT,
     +                    FDLT, FLGT, FDLCT, CSBT, CLBT,
     +                    PART, CSDT, CSFT,
     +                    FSDBT, FSFBT, CSDBT, CSFBT, FSSBT, FSSCBT,
     +                    TBND, EXTA055,
     9                    SHTJ, TFULL, SHJ, DSHJ, DZ, T, O3, Q,
     A                    CO2, CH4, AN2O, F11, F12, F113, F114,
     B                    CLD, RH, ANU, ETA, AERIN, IDIMR, IPAM,
     C                    EXB1,EXB2,EXB3,EXB4,EXB5,EXBT,
     D                    ODB1,ODB2,ODB3,ODB4,ODB5,ODBT,ODBV,
     E                    OFB1,OFB2,OFB3,OFB4,OFB5,OFBT,
     F                    ABB1,ABB2,ABB3,ABB4,ABB5,ABBT,
     G                    EXS1,EXS2,EXS3,EXS4,EXS5,EXST,
     H                    ODS1,ODS2,ODS3,ODS4,ODS5,ODST,ODSV,
     I                    OFS1,OFS2,OFS3,OFS4,OFS5,OFST,
     J                    ABS1,ABS2,ABS3,ABS4,ABS5,ABST,
     K                    ISVAOD,IMAM,
     L                    SALB, CSAL, PRESSG, GT, O3TOP, RMU, SOLV,
     M                    IEXPLVOL, IVTAU, VTAU, TROP,
     +                    SW_EXT_SA, SW_SSA_SA, SW_G_SA, LW_ABS_SA,
     N                    EM0, CLDT, RADJ, WCDW,
     O                    REL_SUB, REI_SUB, CLW_SUB, CIC_SUB, NCLDY,
     +                    SALBT, CSALT, EM0T, GTT, FARET,
     Q                    COLNUM_SW, COLNUM_LW, MAX_SAM, ISEEDROW,
     R                    LCSW, LCLW, IL1, IL2, ILG, LAY, LEV, NXLOC,
     S                    NTILE, SOLAR_C, CLDWATMIN, IVERS, JLAT, KOUNT,
     T                    MCICA, IRADFORCE, IRAD_FLUX_PROFS, IACTIVE_RT,
     +                    ITILERAD)
C
C     * Nov 02/2018 - J. Cole     - Add radiative flux profiles
C     * Nov 01/2017 - J. Cole     - Added CMIP6 stratospheric aerosols.
C     * Feb 10/2015 - M.Lazare/   New version for gcm18:
C     *               J.Cole/     - KSALB and ZSPD not needed, and
C     *               J.LI.         also remove PBLT since never used.
C     *                           - WRKA,WRKB now dimensioned over IB
C     *                             and passed out to physics, etc,
C     *                             ie no longer a work field.
C     *                           - Calls to new subroutines
C     *                             SW_CLD_PROPS4 and LW_CLD_PROPS3.
C     *                           - Compute the direct beam using the
C     *                             unscaled cloud optical properties.
C     * Jun 25,2013 - M.Lazare/   New version for gcm17:
C     *               J.Cole/     - Compute and save the band-mean
C     *               J.Li/         values for FSD,FSF,FSSB,CSD,CSF.
C     *               K.Vonsalzen.- Support for PAM.
C     *                           - Calculate heating rates in usual
C     *                             manner (aggregate into HRS or HRL
C     *                             rather than separate gas band
C     *                             interactions) if not running MAM.
C     *                           - Include semi-direct effect.
C     *                           - Calls to revised new routines:
C     *                             {BCAERO4,OCAERO4,SW_CLD_PROPS3,
C     *                             SATTENU6,LATTENU6}.
C     * MAY 01,2012 - J.LI/       PREVIOUS VERSION RADDRIV8 FOR GCM16:
C     *               J.COLE/     - REVISED CALLS FOR CHANGES TO
C     *               M.LAZARE.     ISCCP SIMULATOR.
C     *                           - REMOVE BACKGROUND AEROSOLS.
C     *                           - ADD DIAGNOSTIC CAPABILITY FOR
C     *                             AEROSOL OPTICAL DEPTH.
C     *                           - MODIFICATIONS TO RADIATIVE FORCING.
C     *                           - ADD SOLAR CONTINUUM, WHICH ALSO
C     *                             REQUIRES CHANGES TO NEW VERSIONS
C     *                             GASOPTS5 AND CKDSW4 CALLS.
C     *                           - REVISED SOLAR VARIABLITY TAKEN
C     *                             FROM LEAN'S RECOMMENDATION.
C     *                           - INTEGRATION OF CORRELATED-K WITH
C     *                             MAM RADIATION.
C     *                           - BUGFIX TO INCLUDE SECOND SUB-BAND
C     *                             OF IB=1 IN PAR.
C     *                           - ISOTHERMAL, NOT EXTRAPOLATE, FOR
C     *                             MOON LAYER (IN THIS ROUTINE AND
C     *                             ALSO NEW PLANCK3).
C     *                           - IMPROVEMENTS TO RADIATIVE FLUXES,
C     *                             WHICH INVOLVE SAVING
C     *                             {FSR,FSRC,OLR,OLRC} AND NOT
C     *                             CALCULATING {FSA,FLA} FROM VERTICALLY
C     *                             INTEGRATING HEATING RATES. AS WELL,
C     *                             ABSORPTION IN MOON LAYER OF L/W
C     *                             RADIATION IS ALSO ADDED.
C     * APR 21,2010 - J.LI/       PREVIOUS VERSION RADDRIV7 FOR GCM15I:
C     *               J.COLE/     - ADD CALCULATION OF "TRUE" DIAGNOSTIC
C     *               M.LAZARE.     FSO INCLUDING SOLAR VARIABILITY.
C     *                             THIS REQUIRES AN ADDITIONAL DATA
C     *                             STATEMENT IN NEW CALLED CKDSW3 FOR
C     *                             INCIDENT SOLAR IN EACH BAND (NEW
C     *                             COMMON BLOCK "SWBAND").
C     *                           - CHANGES TO THE S/W OPTICAL PROPERTIES
C     *                             OF EXPLOSIVE VOLCANOES BY CHANGING
C     *                             PARTICLE SIZES FROM 0.166 MICRONS TO
C     *                             0.35 MICRONS.
C     *                           - ADD L/W EFFECT OF VOLCANOES.
C     *                           - IMPROVEMENTS FOR ACCURACY TO L/W
C     *                             (CHANGES IN CALLS TO NEW ROUTINES:
C     *                             CKDLW4,GASOPTL7,GASOPTLGH6).
C     * FEB 18,2009 - L.SOLHEIM/  PREVIOUS VERSION RADDRIV6 FOR GCM15H:
C     * FEB 18,2009 - M.LAZARE /  - IMPLEMENTATION OF 3D GHG (NEW CALLS).
C     *               J.LI.       - INCLUDE EFFECT OF EXPLOSIVE
C     *                             VOLCANOES UNDER CONTROL OF "IEXPLVOL".
C     *                           - INCLUDE POSSIBILITY OF SOLAR
C     *                             VARIABILITY.
C     *                           - NEW CALLS TO REVISED RAD ROUTINES.
C     *                           - RADCON_T COMMON BLOCK REMOVED AND
C     *                             SOLAR CONSTANT PASSED INSTEAD.
C     * APR 21,2008 - L.SOLHEIM/  - COSMETIC CHANGE TO ADD THREADPRIVATE
C     *               M.LAZARE.     FOR COMMON BLOCK "TRACE", IN SUPPORT
C     *                             OF "RADFORCE" MODEL OPTION.
C     *                           - ADD CALCULATION OF NEW OPTIONAL
C     *                             DIAGNOSTIC RADIATIVE FORCING FIELDS
C     *                             UNDER CONTROL OF NEW PASSED-IN
C     *                             SWITCH "IRADFORCE" (COSMETIC).
C     *                           - USE {IL1G,IL2G} INSTEAD OF {1,LENGATH}
C     *                             IN S/W SECTION, TO AVOID PASSING
C     *                             "1" TO SUBROUTINE(S) CALL (BAD FORM).
C     *                             ALSO COSMETIC.
C     *                           - COMBINED MCICA/NON-MCICA VERSION:
C     *                             WATER CONTENT AND EQUIVILENT
C     *                             RADIUS PASSED IN AND ALL
C     *                             OPTICAL PROPERTIES CALCULATED
C     *                             INSIDE IN NEW ROUTINES.
C     *                           - SURFACE EMISSIVITY FIELD (EM0) NOW
C     *                             PASSED IN AND USED (NON-UNITY ONLY
C     *                             OVER WATER) AND WORK FIELD REMOVED.
C     *                           - GUSTINESS FIELD (GUST) PASSED IN
C     *                             AND USED FOR SURFACE ALBEDO
C     *                             FORMULATION (IFF KSALB=1).
C     *                           - INTERNAL WORK FIELDS REMOVED IN
C     *                             CALLS TO ROUTINES.
C     *                           - METHANE EFFECT ADDED TO SOLAR.
C     *                           - MORE ACCURATE SOLAR CALCULATIONS
C     *                             (CONSTANTS, STRANDN2, STRANUP2).
C     *                           - EFFECT OF BLACK CARBON AND ORGANIC
C     *                             CARBON ADDED.
C     *                           - GUSTINESS FIELD NO LONGER PASSED
C     *                             IN SINCE EFFECT IS NOW ALREADY
C     *                             CONTAINED IN ZSPD FROM PHYSICS.
C     *                           - CALLS TO NEW ROUTINES.
C     * JAN 13/2007 - K.VONSALZEN/  PREVIOUS VERSION RADDRIV4 FOR GCM15F.
C     *               M.LAZARE/
C     *               J.LI/
C     *               P.VAILLANCOURT.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C
      PARAMETER (NBS = 4, NBL = 9, MAXNG = 10)
      REAL, PARAMETER :: M2KM = 0.001
      REAL, PARAMETER :: MIN_TEMP=120.0 ! Minimum temperature in Kelvins
      REAL, PARAMETER :: MAX_TEMP=400.0 ! Maximum temperature in Kelvins

C
C     * OUTPUT ARRAYS.
C
      REAL, DIMENSION(ILG,LAY) ::   HRS, HRL, HRSC, HRLC
      REAL, DIMENSION(ILG,LAY,3) :: HRSH2O
      REAL, DIMENSION(ILG,LAY,9) :: HRSO3, HRLH2O
      REAL, DIMENSION(ILG,LAY,2) :: HRSO2, HRLCO2
      REAL, DIMENSION(ILG,LAY)   :: HRSCO2, HRSCH4, HRLO3, HRLCH4,
     1                              HRLN2O, EXTA055

      REAL, DIMENSION(ILG)     ::   FSG, FSD, FSF, FSV, FSI,
     2                              FDL, FLG, FDLC, CSB, CLB, PAR,
     3                              FSR, FSRC, OLR, OLRC,
     4                              CSD, CSF, FSLO, FSAMOON, FLAMOON,
     5                              FSO
      REAL, DIMENSION(ILG,NBS)  ::   FSDB, FSFB, CSDB, CSFB, FSSB,
     +                               FSSCB, WRKA, WRKB
C
C     * TILED RADIATION ARRAYS.
C
C     * INPUT:
C
      REAL, DIMENSION(ILG,NTILE,NBS) :: SALBT, CSALT
      REAL, DIMENSION(ILG,NTILE) :: FARET, EM0T, GTT
C
C     * OUTPUT:
C
      REAL, DIMENSION(ILG,NTILE) :: FSGT, FSDT, FSFT, FSVT, FSIT,
     1                              FDLT, FLGT, FDLCT, CSBT, CLBT,
     2                              CSDT, CSFT, PART
      REAL, DIMENSION(ILG,NTILE,NBS) :: FSDBT, FSFBT, CSDBT,
     +                                  CSFBT, FSSBT, FSSCBT
C
C     * WORK ARRAYS.
C
      REAL  ,  DIMENSION(ILG,NTILE,LEV)     :: FLXUT,FLXDT
      REAL  ,  DIMENSION(ILG,NTILE,2,LEV)   :: REFLT,TRANT
      REAL  ,  DIMENSION(ILG,NTILE)         :: ALBSURT, CSALGT
      REAL  ,  DIMENSION(ILG,NTILE)         :: BST,EM0TL,GTTL, FAREL
      REAL  ,  DIMENSION(ILG,NTILE,NBS)     :: SALBTL, CSALTL
      INTEGER, DIMENSION(ILG,NTILE)         :: ITILE, ITILEG
C
C     * INPUT ARRAYS.
C
      REAL, DIMENSION(ILG,LEV)       :: SHTJ, TFULL
      REAL, DIMENSION(ILG,LAY)       :: SHJ, DSHJ, DZ, T, O3, Q,
     1                                  CO2, CH4, AN2O,
     2                                  F11, F12, F113, F114,
     3                                  CLD, RH, ANU, ETA, TBND
      REAL, DIMENSION(ILG,LAY,IDIMR) :: AERIN
      REAL, DIMENSION(ILG,NBS)       :: SALB, CSAL
      REAL, DIMENSION(ILG)           :: PRESSG, GT, O3TOP, RMU, EM0,
     1                                  CLDT, RADJ
      REAL, DIMENSION(ILG)           :: VTAU,TROP
      REAL, DIMENSION(ILG,LAY,NBS)   :: SW_EXT_SA, SW_SSA_SA, SW_G_SA
      REAL, DIMENSION(ILG,LAY,NBL)   :: LW_ABS_SA
      REAL, DIMENSION(12)            :: SOLV
      REAL, DIMENSION(ILG)     ::  EXB1,EXB2,EXB3,EXB4,EXB5,EXBT,
     1                             ODB1,ODB2,ODB3,ODB4,ODB5,ODBT,ODBV,
     2                             OFB1,OFB2,OFB3,OFB4,OFB5,OFBT,
     4                             ABB1,ABB2,ABB3,ABB4,ABB5,ABBT
C
      REAL, DIMENSION(ILG)     ::  EXS1,EXS2,EXS3,EXS4,EXS5,EXST,
     1                             ODS1,ODS2,ODS3,ODS4,ODS5,ODST,ODSV,
     2                             OFS1,OFS2,OFS3,OFS4,OFS5,OFST,
     4                             ABS1,ABS2,ABS3,ABS4,ABS5,ABST
C
C     * INPUT ARRAYS FOR MCICA.
C
      REAL, DIMENSION(ILG,LAY,NXLOC) :: REL_SUB, REI_SUB, CLW_SUB,
     1                                  CIC_SUB
      INTEGER, DIMENSION(ILG) :: NCLDY
C
C     * INPUT ARRAY FOR BC CALCULATIONS.
C
      REAL, DIMENSION(ILG,LAY)       :: WCDW
C
C     * RADIATIVE FLUX PROFILES UNDER CONTROL OF "IRAD_FLUX_PROFS"
C
      !--- Solar (FS*) and thermal (FL*) all-sky and clear-sky
      !--- upward and downward
      REAL, DIMENSION(ILG,LAY+2) :: FSAU,FSAD,FLAU,FLAD,
     1                              FSCU,FSCD,FLCU,FLCD

C INOUT arrays used for random number generator
      INTEGER, DIMENSION(ILG,4) :: ISEEDROW

! Arrays containing the indices of the subcolumns used for the
! SW and LW radiative transfer calculations

      REAL, DIMENSION(ILG,MAX_SAM) :: COLNUM_SW, COLNUM_LW
      INTEGER :: ICOLUMN
C
C****************************
C     * INTERNAL WORK ARRAYS:
C****************************
C
C     * GENERAL WORK ARRAYS.
C
      INTEGER, DIMENSION(NBS,MAXNG,2) :: NSAMPLE_SW
      INTEGER, DIMENSION(NBL,MAXNG,2) :: NSAMPLE_LW
      REAL  ,  DIMENSION(ILG,LAY)     :: PG,QG,P,DP,TAUR,TAUG,TAUA,
     1                                   F1,F2,TAUOMC,O3G,CO2G,CH4G,
     2                                   O2, DIP, DT, CLDG, O2G,
     3                                   EXTA086, FR1BC, FR2BC, FR1OC,
     4                                   FR2OC, RHOC, T_CLIP
      REAL  ,  DIMENSION(ILG,LEV)     :: FLXU,FLXD,PFULL,TF_CLIP
      REAL  ,  DIMENSION(ILG,2,LEV)   :: REFL,TRAN
      REAL  ,  DIMENSION(ILG,12)      :: A1
      REAL  ,  DIMENSION(ILG)         :: C1,C2,BS
      REAL  ,  DIMENSION(ILG,LAY,NBS) :: VSTAU,VSSSA,VSG
      REAL  ,  DIMENSION(ILG,LAY,NBL) :: VSABS
      REAL  ,  DIMENSION(ILG)         :: SUMSIG
      REAL  ,  DIMENSION(4)           :: SWEIGHT, VSSSCA, VSCGA
      REAL  ,  DIMENSION(9)           :: VOLABS
      INTEGER, DIMENSION(ILG)         :: MTOP
      REAL  ,  DIMENSION(12)          :: SVAR
      REAL  , DIMENSION(ILG,LAY) ::  AM3LOAD, REAM3, VEAM3, FR1, FR2,
     +                               SSLOAD, RESS, VESS, DSLOAD, REDS,
     +                               VEDS, BCLOAD, REBC, VEBC, OCLOAD,
     +                               REOC, VEOC, SO4LOAD, MSO4LOAD,
     +                               SSLOAD1, SSLOAD2, DSLOAD1,
     +                               DSLOAD2, BCOLOAD, BCYLOAD,
     +                               OCOLOAD, OCYLOAD, MOCOLOAD,
     +                               MOCYLOAD, RESO4, VESO4,
     +                               REBCY, VEBCY, REOCY, VEOCY,
     +                               REBCO, VEBCO, REOCO, VEOCO
C
      REAL  ,  DIMENSION(ILG,LAY,4):: TAUA055, TAUA086, ABSA055
C
C     * GATHERED AND OTHER WORK ARRAYS USED GENERALLY BY SOLAR.
C
      REAL  ,  DIMENSION(ILG,LAY,NBS) :: EXTA,EXOMA,EXOMGA,FA
      REAL  ,  DIMENSION(ILG,LAY,NBS) :: EXTB,ODFB,SSAB,ABSB
      REAL  ,  DIMENSION(ILG,LAY)     :: OMA055,GA055,
     1                                   EXTB055,ODF055,SSA055,ABS055
      REAL  ,  DIMENSION(ILG,LAY)     :: TAUCSG,TAUOMGC,CLDMG,
     1                                   TG,TAUOMA,TAUOMGA
      REAL  ,  DIMENSION(ILG,LEV)     :: PFULLG
      REAL  ,  DIMENSION(ILG,2,LEV)   :: CUMDTR
      REAL  ,  DIMENSION(ILG,12)      :: A1G
      REAL  ,  DIMENSION(ILG)         :: O3TOPG,ALBSUR,RMUG,CSALG,
     1                                   WRKAG,WRKBG
      INTEGER, DIMENSION(ILG,LAY)     :: INPTG,INPTMG,NBLK,NBLKG
      INTEGER, DIMENSION(ILG)         :: ISUN,NCTG
C
      REAL, DIMENSION(ILG,LAY)        :: BXT1ROL,BXT2ROL,BXT3ROL,
     1                                   BXT4ROL,BXT5ROL,BXTTROL
      REAL, DIMENSION(ILG,LAY)        :: BOD1ROL,BOD2ROL,BOD3ROL,
     1                                   BOD4ROL,BOD5ROL,BODTROL
      REAL, DIMENSION(ILG,LAY)        :: BOF1ROL,BOF2ROL,BOF3ROL,
     1                                   BOF4ROL,BOF5ROL,BOFTROL
      REAL, DIMENSION(ILG,LAY)        :: BAB1ROL,BAB2ROL,BAB3ROL,
     1                                   BAB4ROL,BAB5ROL,BABTROL
C
      REAL, DIMENSION(ILG,LAY)        :: SXT1ROL,SXT2ROL,SXT3ROL,
     1                                   SXT4ROL,SXT5ROL,SXTTROL
      REAL, DIMENSION(ILG,LAY)        :: SOD1ROL,SOD2ROL,SOD3ROL,
     1                                   SOD4ROL,SOD5ROL,SODTROL
      REAL, DIMENSION(ILG,LAY)        :: SOF1ROL,SOF2ROL,SOF3ROL,
     1                                   SOF4ROL,SOF5ROL,SOFTROL
      REAL, DIMENSION(ILG,LAY)        :: SAB1ROL,SAB2ROL,SAB3ROL,
     1                                   SAB4ROL,SAB5ROL,SABTROL
C
C     * WORK ARRAYS USED GENERALLY BY LONGWAVE.
C
      REAL  ,  DIMENSION(ILG,LAY,NBL) :: ABSA
      REAL  ,  DIMENSION(ILG,LAY)     :: TAUCI,OMCI,GCI,CLDM,URBF
      REAL  ,  DIMENSION(ILG,LEV)     :: BF
      INTEGER, DIMENSION(ILG,LAY)     :: INPT,INPTM,INPR,NCD,NCU
      REAL, DIMENSION(ILG,2)          :: TOA_OLR
      INTEGER, DIMENSION(ILG)         :: NCT
      INTEGER, DIMENSION(LAY)         :: NCUM,NCDM
C
C     * WORK ARRAYS USED BY MCICA
C
      REAL(8), DIMENSION(ILG)         :: RAN_NUM
      REAL, DIMENSION(ILG,LAY)        :: CF_MCICA,LWC,IWC,REL,REI
      INTEGER, DIMENSION(ILG)         :: RAN_INDEX, TMP_COUNT
C
C     * WORK ARRAYS USED FOR DIAGNOSIS OF AEROSOL OPTICAL PROPERTIES.
C
      REAL, DIMENSION (ILG,6):: AODS, EXTS, ABSS, MASS, ANGSTR, AOD86
C
C     * BAND INFORMATION.
C
      REAL   SFINPTL(NBL)
      INTEGER KGS(NBS), KGSGH(NBS), KGL(NBL), KGLGH(NBL)
C
      LOGICAL LCSW, LCLW
      LOGICAL GH, RADUST
C
      COMMON /ECCENT/ RRSQ,COSD,SIND
      COMMON /ITOPLW/ LEV1
C
      DATA KGS   / 6, 4, 6, 4 /
      DATA KGL   / 1, 1, 2, 5, 2, 3, 3, 6, 4 /
C
C     * FOR MAXIMUM HEIGHT ABOUT 0.005 MB (100 KM).
C
      DATA KGSGH / 3, 4, 4, 9 /
      DATA KGLGH / 5, 1, 3, 5, 4, 0, 7, 3, 6 /
C
C----------------------------------------------------------------------C
C     THIS CODE CAN BE EXTENDED TO ABOUT 100 KM, IF THE MODEL TOP LEVELC
C     IS LOWER THAN THE MAXIMUM HEIGHT, THE CALCULATION CAN BE         C
C     SIMPLIFIED WITH LESS NUMBERS OF KGSGH AND KGLGH ACCOUNTED        C
C     IF MODEL TOP LEVEL IS CLOSE TO 1 MB:                             C
C                                                                      C
C     DATA KGSGH / 3, 3, 3, 6 /                                        C
C     DATA KGLGH / 2, 1, 2, 4, 3, 0, 6, 2, 3 /                         C
C----------------------------------------------------------------------C
C
C     * FRACTION OF INCIDENT SOLAR RADIATION IN VISIBLE AND NEAR-IR
C     * SPECTRAL INTERVALS (CONSISTENT WITH THOSE USED IN NEW RADIATION)
C
      REAL SUN(NBS)
      DATA SUN /0.46288996, 0.32269499, 0.18182821, 0.032586781/

      COMMON /SWBAND/ BAND_GW(4)
C
C...  WEIGHT = SPECIFIC EXTINCTION FOR EACH BAND / THAT AT 0.55 UM
      DATA (SWEIGHT(I), I = 1, 4) /  1.0704744, 0.6936452,
     1                               0.2615633, 0.0870727  /
C
C...  SINGLE SCATTERING ALBEDO / ASYMMETRY FACTOR FOR VOCANO
C     STYRATOSPHERIC AEROSOL, ASSUMING THE MAIN COMPONENT IS H2SO4 AND
C     THE DRY SIZE 0.35 UM (DATA FROM LI ET AL 2002 JAS)
C
      DATA (VSSSCA(I), I = 1, 4) /
     1                     0.9999998,0.9999966,0.9976643,0.4911379  /
      DATA (VSCGA(I), I = 1, 4) /
     1                     0.7328027,0.7189081,0.6573516,0.5193004  /
C...  L/W DATA FOR EXPLOSIVE VOLCANOES
      DATA (VOLABS(I), I = 1, 9) /0.0573556, 0.0552949, 0.0532003,
     1                            0.1229232, 0.0880285, 0.0417163,
     2                            0.0282683, 0.0169088, 0.0/
C
      DATA RADUST /.TRUE./
C
C----------------------------------------------------------------------C
C     FOR HRCOEF, 9.80665 / 1004.64 / 100 = 9.761357E-05, IN (K / SEC),C
C     SINCE WE USE DP (DIFF IN PRESSURE) INSTEAD OF DIFF IN METER,     C
C     THERE IS A FACTOR 1.02. THUS 9.761357E-05 * 1.02 = 9.9565841E-05 C
C     UU3 = 3 * U * U, U = 1 / e0.5                                   C
C----------------------------------------------------------------------C
C
      DATA HRCOEF, UU3, CUT / 9.9565841E-05, 1.1036383, 0.001 /
C     THE VALUE FOR SCALE_SOLAR WAS TAKEN FROM LEAN'S RECOMMENDATION  C
C     IN THE CMIP5 DOCUMENTATION.                                     C
      DATA SPECIRR,SCALE_SOLAR /1367.9396,0.9965/
      DATA QMIN     /2.E-6 /
C=======================================================================
      PI                        =  3.1415926
C
C----------------------------------------------------------------------C
C     SCALE MEAN (ANNUAL) VALUE OF SOLAR CONSTANT BY RRSQ ACCOUNTING   C
C     FOR ECCENTRICITY (PASSED THROUGH COMMON BLOCK "ECCENT" - SEE     C
C     ROUTINE SDET4). THE SPECTRAL IRRADIANCE FOR MODEL IS 1367.9396   C
C     W/m2  WHICH IS THE SOLAR ENERGY CONTAINED IN THE SPECTRAL       C
C     REGION 0.2 - 10 UM (50000 - 1000 CM).                         C
C     FOR LONGWAVE, FROM BAND1 TO BAND4, THE SOLAR AND INFRARED        C
C     INTERACTION IS CONSIDERED. THE TOTAL SOLAR ENERGY CONSIDERED IN  C
C     THE INFRARED REGION IS 11.9096 W / M2. SFINPTL IS THE INPUT     C
C     SOLAR FLUX IN EACH LONGWAVE BAND                                 C
C     THE SOLAR INPUT IN SHORTWAVE REGION IS 1367.9396 - 11.9096 =     C
C     1356.0300, THE SOLAR FRACTIONS FOR EACH BAND ARE SET IN GASOPTS2 C
C----------------------------------------------------------------------C
C
      SOLARC                    =  SOLAR_C
      FRACS                     =  RRSQ * SOLARC / SPECIRR
      X                         =  FRACS / PI
C
      SFINPTL(1)                =  3.67839 * X
      SFINPTL(2)                =  2.79694 * X
      SFINPTL(3)                =  3.20284 * X
      SFINPTL(4)                =  1.13984 * X
      SFINPTL(5)                =  0.31893 * X
      SFINPTL(6)                =  0.35404 * X
      SFINPTL(7)                =  0.29578 * X
      SFINPTL(8)                =  0.99624E-01 * X
      SFINPTL(9)                =  0.23220E-01 * X
C
C----------------------------------------------------------------------C
C     INITIALIZATION                                                   C
C----------------------------------------------------------------------C
      DO 10 I = 1, 12
        SVAR(I)                 = SOLV(I) * RRSQ * SCALE_SOLAR
   10 CONTINUE
C
      DO 20 I = IL1, IL2
        PFULL(I,LEV)            =  0.01 * PRESSG(I) * SHTJ(I,LEV)
        MTOP(I)                 =  0
        ISUN(I)                 =  1
   20 CONTINUE
C------------------------------------------------------------------------C
C Create mask indicating if calculation should be performed for current  C
C tile (1) or not (0).  If ITILERAD=0 (no tiles) it defaults to only one C
C tile, if ITILERAD=1 it uses FARET to decide if calculation should      C
C be performed and if ITILERAD=2 then tiles will be randomly sampled,    C
C as in McICA.                                                           C
C------------------------------------------------------------------------C

      IF (ITILERAD .EQ. 0) THEN
! JNSC Add computation of the mean surface emissivity, albedo and skin temperature
         ITILE(IL1:IL2,1) = 1
         DO K = 1, NTILE
            DO I = IL1, IL2
               IF (K .EQ. 1) THEN
                  ITILE(I,K) = 1
                  FAREL(I,K) = 1.0
               ELSE
                  ITILE(I,K) = 0
                  FAREL(I,K) = 0.0
               END IF
               EM0TL(I,K) = EM0(I)
               GTTL(I,K)  = GT(I)
               DO IB = 1, NBS
                  SALBTL(I,K,IB) = SALB(I,IB)
                  CSALTL(I,K,IB) = CSAL(I,IB)
               END DO ! IB
            END DO ! I
         END DO ! K
      ELSEIF (ITILERAD .EQ. 1) THEN
         DO K = 1, NTILE
            DO I = IL1, IL2
               IF (FARET(I,K) .GT. 0.0) THEN
                  ITILE(I,K) = 1
               ELSE
                  ITILE(I,K) = 0
               END IF
               EM0TL(I,K) = EM0T(I,K)
               GTTL(I,K)  = GTT(I,K)
               FAREL(I,K) = FARET(I,K)
               DO IB = 1, NBS
                  SALBTL(I,K,IB) = SALBT(I,K,IB)
                  CSALTL(I,K,IB) = CSALT(I,K,IB)
               END DO ! IB
            END DO ! I
         END DO ! K
      END IF
C
      DO 30 K = 1, LAY
      DO 30 I = IL1, IL2
        X                       =  0.01 * PRESSG(I)
        P(I,K)                  =  SHJ (I,K) * X
        PFULL(I,K)              =  SHTJ(I,K) * X
C
C----------------------------------------------------------------------C
C     O2 MASS MIXING RATIO, NOT TREATED AS GHG BUT MAKE IT 3D          C
C     CONSISTENT WITH OTHER GASES                                      C
C----------------------------------------------------------------------C
C
        O2(I,K)                 =  0.2315
   30 CONTINUE

      DO IB = 1, NBS
         DO K = 1, LAY
            DO I = IL1, IL2
               VSTAU(I,K,IB) = 0.0
               VSSSA(I,K,IB) = 0.0
               VSG(I,K,IB)   = 0.0
            END DO ! I
         END DO ! K
      END DO ! IB

      DO IB = 1, NBL
         DO K = 1, LAY
            DO I = IL1, IL2
               VSABS(I,K,IB) = 0.0
            END DO ! I
         END DO ! K
      END DO ! IB

C
C----------------------------------------------------------------------C
C     DETERMINATION OF THE HIGHEST PRESSURE LEVEL FOR CONTINUUM        C
C     CALCULATIONS (> 200 MB). REUSING SPACES OF MTOP AND ISUN.        C
C----------------------------------------------------------------------C
C
      DO 35 K = 1, LEV
      DO 35 I = IL1, IL2
        IF (PFULL(I,K) .GE. 200.)                                   THEN
          MTOP(I)               =  MTOP(I) + 1
          IF (MTOP(I) .EQ. 1) ISUN(I) =  K
        ENDIF
   35 CONTINUE
C
      MCONT = LEV
C
      DO 36 I = IL1, IL2
        MCONT                   =  MIN (ISUN(I), MCONT)
   36 CONTINUE
      MCONT = MCONT - 1

C
C----------------------------------------------------------------------C
C Clip the temperature profiles as needed so they are limited to be    C
C within the allowed range.                                            C
C----------------------------------------------------------------------C
C

      DO K = 1, LAY
         DO I = IL1, IL2
            TBND(I,K) =  0.0
         END DO ! I
      END DO ! K

      DO K = 1, LEV
         DO  I = IL1, IL2

            TF_CLIP(I,K) = TFULL(I,K)

            IF (TF_CLIP(I,K) .LE. MIN_TEMP) THEN
               TF_CLIP(I,K) = MIN_TEMP
            ENDIF

            IF (TF_CLIP(I,K) .GE. MAX_TEMP) THEN
               TF_CLIP(I,K) = MAX_TEMP
            END IF

         END DO ! I
      END DO ! K

      DO K = 1, LAY
         DO  I = IL1, IL2

            T_CLIP(I,K) = T(I,K)

            IF (T_CLIP(I,K) .LE. MIN_TEMP) THEN
               T_CLIP(I,K) = MIN_TEMP
               TBND(I,K) = TBND(I,K) + 1.0
            ENDIF

            IF (T_CLIP(I,K) .GE. MAX_TEMP) THEN
               T_CLIP(I,K) = MAX_TEMP
               TBND(I,K) = TBND(I,K) + 1.0
            END IF

         END DO ! I
      END DO ! K

C
C     * AEROSOL INPUT INFORMATION.
C
      IF ( IPAM.EQ.0 ) THEN
        SO4LOAD=AERIN(:,:,1)
        SSLOAD1=AERIN(:,:,2)
        SSLOAD2=AERIN(:,:,3)
        DSLOAD1=AERIN(:,:,4)
        DSLOAD2=AERIN(:,:,5)
        BCOLOAD=AERIN(:,:,6)
        BCYLOAD=AERIN(:,:,7)
        OCOLOAD=AERIN(:,:,8)
        OCYLOAD=AERIN(:,:,9)
      ELSE
        AM3LOAD=AERIN(:,:,1)
        REAM3  =AERIN(:,:,2)
        VEAM3  =AERIN(:,:,3)
        FR1    =AERIN(:,:,4)
        FR2    =AERIN(:,:,5)
        SSLOAD =AERIN(:,:,6)
        RESS   =AERIN(:,:,7)
        VESS   =AERIN(:,:,8)
        DSLOAD =AERIN(:,:,9)
        REDS   =AERIN(:,:,10)
        VEDS   =AERIN(:,:,11)
        BCLOAD =AERIN(:,:,12)
        REBC   =AERIN(:,:,13)
        VEBC   =AERIN(:,:,14)
        OCLOAD =AERIN(:,:,15)
        REOC   =AERIN(:,:,16)
        VEOC   =AERIN(:,:,17)
      ENDIF
C
C     * CALCULATIONS RELATIVING TO EXPLOSIVE VOLCANOES.
C     * NOTE THAT IF THIS SWITCH IS NOT SET (IE IEXPLVOL=0), THEN
C     * THE INITIALIZED ZERO VALUE OF VSTAU IS USED IN THE OPTICAL
C     * DEPTH CALCULATIONS. WE DO THIS SO THAT THE CODE DOES NOT
C     * GET TOO CONFUSING WITH TOO MANY REFERENCES TO THE
C     * "IEXPLVOL" SWITCH.
C
      IF(IEXPLVOL.NE.0)                                             THEN
         IF (IVTAU .EQ. 1) THEN
         ! Sato-like dataset providing vertically integrated optical thickness
         ! at visible wavelength.  Total optical thickness is spread across
         ! each layer and optical properties at other wavelengths are
         ! determined.

            SUMSIG(IL1:IL2)=0.
            DO K = 1, LAY
               DO I = IL1, IL2
                  TROPCRIT = 0.01 * TROP(I)
                  PFULLP = 0.01 * PRESSG(I) * SHTJ(I,K+1)
                  IF (PFULL(I,K) .GE. 10. .AND.
     &                PFULLP .LE. TROPCRIT) THEN
                     SUMSIG(I) =  SUMSIG(I) + PFULLP - PFULL(I,K)
                  ENDIF
               ENDDO ! I
            ENDDO ! K

            DO IB = 1, NBS
               DO  K = 1, LAY
                  DO  I = IL1, IL2
                     TROPCRIT = 0.01 * TROP(I)
                     IF (PFULL(I,K) .GE. 10. .AND.
     &                   PFULL(I,K+1) .LE. TROPCRIT) THEN
                        VSTAU(I,K,IB) = VTAU(I)
     &                                * (PFULL(I,K+1) - PFULL(I,K))
     &                                   /MAX(SUMSIG(I), 0.01)
     &                                * SWEIGHT(IB)
                        VSSSA(I,K,IB) = VSSSCA(IB)
                        VSG(I,K,IB)   = VSCGA(IB)
                     ENDIF ! PFULL
                  END DO ! I
               END DO ! K
            END DO ! IB

            DO IB = 1, NBL
               DO  K = 1, LAY
                  DO  I = IL1, IL2
                     TROPCRIT = 0.01 * TROP(I)
                     IF (PFULL(I,K) .GE. 10. .AND.
     &                   PFULL(I,K+1) .LE. TROPCRIT) THEN
                        VSABS(I,K,IB) = VTAU(I)
     &                                * (PFULL(I,K+1) - PFULL(I,K))
     &                                   /MAX(SUMSIG(I), 0.01)
     &                                * VOLABS(IB)
                     ENDIF ! PFULL
                  END DO ! I
               END DO ! K
            END DO ! IB

         ELSE IF (IVTAU .EQ. 2) THEN

         ! ETHZ dataset that provides vertically and wavelength resolved data.
            DO IB = 1, NBS
               DO  K = 1, LAY
                  DO  I = IL1, IL2
                     VSTAU(I,K,IB) = SW_EXT_SA(I,K,IB)*DZ(I,K)*M2KM
                     VSSSA(I,K,IB) = SW_SSA_SA(I,K,IB)
                     VSG(I,K,IB)   = SW_G_SA(I,K,IB)
                  END DO ! I
               END DO ! K
            END DO ! IB

            DO IB = 1, NBL
               DO  K = 1, LAY
                  DO  I = IL1, IL2
                     VSABS(I,K,IB) = LW_ABS_SA(I,K,IB)*DZ(I,K)*M2KM
                  END DO ! I
               END DO ! K
            END DO ! IB

         ENDIF ! IVTAU
      ENDIF ! IEXPLVOL

C
C----------------------------------------------------------------------C
C DEFINE THE SPECTRAL SAMPLING FOR THE SHORTWAVE AND LONGWAVE.         C
C----------------------------------------------------------------------C
C
      CALL INITSPECSAMPL2(NSAMPLE_SW, NSAMPLE_LW, NBS, NBL,
     1                    MAXNG, IVERS, IRADFORCE, MAX_SAM)
C
C----------------------------------------------------------------------C
C     INITIALIZE THE BAND-DEPENDANT OPTICAL PROPERTY ARRAYS            C
C----------------------------------------------------------------------C
C
      IF ( ISVAOD.NE.0 ) THEN
       DO 38 I = IL1, IL2
C
        EXB1(I) = 0.0
        EXB2(I) = 0.0
        EXB3(I) = 0.0
        EXB4(I) = 0.0
        EXB5(I) = 0.0
        EXBT(I) = 0.0
        ODB1(I) = 0.0
        ODB2(I) = 0.0
        ODB3(I) = 0.0
        ODB4(I) = 0.0
        ODB5(I) = 0.0
        ODBT(I) = 0.0
        ODBV(I) = 0.0
        OFB1(I) = 0.0
        OFB2(I) = 0.0
        OFB3(I) = 0.0
        OFB4(I) = 0.0
        OFB5(I) = 0.0
        OFBT(I) = 0.0
        ABB1(I) = 0.0
        ABB2(I) = 0.0
        ABB3(I) = 0.0
        ABB4(I) = 0.0
        ABB5(I) = 0.0
        ABBT(I) = 0.0
C
        EXS1(I) = 0.0
        EXS2(I) = 0.0
        EXS3(I) = 0.0
        EXS4(I) = 0.0
        EXS5(I) = 0.0
        EXST(I) = 0.0
        ODS1(I) = 0.0
        ODS2(I) = 0.0
        ODS3(I) = 0.0
        ODS4(I) = 0.0
        ODS5(I) = 0.0
        ODST(I) = 0.0
        ODSV(I) = 0.0
        OFS1(I) = 0.0
        OFS2(I) = 0.0
        OFS3(I) = 0.0
        OFS4(I) = 0.0
        OFS5(I) = 0.0
        OFST(I) = 0.0
        ABS1(I) = 0.0
        ABS2(I) = 0.0
        ABS3(I) = 0.0
        ABS4(I) = 0.0
        ABS5(I) = 0.0
        ABST(I) = 0.0
  38   CONTINUE
      ENDIF
C
      DO 40 IB = 1, NBS
      DO 40 K  = 1, LAY
      DO 40 I  = IL1, IL2
        EXTA  (I,K,IB)          =  1.0E-20
        EXOMA (I,K,IB)          =  1.0E-20
        EXOMGA(I,K,IB)          =  1.0E-20
        FA    (I,K,IB)          =  0.0
C
        EXTB  (I,K,IB)          =  0.0
        ODFB  (I,K,IB)          =  0.0
        SSAB  (I,K,IB)          =  0.0
        ABSB  (I,K,IB)          =  0.0
  40  CONTINUE
C
      DO 45 K  = 1, LAY
      DO 45 I  = IL1, IL2
        EXTA055(I,K) = 1.0E-20
        OMA055(I,K)  = 1.0E-20
        GA055(I,K)   = 0.0
        EXTB055(I,K)          =  0.0
        ODF055(I,K)           =  0.0
        SSA055(I,K)           =  0.0
        ABS055(I,K)           =  0.0
  45  CONTINUE

      DO 60 IB = 1, NBL
      DO 60 K  = 1, LAY
      DO 60 I  = IL1, IL2
        ABSA  (I,K,IB)          =  1.0E-20
  60  CONTINUE
C
C----------------------------------------------------------------------C
C     CALCULATE THE CLOUD PARAMETERS FOR SWTRAN AND LWTRAN             C
C     REUSING INPTG, INPTMG, TAUOMGC SPACE                             C
C----------------------------------------------------------------------C
C
      CALL CLDIFM5(CLDM, TAUOMGC, CLD, A1, NCD, NCU, NBLK, NCT,
     1             NCUM, NCDM, INPTMG, CUT, MAXC,
     2             IL1, IL2, ILG, LAY, LEV)
C
C----------------------------------------------------------------------C
C     CALCULATE THE SULFATE AEROSOL OPTICAL PROPERTIES, THE AEROSOL IS C
C     HYGROSCOPIC GROWTH DEPENDENT                                     C
C----------------------------------------------------------------------C
C
      DO 65 K = 1, LAY
        KP1 = K + 1
        DO 62 I = IL1,IL2
          DP(I,K)               =  0.0102 * PRESSG(I) *
     1                            (SHTJ(I,KP1) - SHTJ(I,K))
  62    CONTINUE
  65  CONTINUE
C
      IF ( ISVAOD.NE.0 ) THEN
C
C----------------------------------------------------------------------C
C     DIAGONOSTIC RESULT OF 0.55 AND ANGSTROM PARAMETER: 1, 2, 3, 4, 5 C
C     ARE FOR SULFATE, SSALT, DUST, BC AND OC, 6 IS FOR THE TOTAL      C
C     AODS IS VERTICAL INTEGRATED AOD AT 0.55 UM                       C
C     EXTS IS VERTICAL MASS WEIGHTED MEAN EXTINCTION COEFFICIENT       C
C     ABSS IS ABSORPTANCE = AOD (1 - SSA)                              C
C----------------------------------------------------------------------C
C
        IF (IPAM .EQ. 1)                                           THEN
          DO 67 I = IL1, IL2
          DO 67 J = 1, 6
              MASS(I,J)  =  0.0
              AODS(I,J)  =  0.0
              EXTS(I,J)  =  0.0
              ABSS(I,J)  =  0.0
              AOD86(I,J) =  0.0
  67      CONTINUE
        ELSE
          DO 70 K = 1, LAY
          DO 70 I = IL1,IL2
              BOD1ROL(I,K)=0.0
              BOD2ROL(I,K)=0.0
              BOD3ROL(I,K)=0.0
              BOD4ROL(I,K)=0.0
              BOD5ROL(I,K)=0.0
              BODTROL(I,K)=0.0
C
              BXT1ROL(I,K)=0.0
              BXT2ROL(I,K)=0.0
              BXT3ROL(I,K)=0.0
              BXT4ROL(I,K)=0.0
              BXT5ROL(I,K)=0.0
              BXTTROL(I,K)=0.0
C
              BOF1ROL(I,K)=0.0
              BOF2ROL(I,K)=0.0
              BOF3ROL(I,K)=0.0
              BOF4ROL(I,K)=0.0
              BOF5ROL(I,K)=0.0
              BOFTROL(I,K)=0.0
C
              BAB1ROL(I,K)=0.0
              BAB2ROL(I,K)=0.0
              BAB3ROL(I,K)=0.0
              BAB4ROL(I,K)=0.0
              BAB5ROL(I,K)=0.0
              BABTROL(I,K)=0.0
C
              SOD1ROL(I,K)=0.0
              SOD2ROL(I,K)=0.0
              SOD3ROL(I,K)=0.0
              SOD4ROL(I,K)=0.0
              SOD5ROL(I,K)=0.0
              SODTROL(I,K)=0.0
C
              SXT1ROL(I,K)=0.0
              SXT2ROL(I,K)=0.0
              SXT3ROL(I,K)=0.0
              SXT4ROL(I,K)=0.0
              SXT5ROL(I,K)=0.0
              SXTTROL(I,K)=0.0
C
              SOF1ROL(I,K)=0.0
              SOF2ROL(I,K)=0.0
              SOF3ROL(I,K)=0.0
              SOF4ROL(I,K)=0.0
              SOF5ROL(I,K)=0.0
              SOFTROL(I,K)=0.0
C
              SAB1ROL(I,K)=0.0
              SAB2ROL(I,K)=0.0
              SAB3ROL(I,K)=0.0
              SAB4ROL(I,K)=0.0
              SAB5ROL(I,K)=0.0
              SABTROL(I,K)=0.0
  70      CONTINUE
        ENDIF
      ENDIF
C
      IF ( IPAM.EQ.0 ) THEN
        DO K = 1, LAY
        DO I = IL1, IL2
          REAM3(I,K)   =  0.166
          VEAM3(I,K)   =  0.617
          REBCY(I,K)   =  0.166
          VEBCY(I,K)   =  0.617
          REOCY(I,K)   =  0.166
          VEOCY(I,K)   =  0.617
C
          RESO4(I,K)   =  0.166
          VESO4(I,K)   =  0.617
          REBCO(I,K)   =  0.05
          VEBCO(I,K)   =  0.617
          REOCO(I,K)   =  0.1
          VEOCO(I,K)   =  0.617
C
C----------------------------------------------------------------------C
C  SCALE SO4 TO (NH)4SO4, ORGANIC AEROSOL MASS = 1.4 * OCYLOAD, FOR BLKC
C----------------------------------------------------------------------C
C
          MSO4LOAD(I,K)=  4.12 * SO4LOAD(I,K)
          MOCOLOAD(I,K)=  1.40 * OCOLOAD(I,K)
          MOCYLOAD(I,K)=  1.40 * OCYLOAD(I,K)
        ENDDO
        ENDDO
C
        CALL AEROMX3OP (EXTA, EXOMA, EXOMGA, FA, ABSA, TAUA055,
     1                  ABSA055, TAUA086, RH, MSO4LOAD, BCOLOAD,
     2                  BCYLOAD, MOCOLOAD, MOCYLOAD, REAM3, VEAM3,
     3                  RESO4, VESO4, REBCY, VEBCY, REOCY, VEOCY,
     4                  REBCO, VEBCO, REOCO, VEOCO, IL1, IL2, ILG,
     5                  LAY)
C
C       *  SULFATE/BC/OC AOD
C
        IF ( ISVAOD.NE.0 ) THEN
          DO 75 K = 1, LAY
          DO 75 I = IL1,IL2
            SOD1ROL(I,K) = TAUA055(I,K,1)
            SOD4ROL(I,K) = TAUA055(I,K,2)
            SOD5ROL(I,K) = TAUA055(I,K,3)
            SODTROL(I,K) = SODTROL(I,K) + TAUA055(I,K,4)
            SAB1ROL(I,K) = ABSA055(I,K,1)
            SAB4ROL(I,K) = ABSA055(I,K,2)
            SAB5ROL(I,K) = ABSA055(I,K,3)
            SABTROL(I,K) = SABTROL(I,K) + ABSA055(I,K,4)
  75      CONTINUE
        ENDIF
C
C----------------------------------------------------------------------C
C     CALCULATE THE SEA SALT (SS) AEROSOL OPTICAL PROPERTIES,          C
C     THE AEROSOL IS HYGROSCOPIC GROWTH DEPENDENT.                     C
C     2 MODE: 1 (FINE) AND 2 (COARSE)                                  C
C----------------------------------------------------------------------C
C
        CALL SSALTAERO5 (EXTA, EXOMA, EXOMGA, FA, ABSA,
     1                   EXTA055, OMA055, GA055,
     2                   EXTB, ODFB, SSAB, ABSB, EXTB055,ODF055,
     3                   SSA055, ABS055, RH, SSLOAD1, SSLOAD2,
     4                   IL1, IL2, ILG, LAY)
C
C       *  SEA SALT AOD
C
        IF ( ISVAOD.NE.0 ) THEN
          DO 80 K = 1, LAY
          DO 80 I = IL1,IL2
            BOD2ROL(I,K)=EXTB(I,K,1)*10000.0*
     &                   (SSLOAD1(I,K)+SSLOAD2(I,K))
            BODTROL(I,K)=BODTROL(I,K)+BOD2ROL(I,K)
            BXT2ROL(I,K)=EXTB(I,K,1)
            BOF2ROL(I,K)=ODFB(I,K,1)*10000.0*SSLOAD1(I,K)
            BOFTROL(I,K)=BOFTROL(I,K)+BOF2ROL(I,K)
            BAB2ROL(I,K)=ABSB(I,K,1)
            BABTROL(I,K)=BABTROL(I,K)+ABSB(I,K,1)
            SOD2ROL(I,K)=EXTB055(I,K)*10000.0*
     &                   (SSLOAD1(I,K)+SSLOAD2(I,K))
            SODTROL(I,K)=SODTROL(I,K)+SOD2ROL(I,K)
            SXT2ROL(I,K)=EXTB055(I,K)
            SOF2ROL(I,K)=ODF055(I,K)*10000.0*SSLOAD1(I,K)
            SOFTROL(I,K)=SOFTROL(I,K)+SOF2ROL(I,K)
            SAB2ROL(I,K)=ABS055(I,K)
            SABTROL(I,K)=SABTROL(I,K)+ABS055(I,K)
  80      CONTINUE
        ENDIF
C
C----------------------------------------------------------------------C
C     CALCULATE THE DUST (DS) AEROSOL OPTICAL PROPERTIES,              C
C     2 MODE: 1 (FINE) AND 2 (COARSE)                                  C
C----------------------------------------------------------------------C
C
        CALL DUSTAERO3(EXTA, EXOMA, EXOMGA, FA, ABSA,
     1                 EXTA055, OMA055, GA055,
     2                 EXTB, ODFB, SSAB, ABSB, EXTB055,ODF055,
     3                 SSA055, ABS055, DSLOAD1, DSLOAD2,
     4                 IL1, IL2, ILG, LAY)
C
C       *  DUST AOD
C
        IF ( ISVAOD.NE.0 ) THEN
          DO 85 K = 1, LAY
          DO 85 I = IL1,IL2
            BOD3ROL(I,K)=EXTB(I,K,1)*10000.0*
     &                   (DSLOAD1(I,K)+DSLOAD2(I,K))
            BODTROL(I,K)=BODTROL(I,K)+BOD3ROL(I,K)
            BXT3ROL(I,K)=EXTB(I,K,1)
            BOF3ROL(I,K)=ODFB(I,K,1)*10000.0*DSLOAD1(I,K)
            BOFTROL(I,K)=BOFTROL(I,K)+BOF3ROL(I,K)
            BAB3ROL(I,K)=ABSB(I,K,1)
            BABTROL(I,K)=BABTROL(I,K)+ABSB(I,K,1)
            SOD3ROL(I,K)=EXTB055(I,K)*10000.0*
     &                   (DSLOAD1(I,K)+DSLOAD2(I,K))
            SODTROL(I,K)=SODTROL(I,K)+SOD3ROL(I,K)
            SXT3ROL(I,K)=EXTB055(I,K)
            SOF3ROL(I,K)=ODF055(I,K)*10000.0*DSLOAD1(I,K)
            SOFTROL(I,K)=SOFTROL(I,K)+SOF3ROL(I,K)
            SAB3ROL(I,K)=ABS055(I,K)
            SABTROL(I,K)=SABTROL(I,K)+ABS055(I,K)
  85      CONTINUE
        ENDIF
C
C          * GATHER RESULTS FOR AEROSOL OPTICAL PROPERTIES.
C
        IF ( ISVAOD.NE.0 ) THEN
           DO 100 K = 1, LAY
           DO 100 I = IL1,IL2
             EXB1(I)=EXB1(I)+BXT1ROL(I,K)
             EXB2(I)=EXB2(I)+BXT2ROL(I,K)
             EXB3(I)=EXB3(I)+BXT3ROL(I,K)
             EXB4(I)=EXB4(I)+BXT4ROL(I,K)
             EXB5(I)=EXB5(I)+BXT5ROL(I,K)
             EXBT(I)=EXBT(I)+BXTTROL(I,K)
C
             ODB1(I)=ODB1(I)+BOD1ROL(I,K) * DP(I,K)
             ODB2(I)=ODB2(I)+BOD2ROL(I,K) * DP(I,K)
             ODB3(I)=ODB3(I)+BOD3ROL(I,K) * DP(I,K)
             ODB4(I)=ODB4(I)+BOD4ROL(I,K) * DP(I,K)
             ODB5(I)=ODB5(I)+BOD5ROL(I,K) * DP(I,K)
             ODBT(I)=ODBT(I)+BODTROL(I,K) * DP(I,K)
             ODBV(I)=ODBV(I)+VSTAU(I,K,1)
C
             OFB1(I)=OFB1(I)+BOF1ROL(I,K) * DP(I,K)
             OFB2(I)=OFB2(I)+BOF2ROL(I,K) * DP(I,K)
             OFB3(I)=OFB3(I)+BOF3ROL(I,K) * DP(I,K)
             OFB4(I)=OFB4(I)+BOF4ROL(I,K) * DP(I,K)
             OFB5(I)=OFB5(I)+BOF5ROL(I,K) * DP(I,K)
             OFBT(I)=OFBT(I)+BOFTROL(I,K) * DP(I,K)
C
             ABB1(I)=ABB1(I)+BAB1ROL(I,K) * DP(I,K)
             ABB2(I)=ABB2(I)+BAB2ROL(I,K) * DP(I,K)
             ABB3(I)=ABB3(I)+BAB3ROL(I,K) * DP(I,K)
             ABB4(I)=ABB4(I)+BAB4ROL(I,K) * DP(I,K)
             ABB5(I)=ABB5(I)+BAB5ROL(I,K) * DP(I,K)
             ABBT(I)=ABBT(I)+BABTROL(I,K) * DP(I,K)
C
             EXS1(I)=EXS1(I)+SXT1ROL(I,K)
             EXS2(I)=EXS2(I)+SXT2ROL(I,K)
             EXS3(I)=EXS3(I)+SXT3ROL(I,K)
             EXS4(I)=EXS4(I)+SXT4ROL(I,K)
             EXS5(I)=EXS5(I)+SXT5ROL(I,K)
             EXST(I)=EXST(I)+SXTTROL(I,K)
C
             ODS1(I)=ODS1(I)+SOD1ROL(I,K) * DP(I,K)
             ODS2(I)=ODS2(I)+SOD2ROL(I,K) * DP(I,K)
             ODS3(I)=ODS3(I)+SOD3ROL(I,K) * DP(I,K)
             ODS4(I)=ODS4(I)+SOD4ROL(I,K) * DP(I,K)
             ODS5(I)=ODS5(I)+SOD5ROL(I,K) * DP(I,K)
             ODST(I)=ODST(I)+SODTROL(I,K) * DP(I,K)
C
             OFS1(I)=OFS1(I)+SOF1ROL(I,K) * DP(I,K)
             OFS2(I)=OFS2(I)+SOF2ROL(I,K) * DP(I,K)
             OFS3(I)=OFS3(I)+SOF3ROL(I,K) * DP(I,K)
             OFS4(I)=OFS4(I)+SOF4ROL(I,K) * DP(I,K)
             OFS5(I)=OFS5(I)+SOF5ROL(I,K) * DP(I,K)
             OFST(I)=OFST(I)+SOFTROL(I,K) * DP(I,K)
C
             ABS1(I)=ABS1(I)+SAB1ROL(I,K) * DP(I,K)
             ABS2(I)=ABS2(I)+SAB2ROL(I,K) * DP(I,K)
             ABS3(I)=ABS3(I)+SAB3ROL(I,K) * DP(I,K)
             ABS4(I)=ABS4(I)+SAB4ROL(I,K) * DP(I,K)
             ABS5(I)=ABS5(I)+SAB5ROL(I,K) * DP(I,K)
             ABST(I)=ABST(I)+SABTROL(I,K) * DP(I,K)
 100       CONTINUE
C
           DO 110 I = IL1,IL2
             EXB1(I)=EXB1(I)/LAY
             EXB2(I)=EXB2(I)/LAY
             EXB3(I)=EXB3(I)/LAY
             EXB4(I)=EXB4(I)/LAY
             EXB5(I)=EXB5(I)/LAY
             EXBT(I)=EXBT(I)/LAY
C
             EXS1(I)=EXS1(I)/LAY
             EXS2(I)=EXS2(I)/LAY
             EXS3(I)=EXS3(I)/LAY
             EXS4(I)=EXS4(I)/LAY
             EXS5(I)=EXS5(I)/LAY
             EXST(I)=EXST(I)/LAY
C
             IF(IEXPLVOL.NE.0)THEN
               ODSV(I)=VTAU(I)
             ELSE
               ODSV(I)=0.
             ENDIF
 110       CONTINUE
        ENDIF
      ELSE
C
C----------------------------------------------------------------------C
C     CALCULATE THE AEROSOL OPTICAL PROPERTIES FOR INTERNAL MIXING OF  C
C     SULFATE, BC AND OC. AEROSOL IS HYGROSCOPIC GROWTH DEPENDENT.     C
C     CONTINOUS SIZE DISTRIBUTION                                      C
C----------------------------------------------------------------------C
C
C
C     * MASS MIXING RATIOS OF INTERNALLY MIXED AEROSOL SPECIES.
C
        DO K = 1, LAY
        DO I = IL1, IL2
          FR3          = MIN(MAX((1. - FR1(I,K) - FR2(I,K)), 0.),1.)
          SO4LOAD(I,K) = FR1(I,K) * AM3LOAD(I,K)
          BCYLOAD(I,K) = FR2(I,K) * AM3LOAD(I,K)
          OCYLOAD(I,K) = FR3 * AM3LOAD(I,K)
C
C     * MASS MIXING RATIOS OF EXTERNALLY MIXED AEROSOL SPECIES.
C
          BCOLOAD(I,K) = BCLOAD(I,K)
          OCOLOAD(I,K) = OCLOAD(I,K)
C
C     * EFFECTIVE RADIUS AND VARIANCE OF HYDROPHILIC SPECIES THAT
C     * ARE TREATED AS EXTERNALLY MIXED (UNUSED).
C
          RESO4(I,K) = REAM3(I,K)
          REBCY(I,K) = REAM3(I,K)
          REOCY(I,K) = REAM3(I,K)
          VESO4(I,K) = VEAM3(I,K)
          VEBCY(I,K) = VEAM3(I,K)
          VEOCY(I,K) = VEAM3(I,K)
        ENDDO
        ENDDO
C
        CALL AEROMX3OP (EXTA, EXOMA, EXOMGA, FA, ABSA, TAUA055, ABSA055,
     1                  TAUA086, RH, SO4LOAD, BCOLOAD, BCYLOAD, OCOLOAD,
     2                  OCYLOAD, REAM3, VEAM3, RESO4, VESO4, REBCY,
     3                  VEBCY, REOCY, VEOCY, REBC, VEBC, REOC, VEOC,
     4                  IL1, IL2, ILG, LAY)
C
        IF (ISVAOD .NE. 0)                                          THEN
          DO K = 1, LAY
          DO I = IL1, IL2
            FR3           =  MIN(MAX((1. - FR1(I,K) - FR2(I,K)),0.), 1.)
            AM3LOADS      =  10000. * AM3LOAD(I,K) * DP(I,K)
C
            AODS(I,1)     =  AODS(I,1) +  TAUA055(I,K,1) * DP(I,K)
            AODS(I,4)     =  AODS(I,4) +  TAUA055(I,K,2) * DP(I,K)
            AODS(I,5)     =  AODS(I,5) +  TAUA055(I,K,3) * DP(I,K)
            AODS(I,6)     =  AODS(I,6) +  TAUA055(I,K,4) * DP(I,K)
            AOD86(I,1)    =  AOD86(I,1) + TAUA086(I,K,1) * DP(I,K)
            AOD86(I,4)    =  AOD86(I,4) + TAUA086(I,K,2) * DP(I,K)
            AOD86(I,5)    =  AOD86(I,5) + TAUA086(I,K,3) * DP(I,K)
            AOD86(I,6)    =  AOD86(I,6) + TAUA086(I,K,4) * DP(I,K)
C
            MASS(I,1)     =  MASS(I,1) + AM3LOADS * FR1(I,K)
            MASS(I,4)     =  MASS(I,4) + AM3LOADS * FR2(I,K)
            MASS(I,5)     =  MASS(I,5) + AM3LOADS * FR3
            MASS(I,6)     =  MASS(I,6) + AM3LOADS
C
            ABSS(I,1)     =  ABSS(I,1) + ABSA055(I,K,1) * DP(I,K)
            ABSS(I,4)     =  ABSS(I,4) + ABSA055(I,K,2) * DP(I,K)
            ABSS(I,5)     =  ABSS(I,5) + ABSA055(I,K,3) * DP(I,K)
            ABSS(I,6)     =  ABSS(I,6) + ABSA055(I,K,4) * DP(I,K)
          ENDDO
          ENDDO
        ENDIF
C
C----------------------------------------------------------------------C
C     CALCULATE THE SEA SALT (SS) AEROSOL OPTICAL PROPERTIES, THE      C
C     AEROSOL IS HYGROSCOPIC GROWTH DEPENDENT, CONTINOUS SIZE          C
C     DISTRIBUTION                                                     C
C----------------------------------------------------------------------C
C
        CALL SSALTAEROP(EXTA, EXOMA, EXOMGA, FA, ABSA, EXTA055, EXTA086,
     1                  SSA055, RH, SSLOAD, RESS, VESS,
     2                  IL1, IL2, ILG, LAY)
C
        IF (ISVAOD .NE. 0)                                          THEN
          DO K = 1, LAY
          DO I = IL1, IL2
            SSLOADS       =  10000.0 * SSLOAD(I,K) * DP(I,K)
            X55           =  EXTA055(I,K) * SSLOADS
            X86           =  EXTA086(I,K) * SSLOADS
C
            AOD86(I,2)    =  AOD86(I,2) + X86
            AODS(I,2)     =  AODS(I,2) + X55
            MASS(I,2)     =  MASS(I,2) + SSLOADS
            ABSS(I,2)     =  ABSS(I,2) + X55 * (1.0 - SSA055(I,K))
          ENDDO
          ENDDO
        ENDIF
C
C----------------------------------------------------------------------C
C     CALCULATE THE DUST (DS) AEROSOL OPTICAL PROPERTIES, NO           C
C     HYGROSCOPIC EFFECT, CONTINOUS SIZE DISTRIBUTION                  C
C----------------------------------------------------------------------C
C
        CALL DUSTAEROP(EXTA, EXOMA, EXOMGA, FA, ABSA, EXTA055, EXTA086,
     1                 SSA055, DSLOAD, REDS, VEDS, IL1, IL2, ILG, LAY)
C
        IF (ISVAOD .NE. 0)                                          THEN
          AODRMN=EXP(-10.)
          AODRMX=EXP(10.)
          DO K = 1, LAY
          DO I = IL1, IL2
            DSLOADS       =  10000.0 * DSLOAD(I,K) * DP(I,K)
            X55           =  EXTA055(I,K) * DSLOADS
            X86           =  EXTA086(I,K) * DSLOADS
C
            AOD86(I,3)    =  AOD86(I,3) + X86
            AODS(I,3)     =  AODS(I,3) + X55
            MASS(I,3)     =  MASS(I,3) + DSLOADS
            ABSS(I,3)     =  ABSS(I,3) + X55 * (1.0 - SSA055(I,K))
          ENDDO
          ENDDO
C
          DO I = IL1, IL2
            AODS(I,6)     =  AODS(I,6) + AODS(I,2) + AODS(I,3)
            ABSS(I,6)     =  ABSS(I,6) + ABSS(I,2) + ABSS(I,3)
            MASS(I,6)     =  MASS(I,6) + MASS(I,2) + MASS(I,3)
            AOD86(I,6)    =  AOD86(I,6) + AOD86(I,2) + AOD86(I,3)
            DO J = 1, 6
              IF (MASS(I,J) .GT. 0.0)                               THEN
                EXTS(I,J) =  AODS(I,J) / MASS(I,J)
              ELSE
                EXTS(I,J) =  0.0
              ENDIF
            ENDDO
C----------------------------------------------------------------------C
C     ANGSTROM PARAMETER, BASED ON 0.55 AND 0.865 UM,                  C
C     2.20843 = - 1 / LOG (0.55 / 0.865)                               C
C----------------------------------------------------------------------C
C
            DO J = 1, 6
              IF (      (AODS(I,J)/AOD86(I,J).GT.AODRMN)
     1            .AND. (AODS(I,J)/AOD86(I,J).LT.AODRMX) )          THEN
                ANGSTR(I,J)   =  2.20843 * LOG (AODS(I,J) / AOD86(I,J))
              ELSE
                ANGSTR(I,J)   =  0.0
              ENDIF
            ENDDO
          ENDDO
C
C        * SAVE RESULTS IN OUTPUT ARRAYS.
C
          DO I = IL1, IL2
            ODS1(I)=AODS(I,1)
            ODS2(I)=AODS(I,2)
            ODS3(I)=AODS(I,3)
            ODS4(I)=AODS(I,4)
            ODS5(I)=AODS(I,5)
            ODST(I)=AODS(I,6)
C
C           * SAVE TOTAL AEROSOL OPTICAL DEPTH FOR DAYTIME/CLEAR CONDITIONS.
C
C            ODST(I)=0.
C            IF ( RMU(I) .GT. 0.001 .AND. CLDT(I) < 0.5 ) THEN
C              ODST(I)=AODS(I,6)
C            ENDIF
C
            EXS1(I)=EXTS(I,1)
            EXS2(I)=EXTS(I,2)
            EXS3(I)=EXTS(I,3)
            EXS4(I)=EXTS(I,4)
            EXS5(I)=EXTS(I,5)
            EXST(I)=EXTS(I,6)
C
            ABS1(I)=ABSS(I,1)
            ABS2(I)=ABSS(I,2)
            ABS3(I)=ABSS(I,3)
            ABS4(I)=ABSS(I,4)
            ABS5(I)=ABSS(I,5)
            ABST(I)=ABSS(I,6)
C
C           * SAVE TOTAL ABSORPTION AEROSOL OPTICAL DEPTH FOR
C           * DAYTIME/CLEAR CONDITIONS.
C
C            ABST(I)=0.
C            IF ( RMU(I) .GT. 0.001 .AND. CLDT(I) < 0.5 ) THEN
C              ABST(I)=ABSS(I,6)
C            ENDIF
C
            OFS1(I)=ANGSTR(I,1)
            OFS2(I)=ANGSTR(I,2)
            OFS3(I)=ANGSTR(I,3)
            OFS4(I)=ANGSTR(I,4)
            OFS5(I)=ANGSTR(I,5)
            OFST(I)=ANGSTR(I,6)
C
C           * SAVE ANGSTROM EXPONENT FOR DAYTIME/CLEAR CONDITIONS.
C
C            OFST(I)=0.
C            IF ( RMU(I) .GT. 0.001 .AND. CLDT(I) < 0.5 ) THEN
C              OFST(I)=ANGSTR(I,6)
C            ENDIF
          ENDDO
          DO I = IL1, IL2
            ODB1(I)=0.
            ODB2(I)=0.
            ODB3(I)=0.
            ODB4(I)=0.
            ODB5(I)=0.
            ODBT(I)=0.
C
            EXB1(I)=0.
            EXB2(I)=0.
            EXB3(I)=0.
            EXB4(I)=0.
            EXB5(I)=0.
            EXBT(I)=0.
C
            ABB1(I)=0.
            ABB2(I)=0.
            ABB3(I)=0.
            ABB4(I)=0.
            ABB5(I)=0.
            ABBT(I)=0.
C
            OFB1(I)=0.
            OFB2(I)=0.
            OFB3(I)=0.
            OFB4(I)=0.
            OFB5(I)=0.
            OFBT(I)=0.
          ENDDO
        ENDIF
      ENDIF
C
C----------------------------------------------------------------------C
C     DETERMINATION OF THE INTERPRETATION POINTS IN PRESSURE. INPT FOR C
C     28 REFERENCE LEVELS AND INPTM FOR 18 LEVELS                      C
C----------------------------------------------------------------------C
C
      CALL PREINTP2(INPT, INPTM, DIP, A1(1,12), P, IL1, IL2, ILG, LAY)
C
      IF (LCSW)                                                     THEN
      CSD(IL1:IL2)=0.0
      CSF(IL1:IL2)=0.0
      FSR (IL1:IL2)=0.0
      FSRC(IL1:IL2)=0.0
      CSB(IL1:IL2)=0.0
      FSG(IL1:IL2)=0.0
      FSD(IL1:IL2)=0.0
      FSF(IL1:IL2)=0.0
      FSI(IL1:IL2)=0.0
      FSV(IL1:IL2)=0.0
      PAR(IL1:IL2)=0.0
      FSAMOON(IL1:IL2)=0.0
      CSDB(IL1:IL2,1:NBS)  = 0.0
      CSFB(IL1:IL2,1:NBS)  = 0.0
      FSDB(IL1:IL2,1:NBS)  = 0.0
      FSFB(IL1:IL2,1:NBS)  = 0.0
      FSSB(IL1:IL2,1:NBS)  = 0.0
      FSSCB(IL1:IL2,1:NBS) = 0.0
      WRKA(IL1:IL2,1:NBS) = 0.0
      WRKB(IL1:IL2,1:NBS) = 0.0
      CSDT(IL1:IL2,1:NTILE)=0.0
      CSFT(IL1:IL2,1:NTILE)=0.0
      CSBT(IL1:IL2,1:NTILE)=0.0
      FSGT(IL1:IL2,1:NTILE)=0.0
      FSDT(IL1:IL2,1:NTILE)=0.0
      FSFT(IL1:IL2,1:NTILE)=0.0
      FSIT(IL1:IL2,1:NTILE)=0.0
      FSVT(IL1:IL2,1:NTILE)=0.0
C
      CSDBT (IL1:IL2,1:NTILE,1:NBS) = 0.0
      CSFBT (IL1:IL2,1:NTILE,1:NBS) = 0.0
      FSDBT (IL1:IL2,1:NTILE,1:NBS) = 0.0
      FSFBT (IL1:IL2,1:NTILE,1:NBS) = 0.0
      FSSBT (IL1:IL2,1:NTILE,1:NBS) = 0.0
      FSSCBT(IL1:IL2,1:NTILE,1:NBS) = 0.0
C
      REFLT (IL1:IL2,1:NTILE,1:2,1:LEV)=0.0
      TRANT (IL1:IL2,1:NTILE,1:2,1:LEV)=0.0

C
      DO I = IL1, IL2
        FSLO(I)                 =  11.9096 * RMU(I) * FRACS
C
C       * NOTE: No perturbation in the "overlap" region beyond 4 microns
C       *       and initialize FSO to FSLO.
C
        FSO(I)                  = FSLO(I)
      ENDDO
C
      DO K = 1, LAY
      DO I = IL1, IL2
        HRS(I,K)                =  0.0
        HRL(I,K)                =  0.0
        HRSC(I,K)               =  0.0
        HRLC(I,K)               =  0.0
        HRSH2O(I,K,1)           =  0.0
        HRSH2O(I,K,2)           =  0.0
        HRSH2O(I,K,3)           =  0.0
        HRSO3 (I,K,1)           =  0.0
        HRSO3 (I,K,2)           =  0.0
        HRSO3 (I,K,3)           =  0.0
        HRSO3 (I,K,4)           =  0.0
        HRSO3 (I,K,5)           =  0.0
        HRSO3 (I,K,6)           =  0.0
        HRSO3 (I,K,7)           =  0.0
        HRSO3 (I,K,8)           =  0.0
        HRSO3 (I,K,9)           =  0.0
        HRSCO2(I,K)             =  0.0
        HRSCH4(I,K)             =  0.0
        HRSO2(I,K,1)            =  0.0
        HRSO2(I,K,2)            =  0.0
      ENDDO
      ENDDO

      IF (IRAD_FLUX_PROFS .NE. 0) THEN
! Initialize radiative flux profiles to zero
         DO L = 1, LAY+2
            DO IL = IL1, IL2
               FSAU(IL,L) = 0.0
               FSAD(IL,L) = 0.0
               FLAU(IL,L) = 0.0
               FLAD(IL,L) = 0.0
               FSCU(IL,L) = 0.0
               FSCD(IL,L) = 0.0
               FLCU(IL,L) = 0.0
               FLCD(IL,L) = 0.0
            END DO ! IL
         END DO ! L
      END IF

C
C----------------------------------------------------------------------C
C     DETERMINE WHETHER GRID POINTS ARE IN DAYLIGHT. GATHER THE        C
C     REQUIRED FIELD FOR DAYLIGHT REGION                               C
C----------------------------------------------------------------------C
C
      JYES = 0
      DO 200 I = IL1, IL2
        IF (RMU(I) .GT. 0.001)                                      THEN
          JYES                  =  JYES + 1
          ISUN(JYES)            =  I
        ENDIF
  200 CONTINUE
      LENGATH = JYES
      IL1G=1
      IL2G=LENGATH
      IF(LENGATH.EQ.0) GO TO 499              ! SKIP UNNECESSARY SOLAR!
C
      DO 230 I = IL1G, IL2G
        J = ISUN(I)
        O3TOPG(I)               =  O3TOP(J)
C
C----------------------------------------------------------------------C
C     C1 AND C2 ARE COEFFICIENTS FOR SWTRAN                            C
C----------------------------------------------------------------------C
C
        RMUG(I)                 = (2.0 * RMU(J) + SQRT(498.5225 *
     1                             RMU(J) * RMU(J) + 1.0)) / 24.35
        C1(I)                   =  0.75 * RMUG(I)
        C2(I)                   =  2.0 * C1(I) * RMUG(I)
C
        A1G(I,1)                =  A1(J,1)
        A1G(I,2)                =  A1(J,2)
        A1G(I,3)                =  A1(J,3)
        A1G(I,4)                =  A1(J,4)
        A1G(I,5)                =  A1(J,5)
        A1G(I,6)                =  A1(J,6)
        A1G(I,7)                =  1.0 - A1G(I,1) - A1G(I,2) - A1G(I,3)
        IF (A1G(I,2) .GE . CUT)                                     THEN
          A1G(I,8)              =  A1G(I,4) / A1G(I,2)
        ELSE
          A1G(I,8)              =  0.0
        ENDIF
C
        A1G(I,9)                =  0.0
        A1G(I,10)               =  0.0
        A1G(I,11)               =  0.0
        X                       =  A1G(I,3) + A1G(I,5) + A1G(I,6)
        IF (X .GE . CUT)                                            THEN
          IF (A1G(I,1) .GE . CUT)                                   THEN
            A1G(I,9)            =  A1G(I,6) / (X * A1G(I,1))
          ENDIF
          IF (A1G(I,2) .GE . CUT)                                   THEN
            A1G(I,10)           =  A1G(I,5) / (X * A1G(I,2))
          ENDIF
          A1G(I,11)             =  A1G(I,3) / X
        ENDIF
C
        A1G(I,12)               =  A1(J,12)
        NCTG(I)                 =  NCT(J)
        FLXU(I,LEV)             =  0.0
        FLXD(I,LEV)             =  0.0
        PFULLG(I,LEV)           =  PFULL(J,LEV)
C
C----------------------------------------------------------------------C
C     USING A1(I,3) FOR RMU3                                           C
C----------------------------------------------------------------------C
C
        X                       =  1.0 - RMUG(I)
        A1(I,3)                 =  X * X * X
C
        A1(I,4)                 =  0.0
        A1(I,5)                 =  T_CLIP(J,1) - 250.0
  230 CONTINUE
C
      DO 255 K = 1, LAY
        KP1 = K + 1
        DO 250 I = IL1G, IL2G
          J = ISUN(I)
          FLXU(I,K)             =  0.0
          FLXD(I,K)             =  0.0
          PFULLG(I,K)           =  PFULL(J,K)
C
C----------------------------------------------------------------------C
C     CONVERT FROM SPECIFIC HUMIDITY TO MIXING RATIO (BOUNDED).        C
C     REUSING OMCI FOR DIPG                                            C
C----------------------------------------------------------------------C
C
          QMR                   =  Q(J,K) / (1.0 - Q(J,K))
          QG(I,K)               =  MAX (QMR, QMIN)
C
          CLDMG(I,K)            =  TAUOMGC(J,K)
          CLDG(I,K)             =  CLD(J,K)
          NBLKG(I,K)            =  NBLK (J,K)
C
          O2G(I,K)              =  O2(J,K)
          O3G(I,K)              =  O3(J,K)
          CO2G(I,K)             =  CO2(J,K)
          CH4G(I,K)             =  CH4(J,K)
          TG(I,K)               =  T_CLIP(J,K)
          DT(I,K)               =  TG(I,K) - 250.0
          PG(I,K)               =  P(J,K)
          OMCI(I,K)             =  DIP(J,K)
C
          INPTG(I,K)            =  INPT(J,K)
          INPTMG(I,K)           =  INPTM(J,K)
C
C----------------------------------------------------------------------C
C     HERE DP = DIFP / G = RHO * DZ, WHERE DIFP IS THE LAYER PRESSURE  C
C     DIFFERENCE (IN MB), G IS THE GRAVITY CONSTANT, RHO IS AIR        C
C     DENSITY, AND DZ IS LAYER THICKNESS (IN CM). THEREFORE GAS MIXING C
C     RATIO * DP = GAS MASS * DZ. OR WE CAN CALL DP AS THE AIR MASS    C
C     PATH FOR A MODEL LAYER.                                          C
C     0.0102 = 1.02 * 0.01                                             C
C     1MB = 100 PASCAL = 1000 DYNES / CM2,                            C
C     1.02 = (1000 DYNES / CM2) / (980 CM / (SECOND2)).              C
C     PRESSG, SURFACE PRESSURE IN UNIT PASCAL, SO WITH 0.01 FACTOR     C
C----------------------------------------------------------------------C
C
          DP(I,K)               =  0.0102 * PRESSG(J) *
     1                            (SHTJ(J,KP1) - SHTJ(J,K))
  250   CONTINUE
  255 CONTINUE
        DO M = 1, NTILE
           DO I = IL1G, IL2G
              J = ISUN(I)
              ITILEG(I,M) = ITILE(J,M)
           END DO ! I
        END DO ! M
C
C----------------------------------------------------------------------C
C     SOLAR: 4 BAND FOR CLOUD, AEROSOL, AND RAYLEIGH,                  C
C     20 + 15 (20) MONOCHROMATIC CALCULATIONS FOR GAS AND RADIATIVE    C
C     TRANSFER                                                         C
C                                                                      C
C     FLXU:   ALL SKY SW UPWARD FLUX.                                  C
C     FLXD:   ALL SKY SW DOWNWARD FLUX.                                C
C     FSG:    DOWNWARD FLUX ABSORBED BY GROUND.                        C
C     FSD:    DIRECT DOWNWARD FLUX AT THE SURFACE.                     C
C     FSF:    DIFFUSE DOWNWARD FLUX AT THE SURFACE.                    C
C     FSV:    VISIBLE DOWNWARD FLUX AT THE SURFACE.                    C
C     FSI:    NEAR INFRARED DOWNWARD FLUX AT THE SURFACE.              C
C     PAR:    PHOTOSYNTHETIC ACTIVE RADIATION.                         C
C     ALBPLA: PLANETARY ALBEDO.                                        C
C     FSR:    REFLECTED ALL-SKY SOLAR FLUX AT TOP.                     C
C     FSRC:   REFLECTED CLEAR-SKY SOLAR FLUX AT TOP.                   C
C     FSAMOON:ABSORBED ALL-SKY SOLAR FLUX IN MOON LAYER.               C
C     CSB:    NET CLEAR SKY FLUX AT SURFACE.                           C
C     CSD:    DIRECT CLEAR SKY FLUX AT SURFACE.                        C
C     CSF:    DIFFUSE CLEAR SKY FLUX AT SURFACE.                       C
C     CSDB:   CLEAR SKY DIRECT DOWNWARD AT THE SURFACE FOR EACH BAND.  C
C     CSFB:   CLEAR SKY DIFFUSE DOWNWARD AT THE SURFACE FOR EACH BAND. C
C     FSDB:   ALL SKY DIRECT DOWNWARD AT THE SURFACE FOR EACH BAND.    C
C     FSFB:   ALL SKY DIFFUSE DOWNWARD AT THE SURFACE FOR EACH BAND.   C
C     FSSB:   ALL SKY DOWNWARD AT THE SURFACE FOR EACH BAND.           C
C     FSSCB:  CLEAR SKY DOWNWARD AT THE SURFACE FOR EACH BAND.         C

C     FLXUT:  TILED ALL SKY SW UPWARD FLUX.                            C
C     FLXDT:  TILED ALL SKY SW DOWNWARD FLUX.                          C
C     FSGT:   TILED DOWNWARD FLUX ABSORBED BY GROUND.                  C
C     FSDT:   TILED DIRECT DOWNWARD FLUX AT THE SURFACE.               C
C     FSFT:   TILED DIFFUSE DOWNWARD FLUX AT THE SURFACE.              C
C     FSVT:   TILED VISIBLE DOWNWARD FLUX AT THE SURFACE.              C
C     FSIT:   TILED NEAR INFRARED DOWNWARD FLUX AT THE SURFACE.        C
C     CSBT:   TILED NET CLEAR SKY FLUX AT SURFACE.                     C
C     CSDT:   TILED DIRECT CLEAR SKY FLUX AT SURFACE.                  C
C     CSFT:   TILED DIFFUSE CLEAR SKY FLUX AT SURFACE.                 C
C     CSDBT:  TILED CLEAR SKY DIRECT DOWNWARD AT THE SURFACE           C
C             FOR EACH BAND.                                           C
C     CSFBT:  TILED CLEAR SKY DIFFUSE DOWNWARD AT THE SURFACE          C
C             FOR EACH BAND.                                           C
C     FSDBT:  TILED ALL SKY DIRECT DOWNWARD AT THE SURFACE             C
C             FOR EACH BAND.                                           C
C     FSFBT:  TILED ALL SKY DIFFUSE DOWNWARD AT THE SURFACE            C
C             FOR EACH BAND.                                           C
C     FSSBT:  TILED ALL SKY DOWNWARD AT THE SURFACE FOR EACH BAND.     C
C     FSSCBT: TILED CLEAR SKY DOWNWARD AT THE SURFACE FOR EACH BAND.   C
C----------------------------------------------------------------------C
C
      ICOLUMN = 1
C
C     * INITIALIZE TILED WORK FIELDS.
C
      DO 275 K = 1, LEV
      DO 275 M = 1,NTILE
      DO 275 I = IL1G, IL2G
        FLXUT(I,M,K) = 0.
        FLXDT(I,M,K) = 0.
  275 CONTINUE

      DO 480 IB = 1, NBS
C
C----------------------------------------------------------------------C
C     SCALING AEROSOL OPTICAL PROPERTIES. TAUA IS AEROSOL OPTICAL DEPTHC
C----------------------------------------------------------------------C
C
        DO K = 1, LAY
        DO I = IL1G, IL2G
          J = ISUN(I)
          X1            =  VSTAU(J,K,IB)
          TAUA(I,K)     =  EXTA(J,K,IB) * DP(I,K) + X1

          X2            =  X1 * VSSSA(J,K,IB)
          TAUOMA(I,K)   =  EXOMA(J,K,IB) * DP(I,K) + X2

          X3            =  X2 * VSG(J,K,IB)
          TAUOMGA(I,K)  =  EXOMGA(J,K,IB) * DP(I,K) + X3

          F1(I,K)       =  FA(J,K,IB) * DP(I,K) + X3 * VSG(J,K,IB)
C
        ENDDO
        ENDDO
C
C----------------------------------------------------------------------C
C     RAYLEI, NEAR-IR RAYLEIGH SCATTERING, IT IS INDEPENDENT OF IG.    C
C     REUSING A1(I,1) FOR MOON LAYER ATTENUATION                       C
C----------------------------------------------------------------------C
C
        IF (IB .NE. 1)                                              THEN
          CALL RAYLEI (TAUR, IB, DP, IL1G, IL2G, ILG, LAY)
        ENDIF
C
C
        DO 460 IGTMP = 1, KGS(IB)+KGSGH(IB)

         IF (IGTMP .LE. KGS(IB)) THEN
           GH  = .FALSE.
           IG  = IGTMP
           IGH = 1
         ELSE
           GH  = .TRUE.
           IG  = IGTMP-KGS(IB)
           IGH = 2
         END IF
C
C        * Multiple samples for this IG to reduce McICA noise.
C
         WWW = 1.0/REAL(NSAMPLE_SW(IB,IG,IGH))

         DO 450 ISAMPLE = 1, NSAMPLE_SW(IB,IG,IGH)

C
C          * RANDOMLY SELECT A SUBCOLUMN UPON WHICH TO PERFORM
C          * RADIATIVE TRANSFER FOR THIS TIMESTEP AND INTEGRATION POINT.
C
            CALL KISSVEC(ISEEDROW(1,1),ISEEDROW(1,2),
     1                   ISEEDROW(1,3),ISEEDROW(1,4),
     2                   RAN_NUM,IL1,IL2,ILG)

           DO I = IL1, IL2
              IND1 = FLOOR(RAN_NUM(I)*REAL(NCLDY(I)))+1
              IF (IND1 .GT. NCLDY(I)) IND1 = NCLDY(I)
              IF (IND1 .LT. 1)        IND1 = 1
              RAN_INDEX(I) = IND1
              IF (IRADFORCE .NE. 0) THEN ! Radiative forcing on
                 IF (IACTIVE_RT .EQ. 1) THEN ! Interactive computation
                    COLNUM_SW(I,ICOLUMN) = IND1
                 ELSE           ! Diagnostic computation
                    RAN_INDEX(I) = COLNUM_SW(I,ICOLUMN)
                 END IF
              END IF

           END DO             ! I

           ICOLUMN = ICOLUMN + 1

           DO I = IL1, IL2
              IND1 = RAN_INDEX(I)
              DO K = 1, LAY
                 LWC(I,K)     = CLW_SUB(I,K,IND1)
                 IWC(I,K)     = CIC_SUB(I,K,IND1)
                 REL(I,K)     = REL_SUB(I,K,IND1)
                 REI(I,K)     = REI_SUB(I,K,IND1)
                 IF (LWC(I,K)*DZ(I,K) .GT. CLDWATMIN .OR.
     1               IWC(I,K)*DZ(I,K) .GT. CLDWATMIN) THEN
                    CF_MCICA(I,K) = 1.0
                 ELSE
                    CF_MCICA(I,K) = 0.0
                 END IF
              END DO          ! K
           END DO             ! I
C----------------------------------------------------------------------C
C COMPUTE NCT AND MAXC FOR THE SELECTED SUBCOLUMNS.                    C
C----------------------------------------------------------------------C
           DO I = IL1,IL2
              NCT(I)       = LEV
              TMP_COUNT(I) = 0
           END DO             ! I

           MAXC = LEV

           DO K = 1, LAY
              DO I = IL1, IL2
                 IF (CF_MCICA(I,K) .GT. CUT) THEN
                    TMP_COUNT(I) = TMP_COUNT(I) + 1
                 END IF
                 IF (TMP_COUNT(I) .EQ. 1) THEN
                    NCT(I) = K
                    TMP_COUNT(I) = TMP_COUNT(I) + 1
                 END IF
              END DO          ! I
           END DO             ! K

           DO I = IL1, IL2
              MAXC =  MIN (NCT(I), MAXC)
           END DO
C
C----------------------------------------------------------------------C
C COMPUTE THE SHORTWAVE CLOUD PROPERTIES FOR THIS WAVELENGTH BAND
C AND CLOUDY SUB-COLUMN
C----------------------------------------------------------------------C
C
          CALL SW_CLD_PROPS4(CLDG, TAUCSG, TAUOMC, TAUOMGC, F2,          ! OUTPUT
     1                      WRKAG,WRKBG,                           ! OUTPUT
     2                      DZ, CF_MCICA, LWC, IWC, REL, REI,            ! INPUT
     3                      TAUA, TAUOMA, TAUOMGA, F1, ANU,              ! INPUT
     4                      CLDT, RMUG, ETA, ISUN,                       ! INPUT
     5                      MCICA, CUT, LENGATH, IL1, IL2, ILG, LAY, IB) ! CONTROL
C
          DO K = 1, NTILE
             DO I = IL1G, IL2G
                J = ISUN(I)
                ALBSURT(I,K) = SALBTL(J,K,IB)
                CSALGT (I,K) = CSALTL(J,K,IB)
             ENDDO
          END DO
C
           DO I = IL1G, IL2G
             NCTG(I)   =  NCT(ISUN(I))
           ENDDO
C
         IF (.NOT. GH) THEN ! Perform integration over these ig values

          IF (IB .EQ. 1)                                            THEN
C
C----------------------------------------------------------------------C
C     RAYLEV, VISIBLE RAYLEIGH SCATTERING, IT IS DEPENDANT ON IG.      C
C----------------------------------------------------------------------C
C
            CALL RAYLEV2(TAUR, IG, DP, A1(1,3), IL1G, IL2G, ILG, LAY)
C
C----------------------------------------------------------------------C
C     SOLAR ATTENUATION ABOVE THE MODEL TOP LAY. ONLY APPLY TO BAND    C
C     ONE FOR O3 AND O2. THIS IS TRUE ONLY FOR MODEL TOP LEVEL ABOVE   C
C     ABOUT 1 MB, WATER VAPOR CONTRIBUTION IS SMALL.                   C
C----------------------------------------------------------------------C
C
            CALL SATTENU6(A1, IB, IG, RMUG, O3TOPG, CO2G, CH4G, O2G,
     1                    PFULLG, A1G(1,12), DT, A1(1,5), INPTG, GH,
     2                    IL1G, IL2G, ILG)
          ELSE
            DO 320 I = IL1G, IL2G
              A1(I,1)           =  1.0
  320       CONTINUE
          ENDIF
C
C----------------------------------------------------------------------C
C     DOWNWARD FLUX ABOVE 1 MB, FURTHER FLUX ATTENUATION FACTOR FOR    C
C     THE LOWER REGION                                                 C
C----------------------------------------------------------------------C
C
          IF (LEV1 .GT. 1)                                          THEN
            CALL STRANDN3(TRANT, BS, A1, RMUG, DP, DT, O3G, O2G,
     1                    A1(1,3),ITILEG,
     2                    IB, IG, LEV1, IL1G, IL2G, ILG, LAY, LEV,
     3                    NTILE)
          ELSE
            DO 330 I = IL1G, IL2G
              BS(I)             =  A1(I,1)
  330       CONTINUE
          ENDIF
C
          CALL GASOPTS5(TAUG, GW, DP, IB, IG, O3G, QG, CO2G, CH4G, O2G,
     1                  INPTMG, MCONT, OMCI, DT, A1(1,3), LEV1, GH,
     2                  IL1G, IL2G, ILG, LAY)
C
          CALL SWTRAN_MCICA(REFLT, TRANT, ITILEG,
     1                      CUMDTR, BS, TAUA, TAUR, TAUG,
     2                      TAUOMA, TAUOMGA, F1, F2, TAUCSG,
     3                      TAUOMC, TAUOMGC, CLDG,
     4                      RMUG, C1, C2, ALBSURT, CSALGT, NCTG,
     5                      CUT, LEV1, IL1G, IL2G, ILG, LAY, LEV,
     6                      NTILE)
C
          IF (LEV1 .GT. 1)                                          THEN
            CALL STRANUP3(REFLT, ITILEG, DP, DT, O3G, O2G, IB, IG,
     1                    LEV1, IL1G, IL2G, ILG, LAY, LEV, NTILE)
          ENDIF
C
C----------------------------------------------------------------------C
C Compute the gridbox mean REFL and TRAN using REFLT, TRANT and the    C
C fraction covered by each tile.                                       C
C----------------------------------------------------------------------C
C
          REFL(IL1:IL2,1:2,1:LEV) = 0.0
          TRAN(IL1:IL2,1:2,1:LEV) = 0.0

          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1G, IL2G
                   J = ISUN(I)
                   IF (ITILEG(I,M) .GT. 0) THEN
                      REFL(I,1,K) = REFL(I,1,K)
     1                            + REFLT(I,M,1,K)*FAREL(J,M)
                      REFL(I,2,K) = REFL(I,2,K)
     1                            + REFLT(I,M,2,K)*FAREL(J,M)
                      TRAN(I,1,K) = TRAN(I,1,K)
     1                            + TRANT(I,M,1,K)*FAREL(J,M)
                      TRAN(I,2,K) = TRAN(I,2,K)
     1                            + TRANT(I,M,2,K)*FAREL(J,M)
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K

C         * FOR THE TOTAL SKY FLUXES WEIGHT THE CLEAR AND CLOUDY SKY
C         * FLUXES BY THE TOTAL VERTICALLY PROJECTED CLOUD FRACTION.
C
          DO K = 1, LEV
          DO I = IL1G, IL2G
            J = ISUN(I)
            IF (CLDT(J).LT.1.0) THEN
              REFL(I,2,K) = (1.0 - CLDT(J)) * REFL(I,1,K) +
     1                      CLDT(J)  * REFL(I,2,K)
              TRAN(I,2,K) = (1.0 - CLDT(J)) * TRAN(I,1,K) +
     1                      CLDT(J)  * TRAN(I,2,K)
           END IF
          END DO ! I
          END DO ! K
C
          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1G, IL2G
                   J = ISUN(I)
                   IF (ITILEG(I,M) .GT. 0) THEN
                      IF (CLDT(J).LT.1.0) THEN
                         REFLT(I,M,2,K) = (1.0-CLDT(J))*REFLT(I,M,1,K)
     1                                  + CLDT(J)*REFLT(I,M,2,K)
                         TRANT(I,M,2,K) = (1.0-CLDT(J))*TRANT(I,M,1,K)
     1                                  + CLDT(J)*TRANT(I,M,2,K)
                      END IF
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K
C
C----------------------------------------------------------------------C
C     GATHER BACK THE REQUIRED FIELDS                                  C
C----------------------------------------------------------------------C
C
          IF (IB .EQ. 1)                                            THEN
            X       = SVAR(10 - IG)
            SCALE_X = 1.0 ! Not needed for band 1
          ELSE
            X       = SVAR(8 + IB)
            SCALE_X = GW/BAND_GW(IB)
          ENDIF
          RGW = (GW*FRACS + SCALE_X*X) * WWW
          DO I = IL1, IL2
            FSO(I) = FSO(I) + RGW*RMU(I)
          END DO ! I
          DO 350 I = IL1G, IL2G
            J = ISUN(I)
            X                 =  (1.0 - CLDT(J)) * CUMDTR(I,1,LEV) +
     1                           CLDT(J) * CUMDTR(I,2,LEV)
            A1(I,2)             =  RGW * RMU(J)
            WRKA(J,IB)          =  WRKA(J,IB)+WRKAG(I)*WWW/REAL(KGS(IB))
            WRKB(J,IB)          =  WRKB(J,IB)+WRKBG(I)*WWW/REAL(KGS(IB))
            FSD(J)              =  FSD(J) + X * BS(I) * A1(I,2)
            FSR(J)          =  FSR (J) + REFL(I,2,1) * A1(I,1) * A1(I,2)
            FSRC(J)         =  FSRC(J) + REFL(I,1,1) * A1(I,1) * A1(I,2)
            CSB(J)              =  CSB(J) + (TRAN(I,1,LEV) -
     1                             REFL(I,1,LEV)) * A1(I,2)
C
            CSD(J)              =  CSD(J) +
     1                             CUMDTR(I,1,LEV) * BS(I) * A1(I,2)
            CSF(J)              =  CSF(J) + TRAN(I,1,LEV) * A1(I,2)
            FLXU(I,1)           =  FLXU(I,1) + REFL(I,2,1) * A1(I,2)
            FLXD(I,1)           =  FLXD(I,1) + TRAN(I,2,1) * A1(I,2)
C
            FSDB(J,IB)          =  FSDB(J,IB) + X * BS(I) * A1(I,2)
            FSSB(J,IB)          =  FSSB(J,IB) + TRAN(I,2,LEV)* A1(I,2)
            CSDB(J,IB)          =  CSDB(J,IB) +
     1                             CUMDTR(I,1,LEV) * BS(I) * A1(I,2)
            FSSCB(J,IB)         =  FSSCB(J,IB) + TRAN(I,1,LEV)* A1(I,2)
 350     CONTINUE

         IF (IRAD_FLUX_PROFS .NE. 0) THEN
            DO I = IL1G, IL2G
               J = ISUN(I)
               A1(I,2) =  RGW * RMU(J)
! Top of atmosphere
               FSCU(J,1) = FSCU(J,1) + REFL(I,1,1)*A1(I,1)*A1(I,2)
               FSCD(J,1) = FSCD(J,1) + A1(I,2)
               FSAU(J,1) = FSAU(J,1) + REFL(I,2,1)*A1(I,1)*A1(I,2)
               FSAD(J,1) = FSAD(J,1) + A1(I,2)
! Top of model
               FSCU(J,2) = FSCU(J,2) + REFL(I,1,1)*A1(I,2)
               FSCD(J,2) = FSCD(J,2) + TRAN(I,1,1)*A1(I,2)
            END DO ! I
         END IF ! IRAD_FLUX_PROFS

C
C----------------------------------------------------------------------C
C     HEATING RATE CALCULATION, FOR STABILITY IN CALCULATION, EACH IG  C
C     IS DONE SEPARATELY. HEATING RATE IN (K / SEC),                   C
C----------------------------------------------------------------------C
C
          DO 375 K = 1, LAY
            KP1 = K + 1
            DO 370 I = IL1G, IL2G
              J = ISUN(I)
              FLXU(I,KP1)       =  FLXU(I,KP1) + REFL(I,2,KP1) * A1(I,2)
              FLXD(I,KP1)       =  FLXD(I,KP1) + TRAN(I,2,KP1) * A1(I,2)
C
C             * FOR NON-MAM APPLICATIONS, CALCULATE HEATING RATE IN
C             * USUAL MANNER.
C
              IF(IMAM.EQ.0)                      THEN
                 DFNET          = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
                 HRS(J,K)       =  HRS(J,K) + HRCOEF * MAX(DFNET,0.) /
     1                             DP(I,K)
                 DFNETC         = (TRAN(I,1,K) - TRAN(I,1,KP1) -
     1                             REFL(I,1,K) + REFL(I,1,KP1)) *
     2                             A1(I,2)
                 HRSC(J,K)      =  HRSC(J,K) + HRCOEF * MAX(DFNETC,0.) /
     1                             DP(I,K)
              ENDIF
C
              IF (IRAD_FLUX_PROFS .NE. 0) THEN
                 FSCU(J,KP1+1) = FSCU(J,KP1+1) + REFL(I,1,KP1)*A1(I,2)
                 FSCD(J,KP1+1) = FSCD(J,KP1+1) + TRAN(I,1,KP1)*A1(I,2)
              END IF

  370       CONTINUE
  375     CONTINUE
C
          DO M = 1, NTILE
             DO I = IL1G, IL2G
                J = ISUN(I)
                IF (ITILEG(I,M) .GT. 0) THEN
                   X              = (1.0 - CLDT(J)) * CUMDTR(I,1,LEV) +
     1                              CLDT(J) * CUMDTR(I,2,LEV)
                   FLXUT(I,M,LEV) = FLXUT(I,M,LEV) + REFLT(I,M,2,LEV)
     1                            * A1(I,2)
                   FLXDT(I,M,LEV) = FLXDT(I,M,LEV) + TRANT(I,M,2,LEV) *
     1                              A1(I,2)
                   FSDT(J,M)      = FSDT(J,M) + X * BS(I) * A1(I,2)
                   CSBT(J,M)      = CSBT(J,M) + (TRANT(I,M,1,LEV) -
     1                              REFLT(I,M,1,LEV)) * A1(I,2)
                   CSDT(J,M)      = CSDT(J,M) +
     1                              CUMDTR(I,1,LEV) * BS(I) * A1(I,2)
                   CSFT(J,M)      = CSFT(J,M) + TRANT(I,M,1,LEV)*A1(I,2)

                   FSDBT (J,M,IB) = FSDBT (J,M,IB) + X * BS(I) * A1(I,2)
                   FSSBT (J,M,IB) = FSSBT (J,M,IB) + TRANT(I,M,2,LEV) *
     1                              A1(I,2)
                   CSDBT (J,M,IB) = CSDBT (J,M,IB) +
     1                              CUMDTR(I,1,LEV) * BS(I) * A1(I,2)
                   FSSCBT(J,M,IB) = FSSCBT(J,M,IB) + TRANT(I,M,1,LEV) *
     1                              A1(I,2)
                END IF
             ENDDO
          ENDDO
C
C----------------------------------------------------------------------C
C     FSAMOON IS THE ENERGY ABSORBED BETWEEN TOA AND MODEL TOP LEVEL.  C
C     A1(I,4) IS THE ADJUSTMENT FOR UPWARD FLUX FROM MODEL TOP LEVEL   C
C     TO TOA USED FOR PLANETARY ALBEDO                                 C
C----------------------------------------------------------------------C
C
          IF (IB .EQ. 1)                                            THEN
            DO 380 I = IL1G, IL2G
              J = ISUN(I)
              X                 = (1.0 - A1(I,1)) * A1(I,2)
              FSAMOON(J)        =  FSAMOON(J) + X * (1.0 + REFL(I,2,1))
              A1(I,4)           =  A1(I,4) - X * REFL(I,2,1)
  380       CONTINUE
          ENDIF
          IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     SEPARATE THE HEATING RATE FOR EACH GAS FOR MAM PROJECT           C
C     SINCE THE INTERACTION BETWEEN EACH GASES, THE HEATING RATE COULD C
C     NOT BE EXACTLY SEPARATED. FOR MAIN K INTERVALS, HRSO3 ARE        C
C     OBTAINED FOR EACH K OF VISIBLE BAND                              C
C                                                                      C
C     O3(1), O2(1),  50000 - 42000 (CM-1)                              C
C     O3(2),         42000 - 37400                                     C
C     O3(3),         37400 - 35700                                     C
C     O3(4),         35700 - 34000                                     C
C     O3(5),         34000 - 32185                                     C
C     O3(6),         32185 - 30300                                     C
C     O3(7),         30300 - 25000                                     C
C     O3(8),         25000 - 20000                                     C
C     O3(9),  O2(2), 20000 - 14500 -> H2O IS SET TO 0 ABOVE 100 MB     C
C     H2O(1), O2(2), 8400 - 14500  -> BELOW IT H2O AND O3 ARE COMBINED C
C     H2O(2), CO2,   4200 - 8400                                       C
C     H2O(3), CO2,   8400 - 2500                                       C
C     HRSCO2 + HRLCO2(1) (IN LONGWAVE) ARE 4.3 UM CO2 BAND             C
C     HRSCH4         8400 - 2500  ONLY CONSIDER THAT IN GH             C
C----------------------------------------------------------------------C
C
          IF (IB .EQ. 1)                                            THEN
            L = 10 - IG
            DO K = 1, LAY
              KP1 = K + 1
              DO I = IL1G, IL2G
                J = ISUN(I)
                DFNET           = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
                HRSO3(J,K,L)    =  HRCOEF * MAX(DFNET,0.) / DP(I,K)
              ENDDO
            ENDDO
C
            IF (IG .EQ. 1)                                          THEN
C
C----------------------------------------------------------------------C
C     SEPARATE THE O2 HEATING OUT ABOVE 1 MB                           C
C----------------------------------------------------------------------C
C
              DO I = IL1G, IL2G
                TAU             = (0.511E-04 - 0.881E-05 * A1(I,3)) *
     1                             O2(I,1) * A1G(1,12)
                TRAN(I,2,1)     =  EXP( - TAU / RMUG(I))
              ENDDO
C
              DO K = 1, LEV1 + 1
                KP1 = K + 1
                DO I = IL1G, IL2G
                  TAU           = (0.511E-04 - 0.881E-05 * A1(I,3)) *
     1                             O2(I,K) * DP(I,K)
                  TRAN(I,2,KP1) =  TRAN(I,2,K) * EXP( - TAU / RMUG(I))
                ENDDO
              ENDDO
C
              DO K = 1, LEV1 + 1
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         = (TRAN(I,2,K) - TRAN(I,2,KP1)) *
     1                             A1(I,2)
                  X             =  HRCOEF * MAX(DFNET,0.) / DP(I,K)
                  HRSO3(J,K,L)  =  HRSO3(J,K,L) - X
                  HRSO2(J,K,2)  =  X
                ENDDO
              ENDDO
            ENDIF
C
          ELSE IF (IB .EQ. 2)                                       THEN
            IF (IG .LE. 2)                                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
                  HRSO2(J,K,2)  =  HRSO2(J,K,2) +
     1                             HRCOEF * MAX(DFNET,0.) / DP(I,K)
                ENDDO
              ENDDO
            ELSE
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
C
C----------------------------------------------------------------------C
C     ZERO THE H2O ABOVE 100 MB FOR NON-LTE CONSIDERATION              C
C----------------------------------------------------------------------C
C
                  IF (PFULL(I,K) .GT. 100.0)                        THEN
                    HRSH2O(J,K,1)
     1                          =  HRSH2O(J,K,1) +
     2                             HRCOEF * MAX(DFNET,0.) / DP(I,K)
                  ENDIF
                ENDDO
              ENDDO
            ENDIF
C
          ELSE IF (IB .EQ. 3)                                       THEN
            DO K = 1, LAY
              KP1 = K + 1
              DO I = IL1G, IL2G
                J = ISUN(I)
                DFNET           = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
                HRSH2O(J,K,2)   =  HRSH2O(J,K,2) +
     1                             HRCOEF * MAX(DFNET,0.) / DP(I,K)
              ENDDO
            ENDDO
C
          ELSE IF (IB .EQ. 4)                                       THEN
            DO K = 1, LAY
              KP1 = K + 1
              DO I = IL1G, IL2G
                J = ISUN(I)
                DFNET           = (TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)) *
     2                             A1(I,2)
                HRSH2O(J,K,3)   =  HRSH2O(J,K,3) +
     1                             HRCOEF * MAX(DFNET,0.) / DP(I,K)
              ENDDO
            ENDDO
          ENDIF
          ENDIF
C
         ELSEIF (GH) THEN ! Perform integration over these ig points
C
C----------------------------------------------------------------------C
C     IN ACCUMULATED SPACE WITH INTERVAL CLOSE TO 1, THE EXTINCTION    C
C     COEFFICIENTS IS EXTREMELY LARGE, THE CALCULATION PROCESS CAN BE  C
C     SIMPLIFIED BY IGNORING SCATTERING, REFLECTION, CLOUD AND AEROSOL.C
C----------------------------------------------------------------------C
C
          CALL SATTENU6(A1, IB, IG, RMUG, O3TOPG, CO2G, CH4G, O2G,
     1                  PFULLG, A1G(1,12), DT, A1(1,5), INPTG, GH,
     2                  IL1G, IL2G, ILG)
C
          CALL STRANDNGH4(TRANT, ITILEG,
     1                    GWGH, A1, TAUA, TAUOMA, TAUCSG, TAUOMC,
     2                    CLDG, RMUG, DP, O3G, QG, CO2G, CH4G, O2G, IB,
     3                    IG, INPTG, OMCI, DT, LEV1, GH, CUT,
     4                    IL1G, IL2G, ILG, LAY, LEV, NTILE)
C
C----------------------------------------------------------------------C
C Compute the gridbox mean REFL and TRAN using REFLT, TRANT and the    C
C fraction covered by each tile.                                       C
C----------------------------------------------------------------------C
C
          REFL(IL1:IL2,1:2,1:LEV) = 0.0
          TRAN(IL1:IL2,1:2,1:LEV) = 0.0

          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1G, IL2G
                   J = ISUN(I)
                   IF (ITILEG(I,M) .GT. 0) THEN
                      TRAN(I,1,K) = TRAN(I,1,K)
     1                            + TRANT(I,M,1,K)*FAREL(J,M)
                      TRAN(I,2,K) = TRAN(I,2,K)
     1                            + TRANT(I,M,2,K)*FAREL(J,M)
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K

C
C         * FOR THE TOTAL SKY FLUXES WEIGHT THE CLEAR AND CLOUDY SKY
C         * FLUXES BY THE TOTAL VERTICALLY PROJECTED CLOUD FRACTION.
C
          DO K = 1, LEV
          DO I = IL1G, IL2G
            J = ISUN(I)
            IF (CLDT(J).LT.1.0) THEN
              TRAN(I,2,K) = (1.0 - CLDT(J)) * TRAN(I,1,K) +
     1                      CLDT(J)  * TRAN(I,2,K)
            END IF
          END DO ! I
          END DO ! K
C
          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1G, IL2G
                   J = ISUN(I)
                   IF (ITILEG(I,M) .GT. 0) THEN
                      IF (CLDT(J).LT.1.0) THEN
                         TRANT(I,M,2,K) = (1.0-CLDT(J))*TRANT(I,M,1,K)
     1                                  + CLDT(J) * TRANT(I,M,2,K)
                      END IF
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K

          IF (IB .EQ. 1)                                            THEN
            X       =  SVAR(4 - IG)
            SCALE_X = 1.0 ! Not needed for band 1
          ELSE
            X       =  SVAR(8 + IB)
            SCALE_X = GWGH/BAND_GW(IB)
          ENDIF
          RGW = (GWGH*FRACS + SCALE_X*X) * WWW
          DO I = IL1, IL2
            FSO(I) = FSO(I) + RGW*RMU(I)
          END DO ! I
C
          DO 430 I = IL1G, IL2G
            J = ISUN(I)
            A1(I,2)             =  RGW * RMU(J)
            CSB(J)              =  CSB(J) + TRAN(I,1,LEV) * A1(I,2)
C
            FSAMOON(J)          =  FSAMOON(J) +
     1                             A1(I,2) * (1.0 - TRAN(I,2,1))
            FLXD(I,1)           =  FLXD(I,1) + A1(I,2) * TRAN(I,2,1)
C
 430     CONTINUE

         IF (IRAD_FLUX_PROFS .NE. 0) THEN
            DO I = IL1G, IL2G
               J = ISUN(I)
               A1(I,2) =  RGW * RMU(J)
! Top of atmosphere
               FSCD(J,1) = FSCD(J,1) + A1(I,2)
               FSAD(J,1) = FSAD(J,1) + A1(I,2)
! Top of model
               FSCD(J,2) = FSCD(J,2) + TRAN(I,1,1)*A1(I,2)
            END DO ! I
         END IF

          DO 445 K = 1, LAY
            KP1 = K + 1
            DO 440 I = IL1G, IL2G
              J = ISUN(I)
              FLXD(I,KP1)       =  FLXD(I,KP1) + A1(I,2) * TRAN(I,2,KP1)
C
C         * FOR NON-MAM APPLICATIONS, CALCULATE HEATING RATE IN
C         * USUAL MANNER.
C
          IF(IMAM.EQ.0)                      THEN
              DFNET             =  TRAN(I,2,K) - TRAN(I,2,KP1)
              HRS(J,K)          =  HRS(J,K) + HRCOEF * A1(I,2) *
     1                             MAX(DFNET,0.) / DP(I,K)
              DFNETC            =  TRAN(I,1,K) - TRAN(I,1,KP1)
              HRSC(J,K)         =  HRSC(J,K) + HRCOEF * A1(I,2) *
     1                             MAX(DFNETC,0.) / DP(I,K)
          ENDIF
C
          IF (IRAD_FLUX_PROFS .NE. 0) THEN
             FSCD(J,KP1+1) = FSCD(J,KP1+1) + TRAN(I,1,KP1)*A1(I,2)
          END IF

  440       CONTINUE
  445     CONTINUE
C
      DO M = 1, NTILE
         DO I = IL1G, IL2G
            J = ISUN(I)
            IF (ITILEG(I,M) .GT. 0) THEN
               FLXDT(I,M,LEV) =  FLXDT(I,M,LEV) + A1(I,2) *
     1                              TRANT(I,M,2,LEV)
               CSBT(J,M)      =  CSBT(J,M) + TRANT(I,M,1,LEV)*A1(I,2)
            END IF
         ENDDO ! I
      ENDDO ! M
          IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     CALCULATIONS OF HEATING RATES FOR EACH SEPARATE GASES            C
C----------------------------------------------------------------------C
C
          IF (IB .EQ. 1)                                            THEN
            L = 4 - IG
            DO K = 1, LAY
              KP1 = K + 1
              DO I = IL1G, IL2G
                J = ISUN(I)
                DFNET           =  TRAN(I,2,K) - TRAN(I,2,KP1)
                HRSO3(J,K,L)    =  HRCOEF * A1(I,2) * DFNET / DP(I,K)
              ENDDO
            ENDDO
C
            IF (IG .EQ. 3)                                          THEN
C
C----------------------------------------------------------------------C
C     AT IB = 1, IG = 3, O2 AND O3 HEATING RATE ARE COMBINED, FOLLOWINGC
C     CALCULATION FOR O3 HEATING RATE ONLY, THEN USING (O3 + O2) - O3  C
C     TO GET O2 HEATING, THE RESULTS FOR EACH GAS COULD NOT BE AACUATE C
C     BELOW 1 MB SINCE THE STRONG OVERLPAP OF TWO GASES                C
C----------------------------------------------------------------------C
C
              DO I = IL1G, IL2G
                TAU             =  1.02 * (47069.96875 * O3TOPG(I) +
     1                             0.792006E-01 * O2(I,1)) * A1G(I,12)
                TRAN(I,2,1)     =  EXP( - TAU / RMUG(I))
              ENDDO
C
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  TAU           =  47069.96875 * O3G(I,K) * DP(I,K) +
     1                             TAUA(I,K)
                  DTR1          =  EXP( - (TAU - TAUOMA(I,K)) / RMUG(I))
C
                  IF (CLDG(I,K) .LT. CUT)                           THEN
                    TRAN(I,2,KP1)
     1                          =  TRAN(I,2,K) * DTR1
                  ELSE
                    ABSC        = (1.0 - CLDG(I,K)) * DTR1 + CLDG(I,K) *
     1                             EXP( - (TAU + TAUCSG(I,K) -
     2                             TAUOMC(I,K)) / RMUG(I))
                    TRAN(I,2,KP1)
     1                          =  TRAN(I,2,K) * ABSC
                  ENDIF
                ENDDO
              ENDDO
C
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  X             =  HRCOEF * A1(I,2) * DFNET / DP(I,K)
                  HRSO2(J,K,1)  =  HRSO3(J,K,L) - X
                  HRSO3(J,K,L)  =  X
                ENDDO
              ENDDO
            ENDIF
C
          ELSE IF (IB .EQ. 2)                                       THEN
            IF(IG .EQ. 1)                                           THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
C
C----------------------------------------------------------------------C
C     ZERO THE H2O ABOVE 100 MB FOR NON-LTE CONSIDERATION              C
C----------------------------------------------------------------------C
C
                  IF (PFULL(I,K) .GT. 100.)                         THEN
                    DFNET       =  TRAN(I,2,K) - TRAN(I,2,KP1)
                    HRSH2O(J,K,1)
     1                          =  HRSH2O(J,K,1) + HRCOEF * A1(I,2) *
     2                             DFNET / DP(I,K)
                  ENDIF
                ENDDO
              ENDDO
            ELSE
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSO2(J,K,2)  =  HRSO2(J,K,2) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ENDIF
C
          ELSE IF (IB .EQ. 3)                                       THEN
            IF (IG .LE. 2)                                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSH2O(J,K,2) =  HRSH2O(J,K,2) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ELSE
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSCO2(J,K)   =  HRSCO2(J,K) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ENDIF
C
          ELSE IF (IB .EQ. 4)                                       THEN
            IF (IG .EQ. 6 .OR. IG .EQ. 8)                           THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSH2O(J,K,3) =  HRSH2O(J,K,3) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ELSE IF (IG .EQ. 4)                                     THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSCH4(J,K)   =  HRSCH4(J,K) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ELSE
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1G, IL2G
                  J = ISUN(I)
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1)
                  HRSCO2(J,K)   =  HRSCO2(J,K) + HRCOEF * A1(I,2) *
     1                             DFNET / DP(I,K)
                ENDDO
              ENDDO
            ENDIF
          ENDIF
          ENDIF

C
         END IF !GH
C
  450   CONTINUE
C
        IF (IB .EQ. 1 .AND. IGTMP .LE. 2) THEN
           DO 455 I = IL1G, IL2G
              PAR(ISUN(I))      =  PAR(ISUN(I)) + FLXD(I,LEV)
  455      CONTINUE
        ENDIF
  460  CONTINUE
C
       IF (IB .EQ. 1) THEN
         DO 470 I = IL1G, IL2G
            FSV(ISUN(I))        =  FLXD(I,LEV)
  470    CONTINUE
       ENDIF
C
       IF (IB .EQ. 1) THEN
          DO M = 1, NTILE
             DO I = IL1G, IL2G
                J = ISUN(I)
                IF (ITILEG(I,M) .GT. 0) THEN
                   FSVT(J,M) =  FLXDT(I,M,LEV)
                END IF
             END DO
          END DO
       ENDIF
  480 CONTINUE
      ZERO = 0.0
      DO IB = 1, NBS
         DO I = IL1G, IL2G
            J = ISUN(I)
            FSFB(J,IB) = MAX(FSSB (J,IB)-FSDB(J,IB),ZERO)
            CSFB(J,IB) = MAX(FSSCB(J,IB)-CSDB(J,IB),ZERO)
            IF (FSFB(J,IB) .LT. 0.0 .OR. CSFB(J,IB) .LT. 0.0) THEN
               WRITE(6,*) J,IB,FSFB(J,IB),FSSB(J,IB),FSDB(J,IB),
     1                         CSFB(J,IB),FSSCB(J,IB),CSDB(J,IB)
               CALL XIT('RADDRIV10',-10)
            END IF
         END DO ! I
C
      DO M = 1, NTILE
         DO I = IL1G, IL2G
            J = ISUN(I)
            IF (ITILEG(I,M) .GT. 0) THEN
               FSFBT(J,M,IB) = MAX(FSSBT (J,M,IB)-FSDBT(J,M,IB),ZERO)
               CSFBT(J,M,IB) = MAX(FSSCBT(J,M,IB)-CSDBT(J,M,IB),ZERO)
            END IF
         END DO ! I
      END DO ! M
      END DO ! IB
C
C----------------------------------------------------------------------C
C     GATHER BACK REQUIRED FIELD. FOR PLANETARY ALBEDO THE INCOMING    C
C     ENERGY OF 11.9096 * FRACS IS TOTALLY ABSORBED IN LONGWAVE PART   C
C----------------------------------------------------------------------C
C
      RSOLARC = RRSQ * SOLARC
      DO 490 I = IL1G, IL2G
        J = ISUN(I)
        FSG(J)                  =  FLXD(I,LEV) - FLXU(I,LEV)
        FSI(J)                  =  FLXD(I,LEV) - FSV(J)
        FSF(J)                  =  FLXD(I,LEV) - FSD(J)
C
        CSF(J)                  =  CSF(J) - CSD(J)
  490 CONTINUE
      DO 495 M = 1, NTILE
      DO 495 I = IL1G, IL2G
         J = ISUN(I)
         IF (ITILEG(I,M) .GT. 0) THEN
            FSGT(J,M)               =  FLXDT(I,M,LEV) - FLXUT(I,M,LEV)
            FSIT(J,M)               =  FLXDT(I,M,LEV) - FSVT(J,M)
            FSFT(J,M)               =  FLXDT(I,M,LEV) - FSDT(J,M)
            CSFT(J,M)               =  CSFT(J,M) - CSDT(J,M)
         END IF
  495     CONTINUE
C
      IF (IRAD_FLUX_PROFS .NE. 0) THEN
        DO K = 1, LEV
           DO I = IL1G, IL2G
              J = ISUN(I)
              FSAU(J,K+1) = FLXU(I,K)
              FSAD(J,K+1) = FLXD(I,K)
           END DO ! I
        END DO ! K
      END IF

      IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     ADD THE HEATING RATE FOR EACH GAS TOGETHER                       C
C----------------------------------------------------------------------C
C
         DO K = 1, LAY
            DO I = IL1G, IL2G
               J = ISUN(I)
               HRS(J,K) =  HRSH2O(J,K,1) + HRSH2O(J,K,2)
     1                  +  HRSH2O(J,K,3)
     2                  +  HRSO3 (J,K,1) + HRSO3 (J,K,2)
     3                  +  HRSO3 (J,K,3) + HRSO3 (J,K,4)
     4                  +  HRSO3 (J,K,5) + HRSO3 (J,K,6)
     5                  +  HRSO3 (J,K,7) + HRSO3 (J,K,8)
     6                  +  HRSO3 (J,K,9)
     7                  +  HRSCO2(J,K)   + HRSCH4(J,K)
     8                  +  HRSO2(J,K,1)  + HRSO2(J,K,2)

               HRSC(J,K) = HRS(J,K)
            ENDDO
         ENDDO
      ENDIF ! IMAM

  499 CONTINUE
      ENDIF
C     (LCSW)
C
C----------------------------------------------------------------------C
C     LONGWAVE: 9 BAND FOR CLOUD, AEROSOL, CONTINUUM, AND PLANCK.      C
C     24+22 MONOCHROMATIC CALCULATIONS FOR GAS AND RADIATIVE TRANSFER  C
C                                                                      C
C     FLXU: ALL SKY LW UPWARD FLUX.                                    C
C     FLXD: ALL SKY LW DOWNWARD FLUX.                                  C
C     FDL:  DOWN LW FLUX RECEIVED AT THE GROUND.                       C
C     FDLC: DOWN CLEAR-SKY LW FLUX RECEIVED AT THE GROUND.             C
C     FLG:  NET  LW FLUX RECEIVED AT THE GROUND.                       C
C     OLR:    ALL-SKY UPWARD L/W FLUX AT THE TOP.                      C
C     OLRC:   CLEAR-SKY UPWARD L/W FLUX AT THE TOP.                    C
C     FLAMOON:ABSORBED ALL-SKY LW FLUX IN MOON LAYER.                  C
C     CLB:  NET CLEAR SKY DOWNWARD FLUX AT THE SURFACE.                C
C
C     FLXUT: TILED ALL SKY LW UPWARD FLUX.                             C
C     FLXTD: TILED ALL SKY LW DOWNWARD FLUX.                           C
C     FDLT:  TILED DOWN LW FLUX RECEIVED AT THE GROUND.                C
C     FDLCT: TILED DOWN CLEAR-SKY LW FLUX RECEIVED AT THE GROUND.      C
C     FLGT:  TILED NET  LW FLUX RECEIVED AT THE GROUND.                C
C     CLBT:  TILED NET CLEAR SKY DOWNWARD FLUX AT THE SURFACE.         C
C----------------------------------------------------------------------C
C
      IF (LCLW)                                                     THEN
C
C----------------------------------------------------------------------C
C     CONVERT FROM SPECIFIC HUMIDITY TO MIXING RATIO (BOUNDED) AND     C
C     BOUND TEMPERATURE FOR PLANCK CALCULATION.                        C
C----------------------------------------------------------------------C
C
      DO 505 K = 1, LAY
        KP1 = K + 1
        DO 500 I = IL1, IL2
C
C----------------------------------------------------------------------C
C     CONVERT FROM SPECIFIC HUMIDITY TO MIXING RATIO.                  C
C----------------------------------------------------------------------C
C
          QMR                   =  Q(I,K) / (1.0 - Q(I,K))
          QG(I,K)               =  MAX (QMR, QMIN)
C
          HRL(I,K)              =  0.0
          HRLC(I,K)             =  0.0
          HRLH2O(I,K,1)         =  0.0
          HRLH2O(I,K,2)         =  0.0
          HRLH2O(I,K,3)         =  0.0
          HRLH2O(I,K,4)         =  0.0
          HRLH2O(I,K,5)         =  0.0
          HRLH2O(I,K,6)         =  0.0
          HRLH2O(I,K,7)         =  0.0
          HRLH2O(I,K,8)         =  0.0
          HRLH2O(I,K,9)         =  0.0
          HRLO3(I,K)            =  0.0
          HRLCO2(I,K,1)         =  0.0
          HRLCO2(I,K,2)         =  0.0
          HRLCH4(I,K)           =  0.0
          HRLN2O(I,K)           =  0.0

          DP(I,K)               =  0.0102 * PRESSG(I) *
     1                            (SHTJ(I,KP1) - SHTJ(I,K))
          DT(I,K)               =  T_CLIP(I,K) - 250.0
  500   CONTINUE
  505 CONTINUE
C
      DO 510 I = IL1, IL2
        OLR(I)                  =  0.0
        OLRC(I)                 =  0.0
        TOA_OLR(I,:)            =  0.0
        FLAMOON(I)              =  0.0
        CLB(I)                  =  0.0
        FLG(I)                  =  0.0
        FDLC(I)                 =  0.0
        MTOP(I)                 =  0
        ISUN(I)                 =  1
        A1(I,5)                 =  T_CLIP(I,1) - 250.0
  510 CONTINUE
      DO 520 K = 1, LEV
      DO 520 I = IL1, IL2
        FLXU(I,K)               =  0.0
        FLXD(I,K)               =  0.0
  520 CONTINUE
C
      DO 530 M = 1, NTILE
      DO 530 I = IL1, IL2
        CLBT (I,M)              =  0.0
        FDLCT(I,M)              =  0.0
        FLGT (I,M)              =  0.0
  530 CONTINUE
C
C     * RE-INITIALIZE TILED WORK FIELDS.
C
      DO 550 K = 1, LEV
      DO 550 M = 1, NTILE
      DO 550 I = IL1, IL2
        FLXUT(I,M,K) = 0.
        FLXDT(I,M,K) = 0.
  550 CONTINUE
C
C----------------------------------------------------------------------C
C     DETERMINATION OF THE INTERPRETATION POINTS IN THE RATIO OF CO2   C
C     TO WATER VAPOR FOR TLINEHC. REUSE THE SPACE OF PG FOR DIR        C
C     AND REUSE TAUOMC AS A WORK ARRAY                                 C
C----------------------------------------------------------------------C
C
      CALL PREINTR3 (INPR, PG, QG, CO2, TAUOMC, IL1, IL2, ILG, LAY)
C
      ICOLUMN = 1

      DO 900 IB = 1, NBL
C
C----------------------------------------------------------------------C
C     USING C1 SPACE FOR SLWF WHICH IS THE INPUT SOLAR ENERGY IN THE   C
C     INFRARED REGION. TOTAL 11.9096 W / M2 FROM STANDARD             C
C     CALCULATION                                                      C
C     SCALING CLOUD OPTICAL PROPERTIES FOR IR SCATTERING CALCULATION   C
C----------------------------------------------------------------------C
C
        DO 605 I = IL1, IL2
          IF (RMU(I) .GT. 0.0)                                      THEN
            C1(I)               =  RMU(I) * SFINPTL(IB)
          ELSE
            C1(I)               =  0.0
          ENDIF
 605    CONTINUE
C
C       * CALCULATE AEROSOL OPTICAL DEPTH INCLUDING EXPLOSIVE VOLCANOES.
C
        DO K = 1, LAY
        DO I = IL1, IL2
          TAUA(I,K)             =  ABSA(I,K,IB) * DP(I,K) +
     1                             VSABS(I,K,IB)
        ENDDO
        ENDDO
C
C----------------------------------------------------------------------C
C    REUSING SPACE O3G FOR DBF                                         C
C----------------------------------------------------------------------C
C
        CALL PLANCK3(BF, BST, URBF, A1(1,2), A1(1,3), O3G, TF_CLIP,GTTL,
     1               IB, ITILE,
     2               IL1, IL2, ILG, LAY, LEV, NTILE)
C

        DO 875 IGTMP = 1, KGL(IB)+KGLGH(IB)

         IF (IGTMP .LE. KGL(IB)) THEN
           GH  = .FALSE.
           IG  = IGTMP
           IGH = 1
         ELSE
           GH  = .TRUE.
           IG  = IGTMP-KGL(IB)
           IGH = 2
         ENDIF

         WWW = 1.0/REAL(NSAMPLE_LW(IB,IG,IGH))

         DO 850 ISAMPLE = 1, NSAMPLE_LW(IB,IG,IGH)

C
C          * RANDOMLY SELECT A SUBCOLUMN UPON WHICH TO PERFORM
C          * RADIATIVE TRANSFER FOR THIS TIMESTEP AND INTEGRATION POINT.
C
            CALL KISSVEC(ISEEDROW(1,1),ISEEDROW(1,2),
     1                   ISEEDROW(1,3),ISEEDROW(1,4),
     2                   RAN_NUM,IL1,IL2,ILG)

           DO I = IL1, IL2
              IND1 = FLOOR(RAN_NUM(I)*REAL(NCLDY(I)))+1
              IF (IND1 .GT. NCLDY(I)) IND1 = NCLDY(I)
              IF (IND1 .LT. 1)        IND1 = 1
              RAN_INDEX(I) = IND1
              IF (IRADFORCE .NE. 0) THEN ! Radiative forcing on
                 IF (IACTIVE_RT .EQ. 1) THEN ! Interactive computation
                    COLNUM_LW(I,ICOLUMN) = IND1
                 ELSE           ! Diagnostic computation
                    RAN_INDEX(I) = COLNUM_LW(I,ICOLUMN)
                 END IF
              END IF

           END DO            ! I
      ICOLUMN = ICOLUMN + 1


           DO I = IL1, IL2
              IND1 = RAN_INDEX(I)
              DO K = 1, LAY
                 LWC(I,K) = CLW_SUB(I,K,IND1)
                 IWC(I,K) = CIC_SUB(I,K,IND1)
                 REL(I,K) = REL_SUB(I,K,IND1)
                 REI(I,K) = REI_SUB(I,K,IND1)
                 IF (LWC(I,K)*DZ(I,K) .GT. CLDWATMIN .OR.
     1               IWC(I,K)*DZ(I,K) .GT. CLDWATMIN) THEN
                    CF_MCICA(I,K) = 1.0
                 ELSE
                    CF_MCICA(I,K) = 0.0
                 END IF
              END DO        ! K
           END DO           ! I
C----------------------------------------------------------------------C
C COMPUTE NCT AND MAXC FOR THE SELECTED SUBCOLUMNS.                    C
C----------------------------------------------------------------------C
           DO I = IL1,IL2
              NCT(I)       = LEV
              TMP_COUNT(I) = 0
           END DO            ! I

           MAXC = LEV

           DO K = 1, LAY
              DO I = IL1, IL2
                 IF (CF_MCICA(I,K) .GT. CUT) THEN
                    TMP_COUNT(I) = TMP_COUNT(I) + 1
                 END IF
                 IF (TMP_COUNT(I) .EQ. 1) THEN
                    NCT(I) = K
                    TMP_COUNT(I) = TMP_COUNT(I) + 1
                 END IF
              END DO         ! I
           END DO            ! K

           DO I = IL1, IL2
              MAXC =  MIN (NCT(I), MAXC)
           END DO
C
C          * COMPUTE THE LONGWAVE CLOUD PROPERTIES FOR THIS
C          * WAVELENGTH BAND.
C
           CALL LW_CLD_PROPS3(CLDG, TAUCI, OMCI, GCI, F2,       ! OUTPUT
     1                        DZ, CF_MCICA, LWC, IWC, REL, REI, ! INPUT
     2                        CUT, UU3, IL1, IL2, ILG, LAY, IB) ! CONTROL
C
          IF (.NOT. GH) THEN      ! Perform integration over these ig values
C
           CALL GASOPTL7(TAUG, GW, DP, IB, IG, O3, QG, CO2, CH4, AN2O,
     1                   F11, F12, F113, F114, INPR, INPTM, MCONT, PG,
     2                   DIP, DT, LEV1, GH, IL1, IL2, ILG, LAY)
C
           CALL LWTRAN_MCICA2(REFLT, TRANT,
     1                        C1, TAUCI, OMCI, GCI, F2, TAUA,
     2                        TAUG, BF, BST, URBF, O3G, EM0TL,
     3                        CLDG, NCT, LEV1, CUT, MAXC,
     4                        ITILE,
     5                        IL1, IL2, ILG, LAY, LEV,
     6                        NTILE)

C
C----------------------------------------------------------------------C
C Compute the gridbox mean REFL and TRAN using REFLT, TRANT and the    C
C fraction covered by each tile.                                       C
C----------------------------------------------------------------------C
C
          REFL(IL1:IL2,1:2,1:LEV) = 0.0
          TRAN(IL1:IL2,1:2,1:LEV) = 0.0
          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1, IL2
                   IF (ITILE(I,M) .GT. 0) THEN
                      TRAN(I,1,K) = TRAN(I,1,K)
     1                            + TRANT(I,M,1,K)*FAREL(I,M)
                      TRAN(I,2,K) = TRAN(I,2,K)
     1                            + TRANT(I,M,2,K)*FAREL(I,M)
                      REFL(I,1,K) = REFL(I,1,K)
     1                            + REFLT(I,M,1,K)*FAREL(I,M)
                      REFL(I,2,K) = REFL(I,2,K)
     1                            + REFLT(I,M,2,K)*FAREL(I,M)
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K

C
C          * FOR THE TOTAL SKY FLUXES WEIGHT THE CLEAR AND CLOUDY SKY
C          * FLUXES BY THE TOTAL VERTICALLY PROJECTED CLOUD FRACTION.
C
           DO K = LEV1,LEV
           DO I = IL1,IL2
             IF (CLDT(I).LT.1.0) THEN
               REFL(I,2,K) = (1.0 - CLDT(I)) * REFL(I,1,K) +
     1                        CLDT(I)  * REFL(I,2,K)
               TRAN(I,2,K) = (1.0 - CLDT(I)) * TRAN(I,1,K) +
     1                        CLDT(I)  * TRAN(I,2,K)
             END IF
           END DO ! I
           END DO ! K
C
           DO K = LEV1, LEV
           DO M = 1, NTILE
           DO I = IL1, IL2
              IF (ITILE(I,M) .GT. 0) THEN
                 IF (CLDT(I).LT.1.0) THEN
                    REFLT(I,M,2,K) = (1.0 - CLDT(I)) * REFLT(I,M,1,K) +
     1                               CLDT(I)  * REFLT(I,M,2,K)
                    TRANT(I,M,2,K) = (1.0 - CLDT(I)) * TRANT(I,M,1,K) +
     1                                CLDT(I)  * TRANT(I,M,2,K)
                 END IF
              END IF
           END DO ! I
           END DO ! M
           END DO ! K

           PGW = PI * GW * WWW

           DO 650 K = LEV1, LAY
            KP1 = K + 1
            DO 600 I = IL1, IL2
              FLXU(I,K)         =  FLXU(I,K) + REFL(I,2,K) * PGW
              FLXD(I,K)         =  FLXD(I,K) + TRAN(I,2,K) * PGW
C
C         * FOR NON-MAM APPLICATIONS, CALCULATE HEATING RATE IN
C         * USUAL MANNER.
C
          IF(IMAM.EQ.0)                      THEN
              DFNET             =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
              HRL(I,K)          =  HRL(I,K) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW

              DFNETC            =  TRAN(I,1,K) - TRAN(I,1,KP1) -
     1                             REFL(I,1,K) + REFL(I,1,KP1)
              HRLC(I,K)         =  HRLC(I,K) +
     1                             HRCOEF * DFNETC / DP(I,K) * PGW

          ENDIF
C

          IF (IRAD_FLUX_PROFS .NE. 0) THEN
             FLCU(I,KP1) = FLCU(I,KP1) + REFL(I,1,K)*PGW
             FLCD(I,KP1) = FLCD(I,KP1) + TRAN(I,1,K)*PGW
          END IF

  600       CONTINUE
  650      CONTINUE
C
           DO 660 I = IL1, IL2
            FLXU(I,LEV)         =  FLXU(I,LEV) + REFL(I,2,LEV) * PGW
            FLXD(I,LEV)         =  FLXD(I,LEV) + TRAN(I,2,LEV) * PGW

            IF (IRAD_FLUX_PROFS .NE. 0) THEN
               FLCU(I,LEV+1) = FLCU(I,LEV+1) + REFL(I,1,LEV)*PGW
               FLCD(I,LEV+1) = FLCD(I,LEV+1) + TRAN(I,1,LEV)*PGW
            END IF
C
            OLRC(I)             =  OLRC(I) + (REFL(I,1,LEV1)) * PGW
            OLR(I)              =  OLR(I)  + (REFL(I,2,LEV1)) * PGW
            FLG(I)              =  FLG(I) -
     1                            (REFL(I,2,LEV) - TRAN(I,2,LEV)) * PGW
            FDLC(I)             =  FDLC(I) + TRAN(I,1,LEV) * PGW
            CLB(I)              =  CLB(I) -
     1                            (REFL(I,1,LEV) - TRAN(I,1,LEV)) * PGW
C
  660      CONTINUE
C
             DO 665 M = 1, NTILE
             DO 665 I = IL1, IL2
                IF (ITILE(I,M) .GT. 0) THEN
                   FLXUT(I,M,LEV)      = FLXUT(I,M,LEV)
     1                                 + REFLT(I,M,2,LEV)*PGW
                   FLXDT(I,M,LEV)      = FLXDT(I,M,LEV)
     1                                 + TRANT(I,M,2,LEV)*PGW
                   FLGT (I,M)          = FLGT(I,M) -
     1                               (REFLT(I,M,2,LEV)-TRANT(I,M,2,LEV))
     2                                 * PGW
                   FDLCT(I,M)          = FDLCT(I,M) + TRANT(I,M,1,LEV)
     1                                 * PGW
                   CLBT (I,M)          =  CLBT(I,M) -
     1                               (REFLT(I,M,1,LEV)-TRANT(I,M,1,LEV))
     2                                 * PGW
                END IF
  665        CONTINUE
C
           IF (LEV1 .GT. 1)                                         THEN
            DO 680 K = LEV1 - 1, 1, - 1
              KP1 =  K + 1
              DO 670 I = IL1, IL2
                FLXU(I,K)       =  FLXU(I,K) + REFL(I,2,LEV1) * PGW
                FLXD(I,K)       =  FLXD(I,K) + C1(I) * PGW
C
                IF (IRAD_FLUX_PROFS .NE. 0) THEN
                   FLCU(I,KP1) = FLCU(I,KP1) + REFL(I,1,LEV1)*PGW
                   FLCD(I,KP1) = FLCD(I,KP1) + C1(I)*PGW
                END IF

  670         CONTINUE
  680       CONTINUE
           ENDIF
           IF (IRAD_FLUX_PROFS .NE. 0) THEN
              DO I = IL1, IL2
! Top of atmosphere
                 FLCU(I,1) = FLCU(I,1) + REFL(I,1,LEV1)*PGW
                 FLCD(I,1) = FLCD(I,1) + C1(I)*PGW
                 FLAU(I,1) = FLAU(I,1) + REFL(I,2,LEV1)*PGW
                 FLAD(I,1) = FLAD(I,1) + C1(I)*PGW
              END DO !I
           END IF

          IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     SEPARATE THE COOLING RATE FOR EACH GAS FOR MAM PROJECT           C
C     SINCE THE INTERACTION BETWEEN EACH GASES, THE HEATING RATE COULD C
C     NOT BE EXACTLY SEPARATED. FOR MAIN K INTERVALS, HRSO3 ARE        C
C     OBTAINED FOR EACH K OF VISIBLE BAND                              C
C                                                                      C
C     H2O(1), CO2(1),    2500 - 2200 (CM-1) N2O IS COMBINED WITH CO2   C
C     H2O(2),            2200 - 1900        N2O IS COMBINED WITH H2O(2)C
C     H2O(3),            1900 - 1400                                   C
C     H2O(4), CH4, N2O,  1100 - 1400                                   C
C     H2O(5), O3         980  - 1100                                   C
C     H2O(6), CO2(2),    800  - 980                                    C
C     H2O(7), CO2(2),    540  - 800                                    C
C     H2O(8),            340  - 540                                    C
C     H2O(9),            0    - 340                                    C
C----------------------------------------------------------------------C
C
          IF (IB .NE. 5 .AND. IB .NE. 6)                            THEN
            DO K = LEV1, LAY
              KP1 = K + 1
              DO I = IL1, IL2
                DFNET           =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                HRLH2O(I,K,IB)  =  HRLH2O(I,K,IB) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
              ENDDO
            ENDDO
C
          ELSE IF (IB .EQ. 5)                                       THEN
            DO K = LEV1, LAY
              KP1 = K + 1
              DO I = IL1, IL2
                DFNET           =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                HRLO3(I,K)      =  HRLO3(I,K) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
              ENDDO
            ENDDO
C
          ELSE IF (IB .EQ. 6)                                       THEN
            IF (IG .EQ. 1)                                          THEN
              DO K = LEV1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLCO2(I,K,2) =  HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ELSE
              DO K = LEV1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLH2O(I,K,6) =  HRLH2O(I,K,6) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
          ENDIF
          ENDIF

C
          ELSE IF(GH .AND. IB.NE.6) THEN ! Perform integration over these ig points
C
           CALL GASOPTLGH6(TAUG, GWGH, DP, IB, IG, O3, QG, CO2, CH4,
     1                     AN2O, INPT, MCONT, PG, DIP, DT, LEV1, GH,
     2                     IL1, IL2, ILG, LAY)
C
C----------------------------------------------------------------------C
C     CONSIDER THE ATTENUATION FOR THE DOWNWARD FLUX ABOVE THE MODEL   C
C     TOP LEVEL. THIS IS IMPORTANT TO GET THE CORRECT COOLING RATE. IF C
C     THE MODEL TOP LEVEL PRESSURE IS LOWER THAN 0.01. THIS IS NOT     C
C     NECESSARY                                                        C
C----------------------------------------------------------------------C
C
           CALL LATTENU6(A1, IB, IG, O3TOP, QG, CO2, PFULL, A1(1,12),
     1                   DT, A1(1,5), INPT, IL1, IL2, ILG)
C
           DO 710 I = IL1, IL2
            TRAN0             =  EXP(- A1(I,1))
C
            IF (PFULL(I,1) .GT. 0.001)                            THEN
              X               =  MAX(A1(I,1), 1.E-10)
              UBETA0          =  1.6487213 * A1(I,3) / X
              EPSD0           =  UBETA0 + 1.0
              IF (ABS(EPSD0) .GT. 0.001)                          THEN
                C2(I)         =  C1(I) * TRAN0 + (BF(I,1) - A1(I,2) *
     1                           TRAN0) / EPSD0
              ELSE
                C2(I)         =  C1(I) * TRAN0 + X * A1(I,2) * TRAN0
              ENDIF
            ELSE
              C2(I)           =  C1(I) * TRAN0
            ENDIF
  710      CONTINUE
C
           CALL LWTRAGH4(REFLT, TRANT,
     1                   C2, TAUCI, OMCI, TAUA, TAUG,
     2                   BF, URBF, CLDG, EM0TL, BST, ITILE, CUT,
     3                   IL1, IL2, ILG, LAY, LEV,
     4                   NTILE)
C
C----------------------------------------------------------------------C
C Compute the gridbox mean REFL and TRAN using REFLT, TRANT and the    C
C fraction covered by each tile.                                       C
C----------------------------------------------------------------------C
C
          REFL(IL1:IL2,1:2,1:LEV) = 0.0
          TRAN(IL1:IL2,1:2,1:LEV) = 0.0

          DO K = 1, LEV
             DO M = 1, NTILE
                DO I = IL1, IL2
                   IF (ITILE(I,M) .GT. 0) THEN
                      TRAN(I,1,K) = TRAN(I,1,K)
     1                            + TRANT(I,M,1,K)*FAREL(I,M)
                      TRAN(I,2,K) = TRAN(I,2,K)
     1                            + TRANT(I,M,2,K)*FAREL(I,M)
                      REFL(I,1,K) = REFL(I,1,K)
     1                            + REFLT(I,M,1,K)*FAREL(I,M)
                      REFL(I,2,K) = REFL(I,2,K)
     1                            + REFLT(I,M,2,K)*FAREL(I,M)
                   END IF
                END DO ! I
             END DO ! M
          END DO ! K

C
C          * FOR THE TOTAL SKY FLUXES WEIGHT THE CLEAR AND CLOUDY SKY
C          * FLUXES BY THE TOTAL VERTICALLY PROJECTED CLOUD FRACTION.
C
           DO K = 1, LEV
           DO I = IL1,IL2
             IF (CLDT(I).LT.1.0) THEN
                 REFL(I,2,K) = (1.0 - CLDT(I))*REFL(I,1,K)
     1                          + CLDT(I)  * REFL(I,2,K)
                 TRAN(I,2,K) = (1.0 - CLDT(I))*TRAN(I,1,K)
     1                          + CLDT(I)  * TRAN(I,2,K)
             END IF
           END DO ! I
           END DO ! K
C
           DO K = 1, LEV
           DO M = 1, NTILE
           DO I = IL1, IL2
              IF (ITILE(I,M) .GT. 0) THEN
                 IF (CLDT(I).LT.1.0) THEN
                    REFLT(I,M,2,K) = (1.0 - CLDT(I)) * REFLT(I,M,1,K) +
     1                               CLDT(I)  * REFLT(I,M,2,K)
                    TRANT(I,M,2,K) = (1.0 - CLDT(I)) * TRANT(I,M,1,K) +
     1                                CLDT(I)  * TRANT(I,M,2,K)
                 END IF
              END IF
           END DO ! I
           END DO ! M
           END DO ! K

           PGW = PI * GWGH * WWW
C
           DO 740 K = 1, LAY
            KP1 = K + 1
            DO 730 I = IL1, IL2
              FLXU(I,K)       =  FLXU(I,K) + REFL(I,2,K) * PGW
              FLXD(I,K)       =  FLXD(I,K) + TRAN(I,2,K) * PGW
              IF (IRAD_FLUX_PROFS .NE. 0) THEN
                 FLCU(I,KP1) = FLCU(I,KP1) + REFL(I,1,K)*PGW
                 FLCD(I,KP1) = FLCD(I,KP1) + TRAN(I,1,K)*PGW
              END IF
C
C         * FOR NON-MAM APPLICATIONS, CALCULATE HEATING RATE IN
C         * USUAL MANNER.
C
          IF(IMAM.EQ.0)                      THEN
              DFNET           =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                           REFL(I,2,K) + REFL(I,2,KP1)
              HRL(I,K)        =  HRL(I,K) +
     1                           HRCOEF * DFNET / DP(I,K) * PGW

              DFNETC          =  TRAN(I,1,K) - TRAN(I,1,KP1) -
     1                           REFL(I,1,K) + REFL(I,1,KP1)
              HRLC(I,K)       =  HRLC(I,K) +
     1                           HRCOEF * DFNETC / DP(I,K) * PGW
          ENDIF
C

  730       CONTINUE
  740      CONTINUE
C
C----------------------------------------------------------------------C
C Compute the TOA OLR accounting for the presence of a "moon layer"    C
C with some assumed properties.  This must be consistent with what was C
C used for the downward flux between the TOA and top of model.         C
C----------------------------------------------------------------------C
C

           DO 720 I = IL1, IL2
              TRAN0             =  EXP(- A1(I,1))
C
              IF (PFULL(I,1) .GT. 0.001)                            THEN
                 X               =  MAX(A1(I,1), 1.E-10)
                 UBETA0          =  1.6487213 * A1(I,3) / X
                 EPSD0           =  UBETA0 - 1.0
                 IF (ABS(EPSD0) .GT. 0.001)                         THEN
                    TOA_OLR(I,1)  =  REFL(I,1,1) * TRAN0
     1                            + (BF(I,1)*TRAN0 - A1(I,2))/EPSD0
                    TOA_OLR(I,2)  =  REFL(I,2,1) * TRAN0
     1                            + (BF(I,1)*TRAN0 - A1(I,2))/EPSD0
                 ELSE
                    TOA_OLR(I,1)  =  REFL(I,1,1) * TRAN0
     1                            +  X * BF(I,1) * TRAN0
                    TOA_OLR(I,2)  =  REFL(I,2,1) * TRAN0
     1                            +  X * BF(I,1) * TRAN0
                 ENDIF
              ELSE
                 TOA_OLR(I,1)    =  REFL(I,1,1) * TRAN0
                 TOA_OLR(I,2)    =  REFL(I,2,1) * TRAN0
              ENDIF

              IF (CLDT(I).LT.1.0) THEN
                 TOA_OLR(I,2) = (1.0-CLDT(I))*TOA_OLR(I,1)
     1                        + CLDT(I)*TOA_OLR(I,2)
              END IF
 720       CONTINUE

           DO 750 I = IL1, IL2
            FLXU(I,LEV)       =  FLXU(I,LEV) + REFL(I,2,LEV) * PGW
            FLXD(I,LEV)       =  FLXD(I,LEV) + TRAN(I,2,LEV) * PGW
            IF (IRAD_FLUX_PROFS .NE. 0) THEN
               FLCU(I,LEV+1) = FLCU(I,LEV+1) + REFL(I,1,LEV)*PGW
               FLCD(I,LEV+1) = FLCD(I,LEV+1) + TRAN(I,1,LEV)*PGW
               FLCU(I,1) = FLCU(I,1) + TOA_OLR(I,1)*PGW
               FLAU(I,1) = FLAU(I,1) + TOA_OLR(I,2)*PGW
            END IF

            OLRC(I)           =  OLRC(I) + TOA_OLR(I,1) * PGW
            OLR(I)            =  OLR(I)  + TOA_OLR(I,2) * PGW
            FLG(I)            =  FLG(I) -
     1                           (REFL(I,2,LEV) - TRAN(I,2,LEV)) * PGW
            FDLC(I)           =  FDLC(I) + TRAN(I,1,LEV) * PGW
            FLAMOON(I)        =  FLAMOON(I) + ( (C1(I)-TOA_OLR(I,2)) -
     2                                          (C2(I)-REFL(I,2,1)))*PGW

            CLB(I)            =  CLB(I) -
     1                           (REFL(I,1,LEV) - TRAN(I,1,LEV)) * PGW
C
  750      CONTINUE
             DO 760 M = 1, NTILE
             DO 760 I = IL1, IL2
                IF (ITILE(I,M) .GT. 0) THEN
                   FLXUT(I,M,LEV) = FLXUT(I,M,LEV)
     1                            + REFLT(I,M,2,LEV) * PGW
                   FLXDT(I,M,LEV) = FLXDT(I,M,LEV)
     1                            + TRANT(I,M,2,LEV) * PGW
                   FLGT(I,M)      = FLGT(I,M)
     1                            -(REFLT(I,M,2,LEV)-TRANT(I,M,2,LEV))
     2                            * PGW
                   FDLCT(I,M)     = FDLCT(I,M) + TRANT(I,M,1,LEV)* PGW
                   CLBT(I,M)      = CLBT(I,M)
     1                            -(REFLT(I,M,1,LEV)-TRANT(I,M,1,LEV))
     2                            * PGW
                END IF
 760         CONTINUE
          IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     CALCULATION OF COOLING RATES FOR SEPARATE GASES,                 C
C----------------------------------------------------------------------C
C
            IF (IB.EQ.2 .OR. IB.EQ.3 .OR.
     1          (IB.EQ.4 .AND. IG.EQ.1) .OR.
     2          (IB.EQ.5 .AND. IG.EQ.2) .OR. IB.GE.8)      THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLH2O(I,K,IB)=  HRLH2O(I,K,IB) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
C
            IF (IB .EQ. 5 .AND. IG .NE. 2)                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLO3(I,K)    =  HRLO3(I,K) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
C
            IF (IB .EQ. 1)                                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLCO2(I,K,1) =  HRLCO2(I,K,1) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
C
            IF (IB .EQ. 7)                                          THEN
C
C----------------------------------------------------------------------C
C     15 UM BAND CO2 AND H2O ARE STRONGLY OVERLAPPED. IT IS DIFFICULT  C
C     TO SEPARATE THE COOLING OF THEM IN THE LOWER ATMOSPHERE, WE PUT  C
C     MOST COOLING RATE IN H2O(7) BELOW 100 MB (FOR ALL MAIN K VALUES),C
C     AND ASSUME ALL COOLING IN THE UPPER PART ARE FROM CO2 (FOR ALL   C
C     MINOR K VALUES).                                                 C
C----------------------------------------------------------------------C
C
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLCO2(I,K,2) =  HRLCO2(I,K,2) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
C
            IF (IB .EQ. 4 .AND. IG .GE. 3)                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLCH4(I,K)   =  HRLCH4(I,K) +
     1                             HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
            IF (IB .EQ. 4 .AND. IG .EQ. 2)                          THEN
              DO K = 1, LAY
                KP1 = K + 1
                DO I = IL1, IL2
                  DFNET         =  TRAN(I,2,K) - TRAN(I,2,KP1) -
     1                             REFL(I,2,K) + REFL(I,2,KP1)
                  HRLN2O(I,K)   =  HRCOEF * DFNET / DP(I,K) * PGW
                ENDDO
              ENDDO
            ENDIF
          ENDIF

          ENDIF
  850    CONTINUE
  875   CONTINUE
  900 CONTINUE
C
      DO 950 I = IL1, IL2
        FDL(I)                  =  FLXD(I,LEV)
  950 CONTINUE
C
        DO 975 M = 1, NTILE
        DO 975 I = IL1, IL2
           IF (ITILE(I,M) .GT. 0) THEN
              FDLT(I,M) = FLXDT(I,M,LEV)
           END IF
  975   CONTINUE

        IF(IMAM.NE.0)                      THEN
C
C----------------------------------------------------------------------C
C     ADD THE COOLING RATE FOR EACH GAS TOGETHER                       C
C----------------------------------------------------------------------C
C
           DO K = 1, LAY
              DO I = IL1, IL2
                 HRL(I,K) =  HRLH2O(I,K,1) + HRLH2O(I,K,2)
     1                    + HRLH2O(I,K,3) + HRLH2O(I,K,4)
     2                    + HRLH2O(I,K,5) + HRLH2O(I,K,6)
     3                    + HRLH2O(I,K,7) + HRLH2O(I,K,8)
     4                    + HRLH2O(I,K,9)
     5                    + HRLO3(I,K)
     6                    + HRLCO2(I,K,1) + HRLCO2(I,K,2)
     7                    + HRLCH4(I,K) + HRLN2O(I,K)
              ENDDO
           ENDDO
        ENDIF ! IMAM
C
      IF (IRAD_FLUX_PROFS .NE. 0) THEN
         DO K = 1, LEV
            KP1 = K + 1
            DO I = IL1, IL2
               FLAU(I,KP1) = FLXU(I,K)
               FLAD(I,KP1) = FLXD(I,K)
            END DO ! I
         END DO ! K
      END IF
C
      ENDIF
C     (LCLW)
C
      RETURN
      END
