C     *****************************************************************
C          SUBROUTINE FST88 IS THE MAIN COMPUTATION MODULE OF THE
C     LONG-WAVE RADIATION CODE. IN IT ALL "EMISSIVITY" CALCULATIONS,
C     INCLUDING CALLS TO TABLE LOOKUP SUBROUTINES. ALSO,AFTER CALLING
C     SUBROUTINE "SPA88", FINAL COMBINED HEATING RATES AND GROUND
C     FLUX ARE OBTAINED.
C     *****************************************************************
C              INPUTS:
C        BETINW,BETAWD,AB15WD              BDWIDE
C        BETAD,BO3RND,AO3RND               BANDTA
C        CLDFAC                            CLDCOM
C        QH2O,P,DELP2,DELP,T,VAR1,VAR2,    KDACOM
C        VAR3,VAR4,CNTVAL                  KDACOM
C        TOTVO2,TOTO3,TOTPHI,EMPL,EMX1     KDACOM
C        TPHIO3,EMX2                       KDACOM
C        TEMP,PRESS                        RADISW
C        NCLDS,KTOP,KBTM,CAMT              RADISW
C        IND,INDX2,KMAXV,SOURCE,DSRCE      TABCOM
C        SKC1R,SKC3R,KMAXVM,NREP1,NREP2    TABCOM
C        NST1,NST2,NRP1,NRP2               TABCOM
C        CO2NBL,CO21                       TFCOM
C        CO2SP1,CO2SP2                     TFCOM
C              OUTPUTS:
C        HEATRA,GRNFLX,TOPFLX              LWOUT
C
C          CALLED BY  :    RADMN OR MAIN PGM
C          CALLS      :    CLO88,E1E288,E3V88,SPA88,NLTE
C
C        PASSED VARIABLES:
C              IN E3V88:
C        EMD     =  E3 FUNCTION FOR H2O LINES (0-560,1200-2200 CM-1)
C                     COMPUTED IN E3V88
C        TPL     =  TEMPERATURE INPUT FOR E3 CALCULATION IN E3V88
C        EMPL    =  H2O AMOUNT,INPUT FOR E3 CALCULATION IN E3V88
C                   (COMPUTED IN LWR88; STORED IN KDACOM.H)
C              IN E1E288:
C        E1CTS1  =  E1 FUNCTION FOR THE (I+1)TH LEVEL USING THE
C                   TEMPERATURE OF THE ITH DATA LEVEL,COMPUTED OVER
C                   THE FREQUENCY RANGE 0-560,1200-2200 CM-1. (E1CTS1-
C                   E1CTW1) IS USED IN OBTAINING THE FLUX AT THE TOP
C                   IN THE 0-160,1200-2200 CM-1 RANGE (FLX1E1).
C        E1CTS2  =  E1 FUNCTION FOR THE ITH LEVEL, USING THE TEMP. OF
C                   THE ITH DATA LEVEL,COMPUTED OVER THE FREQUENCY RANGE
C                   0-560,1200-2200 CM-1. (E1CTS2-E1CTW2) IS ALSO USED
C                   IN OBTAINING THE FLUX AT THE TOP IN THE 0-160,.
C                   1200-2200 CM-1 RANGE.
C        E1FLX   =  E1 FCTN. FOR THE ITH LEVEL,USING THE TEMPERATURE AT
C                   THE TOP OF THE ATMOSPHERE. COMPUTED OVER THE FREQ.
C                   RANGE 0-560,1200-2200 CM-1. USED FOR Q(APPROX) TERM.
C                   (IN COMMON BLOCK TFCOM)
C        E1CTW1  =  LIKE E1CTS1,BUT COMPUTED OVER THE 160-560 CM-1 RANGE
C                   AND USED FOR Q(APPROX,CTS) CALCULATION
C        E1CTW2  =  LIKE E1CTS2,BUT COMPUTED OVER THE 160-560 CM-1 RANGE
C                   AND USED FOR Q(APPROX,CTS) CALCULATION
C        FXO     =  TEMPERATURE INDEX USED FOR E1 FUNCTION AND ALSO
C                   USED FOR SOURCE FUNCTION CALC. IN FST88.
C        DT      =  TEMP. DIFF.BETWEEN MODEL TEMPS. AND TEMPS. AT
C                   TABULAR VALUES OF E1 AND SOURCE FCTNS. USED IN
C                   FST88 AND IN E1 FUNCTION CALC.
C        FXOE2   =  TEMPERATURE INDEX USED FOR E2 FUNCTION
C        DTE2    =  TEMP. DIFF. BETWEEN MODEL TEMP. AND TEMPS. AT
C                   TABULAR VALUES OF E2 FUNCTION.
      SUBROUTINE FST88(HEATRA,GRNFLX,TOPFLX,
     1                 QH2O,PRESS,P,DELP,DELP2,TEMP,T,
     2                 CLDFAC,NCLDS,KTOP,KBTM,CAMT,
     3                 CO21,CO2NBL,CO2SP1,CO2SP2,
     4                 VAR1,VAR2,VAR3,VAR4,CNTVAL,
     5                 TOTO3,TPHIO3,TOTPHI,TOTVO2,
     6                 EMX1,EMX2,EMPL)
C
      COMMON/PHYCON/AMOLWT,CSUBP,DIFFCTR,G,GRAVDR,O3DIFCTR,P0,
     *            P0XZP2,P0XZP8,P0X2,RADCON,RGAS,RGASSP,SECPDA
      COMMON/PHYCON/RATCO2MW,RATH2OMW
      COMMON/PHYCON/RADCON1
      COMMON/PHYCON/GINV,P0INV,GP0INV
      COMMON/HCON/HUNDRED,HNINETY,SIXTY,FIFTY,TEN,EIGHT,FIVE,
     *            FOUR,THREE,TWO,ONE,HAF,QUARTR,ZERO
      COMMON/HCON/H83E26,H71E26,H1E15,H1E13,H1E11,H1E8,H4E5,
     *            H165E5,H5725E4,H488E4,H1E4,H24E3,H20788E3,
     *            H2075E3,H1224E3,H5E2,H3082E2,H3E2,H2945E2,
     *            H23E2,H15E2,H35E1,H3P6,H181E1,H18E1,H2P9,H2P8,
     *            H2P5,H1P8,H1P4387,H1P4,H1P25892,HP8,HP518,
     *            HP369,HP1
      COMMON/HCON/H44871M2,H559M3,H1M3,H987M4,H285M4,H1M4,
     *            H6938M5,H394M5,H37412M5,H1439M5,H128M5,H1M5,
     *            H7M6,H4999M6,H25452M6,H1M6,H391M7,H1174M7,
     *            H8725M8,H327M8,H257M8,H1M8,H23M10,H14M10,
     *            H11M10,H1M10,H83M11,H82M11,H8M11,H77M11,
     *            H72M11,H53M11,H48M11,H44M11,H42M11,H37M11,
     *            H35M11,H32M11,H3M11,H28M11,H24M11,H23M11,
     *            H2M11,H18M11,H15M11,H14M11,H114M11,H11M11,
     *            H1M11,H96M12,H93M12,H77M12,H74M12,H65M12,
     *            H62M12,H6M12,H45M12,H44M12,H4M12,H38M12,
     *            H37M12,H3M12,H29M12,H28M12,H24M12,H21M12,
     *            H16M12,H14M12,H12M12,H8M13,H46M13,H36M13,
     *            H135M13,H12M13,H1M13,H3M14,H15M14,H14M14,
     *            H1M17,H1M18,H1M19,H1M20,H1M21,H1M22,H1M23,
     *            H1M24,H26M30,H14M30,H25M31,H21M31,H12M31,
     *            H9M32,H55M32,H45M32,H4M33,H62M34,H1M60
      COMMON/HCON/HMP575,HM13EZ,HM19EZ,HM1E1,HM181E1,HM1E2
      COMMON/HCON/H1E6,H2E6,H1M2,HMP66667,HM6666M2,HP166666,
     *            H41666M2,HMP5,HM2M2,H29316E2,H1226E1,H3116E1,
     *            H9P94,HP6,H625M2,HP228,HP60241,HM1797E1,
     *            H8121E1,H2E2,HM1EZ,H26E2,H44194M2,H1P41819
      COMMON/HCON/HP219,HP144,HP816,H69766E5,H235M3,HP26,
     *            H129M2,H75826M4,H1P082,HP805,H1386E2,
     *            H658M2,H1036E2,H2118M2,H42M2,H323M4,
     *            H67390E2,HP3795,HP5048,H102M5,H451M6
      COMMON/HCON/H16E1,HM161E1,H161E1,H3M3,H101M16,
     *            HM1597E1,H25E2,HP118666,H15M5,H3P5,H18E3,
     *            H6P08108,HMP805,HP602409,HP526315,
     *            H28571M2,H1M16
      COMMON/HCON/H3M4
      COMMON/HCON/HM8E1
      COMMON/HCON/H28E1
C-----------------------------------------------------------------------
       INCLUDE "parmeta"
       INCLUDE "mpp.h"  
#include "sp.h"
C-----------------------------------------------------------------------
C     PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE:
C          IMAX   =  NO. POINTS ALONG THE LAT. CIRCLE USED IN CALCS.
C          L      =  NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL
C***NOTE: THE USER NORMALLY WILL MODIFY ONLY THE IMAX AND L PARAMETERS
C          NBLW   =  NO. FREQ. BANDS FOR APPROX COMPUTATIONS. SEE
C                      BANDTA FOR DEFINITION
C          NBLX   =  NO. FREQ BANDS FOR APPROX CTS COMPUTATIONS
C          NBLY   =  NO. FREQ. BANDS FOR EXACT CTS COMPUTATIONS. SEE
C                      BDCOMB FOR DEFINITION
C          INLTE  =  NO. LEVELS USED FOR NLTE CALCS.
C          NNLTE  =  INDEX NO. OF FREQ. BAND IN NLTE CALCS.
C          NB,KO2 ARE SHORTWAVE PARAMETERS; OTHER QUANTITIES ARE DERIVED
C                    FROM THE ABOVE PARAMETERS.
      PARAMETER (L=LM)
      PARAMETER (IMAX=IM,NCOL=IMAX)
      PARAMETER (NBLW=163,NBLX=47,NBLY=15)
      PARAMETER (NBLM=NBLY-1)
      PARAMETER (LP1=L+1,LP2=L+2,LP3=L+3)
      PARAMETER (LM1=L-1,LM2=L-2,LM3=L-3)
      PARAMETER (LL=2*L,LLP1=LL+1,LLP2=LL+2,LLP3=LL+3)
      PARAMETER (LLM1=LL-1,LLM2=LL-2,LLM3=LL-3)
      PARAMETER (LP1M=LP1*LP1,LP1M1=LP1M-1)
      PARAMETER (LP1V=LP1*(1+2*L/2))
      PARAMETER (LP121=LP1*NBLY)
      PARAMETER (LL3P=3*L+2)
      PARAMETER (NB=12)
      PARAMETER (INLTE=3,INLTEP=INLTE+1,NNLTE=56)
      PARAMETER (LP1I=IMAX*LP1,LLP1I=IMAX*LLP1,LL3PI=IMAX*LL3P)
      PARAMETER (NB1=NB-1)
      PARAMETER (KO2=12)
      PARAMETER (KO21=KO2+1,KO2M=KO2-1)
C     PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE:
C          IMAX   =  NO. POINTS SENT TO RADFS
C          L      =  NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL
C***NOTE: THE USER NORMALLY WILL MODIFY ONLY THE IMAX AND L PARAMETERS
C          NBLW   =  NO. FREQ. BANDS FOR APPROX COMPUTATIONS. SEE
C                      BANDTA FOR DEFINITION
C          NBLX   =  NO. FREQ BANDS FOR APPROX CTS COMPUTATIONS
C          NBLY   =  NO. FREQ. BANDS FOR EXACT CTS COMPUTATIONS. SEE
C                      BDCOMB FOR DEFINITION
C          INLTE  =  NO. LEVELS USED FOR NLTE CALCS.
C          NNLTE  =  INDEX NO. OF FREQ. BAND IN NLTE CALCS.
C          NB,KO2 ARE SHORTWAVE PARAMETERS; OTHER QUANTITIES ARE DERIVED
C                    FROM THE ABOVE PARAMETERS.
C    COMMON BLOCK BANDTA CONTAINS RANDOM BAND PARAMETERS FOR THE LW
C    CALCULATIONS USING 10 CM-1 WIDE BANDS.THE 15 UM CO2 COMPLEX
C    IS 2 BANDS,560-670 AND 670-800 CM-1. OZONE COEFFICIENTS ARE
C    IN 3 BANDS,670-800 (14.1 UM),990-1070 AND 1070-1200 (9.6 UM).
C    THE  (NBLW) BANDS NOW INCLUDE:
C                56 BANDS, 10  CM-1 WIDE    0  -   560  CM-1
C                 2 BANDS, 15 UM COMPLEX  560  -   670  CM-1
C                                         670  -   800  CM-1
C                 3 "CONTINUUM" BANDS     800  -   900  CM-1
C                                         900  -   990  CM-1
C                                        1070  -   1200 CM-1
C                 1 BAND FOR 9.6 UM BAND  990  -   1070 CM-1
C               100 BANDS, 10 CM-1 WIDE  1200  -   2200 CM-1
C                 1 BAND FOR 4.3 UM SRC  2270  -   2380 CM-1
C    THUS NBLW PRESENTLY EQUALS    163
C    ALL BANDS ARE ARRANGED IN ORDER OF INCREASING WAVENUMBER
C
C        ARNDM   =   RANDOM "A" PARAMETER FOR (NBLW) BANDS
C        BRNDM   =   RANDOM "B" PARAMETER FOR (NBLW) BANDS
C        BETAD   =   CONTINUUM COEFFICIENTS FOR (NBLW) BANDS
C        AP,BP   =   CAPPHI COEFFICIENTS FOR (NBLW) BANDS
C        ATP,BTP =   CAPPSI COEFFICIENTS FOR (NBLW) BANDS
C        BANDLO  =   LOWEST FREQUENCY IN EACH OF (NBLW) FREQ. BANDS
C        BANDHI  =   HIGHEST FREQUENCY IN EACH OF (NBLW) FREQ. BANDS
C        AO3RND  =   RANDOM "A" PARAMETER FOR OZONE IN (3) OZONE
C                    BANDS
C        BO3RND  =   RANDOM "B" PARAMETER FOR OZONE IN (3) OZONE
C                    BANDS
C        AB15    =   THE PRODUCT ARNDM*BRNDM FOR THE TWO BANDS
C                    REPRESENTING THE 15 UM BAND COMPLEX OF CO2
C     DATA FOR ARNDM,BRNDM,AP,BP,ATP,BTP,AO3RND,BO3RND ARE OBTAINED BY
C     USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS ARE FROM
C     ROBERTS (1976).
      COMMON / BANDTA / ARNDM(NBLW),BRNDM(NBLW),BETAD(NBLW),AP(NBLW),
     1                  BP(NBLW),ATP(NBLW),BTP(NBLW),BANDLO(NBLW),
     2                  BANDHI(NBLW),AO3RND(3),BO3RND(3),AB15(2)
C
C    COMMON BLOCK BDWIDE CONTAINS RANDOM BAND PARAMETERS FOR SPECIFIC
C    WIDE BANDS. AT PRESENT,THE INFORMATION CONSISTS OF 1) RANDOM
C    MODEL PARAMETERS FOR THE 15 UM BAND,560-800 CM-1; 2) THE
C    CONTINUUM COEFFICIENT FOR THE 800-990,1070-1200 CM-1 BAND
C        SPECIFICALLY:
C        AWIDE       =   RANDOM "A" PARAMETER FOR  BAND
C        BWIDE       =   RANDOM "B" PARAMETER FOR  BAND
C        BETAWD      =   CONTINUUM COEFFICIENTS FOR BAND
C        APWD,BPWD   =   CAPPHI COEFFICIENTS FOR  BAND
C        ATPWD,BTPWD =   CAPPSI COEFFICIENTS FOR BAND
C        BDLOWD      =   LOWEST FREQUENCY IN EACH  FREQ  BAND
C        BDHIWD      =   HIGHEST FREQUENCY IN EACH FREQ  BAND
C        AB15WD      =   THE PRODUCT ARNDM*BRNDM FOR THE ONE BAND
C                        REPRESENTING THE 15 UM BAND COMPLEX OF CO2
C        BETINW      =   CONT.COEFFICIENT FOR A SPECIFIED WIDE
C                        FREQ.BAND (800-990 AND 1070-1200 CM-1).
C        SKO2D       =   1./BETINW, USED IN SPA88 FOR CONT. COEFFS
C        SKC1R       =   BETAWD/BETINW, USED FOR CONT. COEFF. FOR
C                        15 UM BAND IN FST88
C        SKO3R       =   RATIO OF CONT. COEFF. FOR 9.9 UM BAND TO
C                        BETINW, USED FOR 9.6 UM CONT COEFF IN FST88
C     DATA FOR AWIDE,BWIDE,APWD,BPWD,ATPWD,BTPWD,AO3WD,BO3WD ARE
C     OBTAINED BY USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS
C     ARE FROM ROBERTS (1976).
      COMMON / BDWIDE / AWIDE,BWIDE,BETAWD,
     1                  APWD,BPWD,ATPWD,BTPWD,
     2                  BDLOWD,BDHIWD,BETINW,
     3                  AB15WD,SKO2D,SKC1R,SKO3R
C
C    COMMON BLOCK BDCOMB CONTAINS RANDOM BAND PARAMETERS FOR THE LW
C    CALCULATIONS USING COMBINED WIDE FREQUENCY BANDS BETWEEN 160 AND
C    1200 CM-1,AS WELL AS THE 2270-2380 BAND FOR SOURCE CALC.
C        BANDS 1-8: COMBINED WIDE FREQUENCY BANDS FOR 160-560 CM-1
C        BANDS 9-14: FREQUENCY BANDS,AS IN BANDTA (NARROW BANDS)
C                    FOR 560-1200 CM-1
C        BAND  15:  FREQUENCY BAND 2270-2380 CM-1,USED FOR SOURCE
C                   CALCULATION ONLY
C        THUS NBLY PRESENTLY EQUALS   15
C
C        BANDS ARE ARRANGED IN ORDER OF INCREASING WAVENUMBER
C        ACOMB       =   RANDOM "A" PARAMETER FOR (NBLY) BANDS
C        BCOMB       =   RANDOM "B" PARAMETER FOR (NBLY) BANDS
C        BETACM      =   CONTINUUM COEFFICIENTS FOR (NBLY) BANDS
C        APCM,BPCM   =   CAPPHI COEFFICIENTS FOR (NBLY) BANDS
C        ATPCM,BTPCM =   CAPPSI COEFFICIENTS FOR (NBLY) BANDS
C        BDLOCM      =   LOWEST FREQUENCY IN EACH OF (NBLY) FREQ. BANDS
C        BDHICM      =   HIGHEST FREQUENCY IN EACH OF (NBLY) FREQ. BANDS
C        AO3CM       =   RANDOM "A" PARAMETER FOR OZONE IN (3) OZONE
C                        BANDS
C        BO3CM       =   RANDOM "B" PARAMETER FOR OZONE IN (3) OZONE
C                        BANDS
C        AB15CM      =   THE PRODUCT ARNDM*BRNDM FOR THE TWO BANDS
C                        REPRESENTING THE 15 UM BAND COMPLEX OF CO2
C        BETINC      =   CONT.COEFFICIENT FOR A SPECIFIED WIDE
C                        FREQ.BAND (800-990 AND 1070-1200 CM-1).
C        IBAND       =   INDEX NO OF THE 40 WIDE BANDS USED IN
C                        COMBINED WIDE BAND CALCULATIONS. IN OTHER
C                        WORDS,INDEX TELLING WHICH OF THE 40 WIDE
C                        BANDS BETWEEN 160-560 CM-1 ARE INCLUDED IN
C                        EACH OF THE FIRST 8 COMBINED WIDE BANDS
C     DATA FOR ACOMB,BCOMB,APCM,BPCM,ATPCM,BTPCM,AO3CM,BO3CM ARE
C     OBTAINED BY USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS
C     ARE FROM ROBERTS (1976). IBAND INDEX VALUES ARE OBTAINED BY
C     EXPERIMENTATION.
      COMMON / BDCOMB / IBAND(40),ACOMB(NBLY),BCOMB(NBLY),
     1                  BETACM(NBLY),APCM(NBLY),BPCM(NBLY),ATPCM(NBLY),
     2                  BTPCM(NBLY),BDLOCM(NBLY),BDHICM(NBLY),BETINC,
     3                  AO3CM(3),BO3CM(3),AB15CM(2)
C
C     COMMON BLOCK TABCOM CONTAINS QUANTITIES PRECOMPUTED IN SUBROUTINE
C     TABLE FOR USE IN THE LONGWAVE RADIATION PROGRAM:
C          EM1     =  E1 FUNCTION, EVALUATED OVER THE 0-560 AND
C                     1200-2200 CM-1 INTERVALS
C          EM1WDE  =  E1 FUNCTION, EVALUATED OVER THE 160-560 CM-1
C                     INTERVAL
C          TABLE1  =  E2 FUNCTION, EVALUATED OVER THE 0-560 AND
C                     1200-2200 CM-1 INTERVALS
C          TABLE2  =  TEMPERATURE DERIVATIVE OF TABLE1
C          TABLE3  =  MASS DERIVATIVE OF TABLE1
C          EM3     =  E3 FUNCTION, EVALUATED OVER THE 0-560 AND
C                     1200-2200 CM-1 INTERVALS
C          SOURCE  =  PLANCK FUNCTION, EVALUATED AT SPECIFIED TEMPS. FOR
C                     BANDS USED IN CTS CALCULATIONS
C          DSRCE   =  TEMPERATURE DERIVATIVE OF SOURCE
C          IND     =  INDEX, WITH VALUE IND(I)=I. USED IN FST88
C          INDX2   =  INDEX VALUES USED IN OBTAINING "LOWER TRIANGLE"
C                     ELEMENTS OF AVEPHI,ETC.,IN FST88
C          KMAXV   =  INDEX VALUES USED IN OBTAINING "UPPER TRIANGLE"
C                     ELEMENTS OF AVEPHI,ETC.,IN FST88
C          KMAXVM  =  KMAXV(L),USED FOR DO LOOP INDICES
C
      COMMON / TABCOM / IND(IMAX),INDX2(LP1V),KMAXV(LP1),
     1     KMAXVM
      COMMON/TABCOM/EM1(28,180),EM1WDE(28,180),TABLE1(28,180),
     1 TABLE2(28,180),TABLE3(28,180),EM3(28,180),SOURCE(28,NBLY),
     2 DSRCE(28,NBLY)
C
      DIMENSION QH2O(IDIM1:IDIM2,LP1),PRESS(IDIM1:IDIM2,LP1)
      DIMENSION P(IDIM1:IDIM2,LP1),DELP(IDIM1:IDIM2,L),
     & DELP2(IDIM1:IDIM2,L),TEMP(IDIM1:IDIM2,LP1)
      DIMENSION T(IDIM1:IDIM2,LP1),CLDFAC(IDIM1:IDIM2,LP1,LP1),
     & CAMT(IDIM1:IDIM2,LP1)
      DIMENSION NCLDS(IDIM1:IDIM2),KTOP(IDIM1:IDIM2,LP1),
     & KBTM(IDIM1:IDIM2,LP1)
      DIMENSION CO21(IDIM1:IDIM2,LP1,LP1),CO2NBL(IDIM1:IDIM2,L)
      DIMENSION CO2SP1(IDIM1:IDIM2,LP1),CO2SP2(IDIM1:IDIM2,LP1)
      DIMENSION VAR1(IDIM1:IDIM2,L),VAR2(IDIM1:IDIM2,L),
     & VAR3(IDIM1:IDIM2,L),VAR4(IDIM1:IDIM2,L)
      DIMENSION CNTVAL(IDIM1:IDIM2,LP1)
      DIMENSION HEATRA(IDIM1:IDIM2,L),GRNFLX(IDIM1:IDIM2),
     & TOPFLX(IDIM1:IDIM2)
      DIMENSION GXCTS(IDIM1:IDIM2),FLX1E1(IDIM1:IDIM2)
      DIMENSION AVEPHI(IDIM1:IDIM2,LP1),EMISS(IDIM1:IDIM2,LP1),
     & EMISSB(IDIM1:IDIM2,LP1)
C
      DIMENSION TOTO3(IDIM1:IDIM2,LP1),TPHIO3(IDIM1:IDIM2,LP1),
     & TOTPHI(IDIM1:IDIM2,LP1)
      DIMENSION TOTVO2(IDIM1:IDIM2,LP1),EMX1(IDIM1:IDIM2),
     & EMX2(IDIM1:IDIM2),EMPL(IDIM1:IDIM2,LLP1)
C
      DIMENSION EXCTS(IDIM1:IDIM2,L),CTSO3(IDIM1:IDIM2,L),
     & CTS(IDIM1:IDIM2,L),E1FLX(IDIM1:IDIM2,LP1)
      DIMENSION CO2SP(IDIM1:IDIM2,LP1),TO3SPC(IDIM1:IDIM2,L),
     & TO3SP(IDIM1:IDIM2,LP1)
      DIMENSION OSS(IDIM1:IDIM2,LP1),CSS(IDIM1:IDIM2,LP1),
     & SS1(IDIM1:IDIM2,LP1),SS2(IDIM1:IDIM2,LP1),
     1 TC(IDIM1:IDIM2,LP1),DTC(IDIM1:IDIM2,LP1)
      DIMENSION SORC(IDIM1:IDIM2,LP1,NBLY),CSOUR(IDIM1:IDIM2,LP1)
CCC
      DIMENSION   AVVO2(IDIM1:IDIM2,LP1),HEATEM(IDIM1:IDIM2,LP1),
     1            OVER1D(IDIM1:IDIM2,LP1),
     1            TO31D(IDIM1:IDIM2,LP1),CONT1D(IDIM1:IDIM2,LP1),
     2            AVMO3(IDIM1:IDIM2,LP1),AVPHO3(IDIM1:IDIM2,LP1),
     2            C(IDIM1:IDIM2,LLP1),C2(IDIM1:IDIM2,LLP1)
      DIMENSION       ITOP(IDIM1:IDIM2),IBOT(IDIM1:IDIM2),
     &            INDTC(IDIM1:IDIM2)
      DIMENSION
     4 DELPTC(IDIM1:IDIM2),PTOP(IDIM1:IDIM2),PBOT(IDIM1:IDIM2),
     & FTOP(IDIM1:IDIM2),
     5 FBOT(IDIM1:IDIM2) ,EMSPEC(IDIM1:IDIM2,2)
C---DIMENSION OF VARIABLES EQUIVALENCED TO THOSE IN VTEMP---
      DIMENSION VTMP3(IDIM1:IDIM2,LP1),DSORC(IDIM1:IDIM2,LP1)
      DIMENSION ALP(IDIM1:IDIM2,LLP1),CSUB(IDIM1:IDIM2,LLP1),
     &  CSUB2(IDIM1:IDIM2,LLP1)
      DIMENSION FAC1(IDIM1:IDIM2,LP1)
      DIMENSION DELPR1(IDIM1:IDIM2,LP1),DELPR2(IDIM1:IDIM2,LP1)
      DIMENSION EMISDG(IDIM1:IDIM2,LP1),CONTDG(IDIM1:IDIM2,LP1),
     &  TO3DG(IDIM1:IDIM2,LP1)
      DIMENSION FLXNET(IDIM1:IDIM2,LP1)
      DIMENSION IXO(IDIM1:IDIM2,LP1)
      DIMENSION VSUM1(IDIM1:IDIM2,LP1)
      DIMENSION FLXTHK(IDIM1:IDIM2,LP1)
      DIMENSION Z1(IDIM1:IDIM2,LP1)
C---DIMENSION OF VARIABLES PASSED TO OTHER SUBROUTINES---
C   (AND NOT FOUND IN COMMON BLOCKS)
      DIMENSION E1CTS1(IDIM1:IDIM2,LP1),E1CTS2(IDIM1:IDIM2,L)
      DIMENSION E1CTW1(IDIM1:IDIM2,LP1),E1CTW2(IDIM1:IDIM2,L)
      DIMENSION EMD(IDIM1:IDIM2,LLP1),TPL(IDIM1:IDIM2,LLP1)
C   IT IS POSSIBLE TO EQUIVALENCE EMD,TPL TO THE ABOVE VARIABLES,
C   AS THEY GET CALLED AT DIFFERENT TIMES
      DIMENSION FXO(IDIM1:IDIM2,LP1),DT(IDIM1:IDIM2,LP1)
      DIMENSION FXOE2(IDIM1:IDIM2,LP1),DTE2(IDIM1:IDIM2,LP1)
      DIMENSION FXOSP(IDIM1:IDIM2,2),DTSP(IDIM1:IDIM2,2)
C
C     DIMENSION OF LOCAL VARIABLES
      DIMENSION RLOG(IDIM1:IDIM2,L),FLX(IDIM1:IDIM2,LP1)
      DIMENSION TOTEVV(IDIM1:IDIM2,LP1),CNTTAU(IDIM1:IDIM2,LP1)
C
      EQUIVALENCE (ALP,C,CSUB),(CSUB2,C2)
      EQUIVALENCE (FAC1,DSORC,OVER1D,DELPR2,FLXNET)
      EQUIVALENCE (DELPR1,HEATEM)
      EQUIVALENCE (IXO,AVVO2,FLXTHK,TO3DG)
      EQUIVALENCE (Z1,AVMO3,CONTDG)
      EQUIVALENCE (EMISDG,VSUM1,AVPHO3)
      EQUIVALENCE (EMD(IDIM1,1),E1CTS1(IDIM1,1)),
     &   (EMD(IDIM1,LP2),E1CTS2(IDIM1,1))
      EQUIVALENCE (TPL(IDIM1,1),E1CTW1(IDIM1,1)),
     &   (TPL(IDIM1,LP2),E1CTW2(IDIM1,1))
c
C
C          FIRST SECTION IS TABLE LOOKUP FOR SOURCE FUNCTION AND
C     DERIVATIVE (B AND DB/DT).ALSO,THE NLTE CO2 SOURCE FUNCTION
C     IS OBTAINED
C
C---IN CALCS. BELOW, DECREMENTING THE INDEX BY 9
C   ACCOUNTS FOR THE TABLES BEGINNING AT T=100K.
C   AT T=100K.
      DO 101 K=1,LP1
      DO 101 I=MYIS,MYIE
C---TEMP. INDICES FOR E1,SOURCE
      VTMP3(I,K)=AINT(TEMP(I,K)*HP1)
      FXO(I,K)=VTMP3(I,K)-9.
      DT(I,K)=TEMP(I,K)-TEN*VTMP3(I,K)
C---INTEGER INDEX FOR SOURCE (USED IMMEDIATELY)
      IXO(I,K)=FXO(I,K)
101   CONTINUE
      DO 103 k=1,L
      DO 103 I=MYIS,MYIE
C---TEMP. INDICES FOR E2 (KP=1 LAYER NOT USED IN FLUX CALCULATIONS)
      VTMP3(I,K)=AINT(T(I,K+1)*HP1)
      FXOE2(I,K)=VTMP3(I,K)-9.
      DTE2(I,K)=T(I,K+1)-TEN*VTMP3(I,K)
103   CONTINUE
C---SPECIAL CASE TO HANDLE KP=LP1 LAYER AND SPECIAL E2 CALCS.
      DO 105 I=MYIS,MYIE
      FXOE2(I,LP1)=FXO(I,L)
      DTE2(I,LP1)=DT(I,L)
      FXOSP(I,1)=FXOE2(I,LM1)
      FXOSP(I,2)=FXO(I,LM1)
      DTSP(I,1)=DTE2(I,LM1)
      DTSP(I,2)=DT(I,LM1)
105   CONTINUE
C
C---SOURCE FUNCTION FOR COMBINED BAND 1
      DO 4114 I=MYIS,MYIE
      DO 4114 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),1)
        DSORC(I,K)=DSRCE(IXO(I,K),1)
4114   CONTINUE
      DO 4112 K=1,LP1
      DO 4112 I=MYIS,MYIE
      SORC(I,K,1)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4112   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 2
      DO 4214 I=MYIS,MYIE
      DO 4214 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),2)
        DSORC(I,K)=DSRCE(IXO(I,K),2)
4214   CONTINUE
      DO 4212 K=1,LP1
      DO 4212 I=MYIS,MYIE
      SORC(I,K,2)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4212   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 3
      DO 4314 I=MYIS,MYIE
      DO 4314 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),3)
        DSORC(I,K)=DSRCE(IXO(I,K),3)
4314   CONTINUE
      DO 4312 K=1,LP1
      DO 4312 I=MYIS,MYIE
      SORC(I,K,3)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4312   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 4
      DO 4414 I=MYIS,MYIE
      DO 4414 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),4)
        DSORC(I,K)=DSRCE(IXO(I,K),4)
4414   CONTINUE
      DO 4412 K=1,LP1
      DO 4412 I=MYIS,MYIE
      SORC(I,K,4)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4412   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 5
      DO 4514 I=MYIS,MYIE
      DO 4514 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),5)
        DSORC(I,K)=DSRCE(IXO(I,K),5)
4514   CONTINUE
      DO 4512 K=1,LP1
      DO 4512 I=MYIS,MYIE
      SORC(I,K,5)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4512   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 6
      DO 4614 I=MYIS,MYIE
      DO 4614 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),6)
        DSORC(I,K)=DSRCE(IXO(I,K),6)
4614   CONTINUE
      DO 4612 K=1,LP1
      DO 4612 I=MYIS,MYIE
      SORC(I,K,6)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4612   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 7
      DO 4714 I=MYIS,MYIE
      DO 4714 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),7)
        DSORC(I,K)=DSRCE(IXO(I,K),7)
4714   CONTINUE
      DO 4712 K=1,LP1
      DO 4712 I=MYIS,MYIE
      SORC(I,K,7)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4712   CONTINUE
C---SOURCE FUNCTION FOR COMBINED BAND 8
      DO 4814 I=MYIS,MYIE
      DO 4814 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),8)
        DSORC(I,K)=DSRCE(IXO(I,K),8)
4814   CONTINUE
      DO 4812 K=1,LP1
      DO 4812 I=MYIS,MYIE
      SORC(I,K,8)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4812   CONTINUE
C---SOURCE FUNCTION FOR BAND 9 (560-670 CM-1)
      DO 4914 I=MYIS,MYIE
      DO 4914 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),9)
        DSORC(I,K)=DSRCE(IXO(I,K),9)
4914   CONTINUE
      DO 4912 K=1,LP1
      DO 4912 I=MYIS,MYIE
      SORC(I,K,9)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
4912   CONTINUE
C---SOURCE FUNCTION FOR BAND 10 (670-800 CM-1)
      DO 5014 I=MYIS,MYIE
      DO 5014 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),10)
        DSORC(I,K)=DSRCE(IXO(I,K),10)
5014  CONTINUE
      DO 5012 K=1,LP1
      DO 5012 I=MYIS,MYIE
      SORC(I,K,10)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
5012   CONTINUE
C---SOURCE FUNCTION FOR BAND 11 (800-900 CM-1)
      DO 5114 I=MYIS,MYIE
      DO 5114 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),11)
        DSORC(I,K)=DSRCE(IXO(I,K),11)
5114   CONTINUE
      DO 5112 K=1,LP1
      DO 5112 I=MYIS,MYIE
      SORC(I,K,11)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
5112   CONTINUE
C---SOURCE FUNCTION FOR BAND 12 (900-990 CM-1)
      DO 5214 I=MYIS,MYIE
      DO 5214 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),12)
        DSORC(I,K)=DSRCE(IXO(I,K),12)
5214   CONTINUE
      DO 5212 K=1,LP1
      DO 5212 I=MYIS,MYIE
      SORC(I,K,12)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
5212   CONTINUE
C---SOURCE FUNCTION FOR BAND 13 (990-1070 CM-1)
      DO 5314 I=MYIS,MYIE
      DO 5314 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),13)
        DSORC(I,K)=DSRCE(IXO(I,K),13)
5314   CONTINUE
      DO 5312 K=1,LP1
      DO 5312 I=MYIS,MYIE
      SORC(I,K,13)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
5312   CONTINUE
C---SOURCE FUNCTION FOR BAND 14 (1070-1200 CM-1)
      DO 5414 I=MYIS,MYIE
      DO 5414 K=1,LP1
        VTMP3(I,K)=SOURCE(IXO(I,K),14)
        DSORC(I,K)=DSRCE(IXO(I,K),14)
5414   CONTINUE
      DO 5412 K=1,LP1
      DO 5412 I=MYIS,MYIE
      SORC(I,K,14)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
5412   CONTINUE
C
C        THE FOLLOWING SUBROUTINE OBTAINS NLTE SOURCE FUNCTION FOR CO2
C
C
C     CALL NLTE
C
C
C---OBTAIN SPECIAL SOURCE FUNCTIONS FOR THE 15 UM BAND (CSOUR)
C   AND THE WINDOW REGION (SS1)
      DO 131 K=1,LP1
      DO 131 I=MYIS,MYIE
      SS1(I,K)=SORC(I,K,11)+SORC(I,K,12)+SORC(I,K,14)
131   CONTINUE
      DO 143 K=1,LP1
      DO 143 I=MYIS,MYIE
      CSOUR(I,K)=SORC(I,K,9)+SORC(I,K,10)
143   CONTINUE
C
C---COMPUTE TEMP**4 (TC) AND VERTICAL TEMPERATURE DIFFERENCES
C   (OSS,CSS,SS2,DTC). ALL THESE WILL BE USED LATER IN FLUX COMPUTA-
C   TIONS.
C
      DO 901 K=1,LP1
      DO 901 I=MYIS,MYIE
      TC(I,K)=(TEMP(I,K)*TEMP(I,K))**2
901   CONTINUE
      DO 903 K=1,L
      DO 903 I=MYIS,MYIE
      OSS(I,K+1)=SORC(I,K+1,13)-SORC(I,K,13)
      CSS(I,K+1)=CSOUR(I,K+1)-CSOUR(I,K)
      DTC(I,K+1)=TC(I,K+1)-TC(I,K)
      SS2(I,K+1)=SS1(I,K+1)-SS1(I,K)
903   CONTINUE
C
C
C---THE FOLLOWIMG IS A DRASTIC REWRITE OF THE RADIATION CODE TO
C    (LARGELY) ELIMINATE THREE-DIMENSIONAL ARRAYS. THE CODE WORKS
C    ON THE FOLLOWING PRINCIPLES:
C
C          LET K = FIXED FLUX LEVEL, KP = VARYING FLUX LEVEL
C          THEN FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K))
C               OVER ALL KPS, FROM 1 TO LP1.
C
C          WE CAN BREAK DOWN THE CALCULATIONS FOR ALL KS AS FOLLOWS:
C
C          FOR ALL KS K=1 TO LP1:
C              FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K))  (1)
C                      OVER ALL KPS, FROM K+1 TO LP1
C          AND
C              FOR KP FROM K+1 TO LP1:
C                 FLUX(KP) = DELTAB(K)*TAU(K,KP)              (2)
C
C          NOW IF TAU(K,KP)=TAU(KP,K) (SYMMETRICAL ARRAYS)
C          WE CAN COMPUTE A 1-DIMENSIONAL ARRAY TAU1D(KP) FROM
C          K+1 TO LP1, EACH TIME K IS INCREMENTED.
C          EQUATIONS (1) AND (2) THEN BECOME:
C
C             TAU1D(KP) = (VALUES FOR TAU(KP,K) AT THE PARTICULAR K)
C             FLUX(K) = SUM OVER KP : (DELTAB(KP)*TAU1D(KP))   (3)
C             FLUX(KP) = DELTAB(K)*TAU1D(KP)                   (4)
C
C         THE TERMS FOR TAU (K,K) AND OTHER SPECIAL TERMS (FOR
C         NEARBY LAYERS) MUST, OF COURSE, BE HANDLED SEPARATELY, AND
C         WITH CARE.
C
C      COMPUTE "UPPER TRIANGLE" TRANSMISSION FUNCTIONS FOR
C      THE 9.6 UM BAND (TO3SP) AND THE 15 UM BAND (OVER1D). ALSO,
C      THE
C      STAGE 1...COMPUTE O3 ,OVER TRANSMISSION FCTNS AND AVEPHI
C---DO K=1 CALCULATION (FROM FLUX LAYER KK TO THE TOP) SEPARATELY
C   AS VECTORIZATION IS IMPROVED,AND OZONE CTS TRANSMISSIVITY
C   MAY BE EXTRACTED HERE.
      DO 3021 K=1,L
      DO 3021 I=MYIS,MYIE
      AVEPHI(I,K)=TOTPHI(I,K+1)
3021  CONTINUE
C---IN ORDER TO PROPERLY EVALUATE EMISS INTEGRATED OVER THE (LP1)
C   LAYER, A SPECIAL EVALUATION OF EMISS IS DONE. THIS REQUIRES
C   A SPECIAL COMPUTATION OF AVEPHI, AND IT IS STORED IN THE
C   (OTHERWISE VACANT) LP1TH POSITION
C
      DO 803 I=MYIS,MYIE
      AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I)
803   CONTINUE
C   COMPUTE FLUXES FOR K=1
      CALL E1E290(E1CTS1,E1CTS2,E1FLX,E1CTW1,E1CTW2,EMISS,
     1            FXO,DT,FXOE2,DTE2,AVEPHI,TEMP,T)
      DO 302 K=1,L
      DO 302 I=MYIS,MYIE
      FAC1(I,K)=BO3RND(2)*TPHIO3(I,K+1)/TOTO3(I,K+1)
      TO3SPC(I,K)=HAF*(FAC1(I,K)*
     1    (SQRT(ONE+(FOUR*AO3RND(2)*TOTO3(I,K+1))/FAC1(I,K))-ONE))
C   FOR K=1, TO3SP IS USED INSTEAD OF TO31D (THEY ARE EQUAL IN THIS
C   CASE); TO3SP IS PASSED TO SPA90, WHILE TO31D IS A WORK-ARRAY.
      TO3SP(I,K)=EXP(HM1EZ*(TO3SPC(I,K)+SKO3R*TOTVO2(I,K+1)))
      OVER1D(I,K)=EXP(HM1EZ*(SQRT(AB15WD*TOTPHI(I,K+1))+
     1            SKC1R*TOTVO2(I,K+1)))
C---BECAUSE ALL CONTINUUM TRANSMISSIVITIES ARE OBTAINED FROM THE
C  2-D QUANTITY CNTTAU (AND ITS RECIPROCAL TOTEVV) WE STORE BOTH
C  OF THESE HERE. FOR K=1, CONT1D EQUALS CNTTAU
      CNTTAU(I,K)=EXP(HM1EZ*TOTVO2(I,K+1))
      TOTEVV(I,K)=1./CNTTAU(I,K)
302   CONTINUE
      DO 3022 K=1,L
      DO 3022 I=MYIS,MYIE
      CO2SP(I,K+1)=OVER1D(I,K)*CO21(I,1,K+1)
3022  CONTINUE
      DO 3023 K=1,L
      DO 3023 I=MYIS,MYIE
      CO21(I,K+1,1)=CO21(I,K+1,1)*OVER1D(I,K)
3023  CONTINUE
C---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION
      DO 1808 I=MYIS,MYIE
      RLOG(I,1)=OVER1D(I,1)*CO2NBL(I,1)
1808  CONTINUE
C---THE TERMS WHEN KP=1 FOR ALL K ARE THE PHOTON EXCHANGE WITH
C   THE TOP OF THE ATMOSPHERE, AND ARE OBTAINED DIFFERENTLY THAN
C   THE OTHER CALCULATIONS
      DO 305 K=2,LP1
      DO 305 I=MYIS,MYIE
      FLX(I,K)= (TC(I,1)*E1FLX(I,K)
     1          +SS1(I,1)*CNTTAU(I,K-1)
     2          +SORC(I,1,13)*TO3SP(I,K-1)
     3          +CSOUR(I,1)*CO2SP(I,K))
     4          *CLDFAC(I,1,K)
305   CONTINUE
      DO 307 I=MYIS,MYIE
      FLX(I,1)= TC(I,1)*E1FLX(I,1)+SS1(I,1)+SORC(I,1,13)
     1          +CSOUR(I,1)
307   CONTINUE
C---THE KP TERMS FOR K=1...
      DO 303 KP=2,LP1
      DO 303 I=MYIS,MYIE
      FLX(I,1)=FLX(I,1)+(OSS(I,KP)*TO3SP(I,KP-1)
     1                  +SS2(I,KP)*CNTTAU(I,KP-1)
     2                  +CSS(I,KP)*CO21(I,KP,1)
     3                  +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,1)
303   CONTINUE
C          SUBROUTINE SPA88 IS CALLED TO OBTAIN EXACT CTS FOR WATER
C     CO2 AND O3, AND APPROXIMATE CTS CO2 AND O3 CALCULATIONS.
C
      CALL SPA88(EXCTS,CTSO3,GXCTS,SORC,CSOUR,
     1           CLDFAC,TEMP,PRESS,VAR1,VAR2,
     2           P,DELP,DELP2,TOTVO2,TO3SP,TO3SPC,
     3           CO2SP1,CO2SP2,CO2SP)
C
C    THIS SECTION COMPUTES THE EMISSIVITY CTS HEATING RATES FOR 2
C    EMISSIVITY BANDS: THE 0-160,1200-2200 CM-1 BAND AND THE 800-
C    990,1070-1200 CM-1 BAND. THE REMAINING CTS COMTRIBUTIONS ARE
C    CONTAINED IN CTSO3, COMPUTED IN SPA88.
C
      DO 998 I=MYIS,MYIE
      VTMP3(I,1)=1.
998   CONTINUE
      DO 999 K=1,L
      DO 999 I=MYIS,MYIE
      VTMP3(I,K+1)=CNTTAU(I,K)*CLDFAC(I,K+1,1)
999   CONTINUE
      DO 1001 K=1,L
      DO 1001 I=MYIS,MYIE
      CTS(I,K)=RADCON*DELP(I,K)*(TC(I,K)*
     1     (E1CTW2(I,K)*CLDFAC(I,K+1,1)-E1CTW1(I,K)*CLDFAC(I,K,1)) +
     2      SS1(I,K)*(VTMP3(I,K+1)-VTMP3(I,K)))
1001  CONTINUE
C
      DO 1011 K=1,L
      DO 1011 I=MYIS,MYIE
      VTMP3(I,K)=TC(I,K)*(CLDFAC(I,K,1)*(E1CTS1(I,K)-E1CTW1(I,K)) -
     1                  CLDFAC(I,K+1,1)*(E1CTS2(I,K)-E1CTW2(I,K)))
1011  CONTINUE
      DO 1012 I=MYIS,MYIE
      FLX1E1(I)=TC(I,LP1)*CLDFAC(I,LP1,1)*
     1          (E1CTS1(I,LP1)-E1CTW1(I,LP1))
1012  CONTINUE
      DO 1014 K=1,L
      DO 1013 I=MYIS,MYIE
      FLX1E1(I)=FLX1E1(I)+VTMP3(I,K)
1013  CONTINUE
1014  CONTINUE
C
C---NOW REPEAT FLUX CALCULATIONS FOR THE K=2..LM1  CASES.
C   CALCULATIONS FOR FLUX LEVEL L AND LP1 ARE DONE SEPARATELY, AS ALL
C   EMISSIVITY AND CO2 CALCULATIONS ARE SPECIAL CASES OR NEARBY LAYERS.
C
      DO 321 K=2,LM1
      KLEN=K
C
      DO 3218 KK=1,LP1-K
      DO 3218 I=MYIS,MYIE
      AVEPHI(I,KK+K-1)=TOTPHI(I,KK+K)-TOTPHI(I,K)
3218  CONTINUE
      DO 1803 I=MYIS,MYIE
      AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I)
1803   CONTINUE
C---COMPUTE EMISSIVITY FLUXES (E2) FOR THIS CASE. NOTE THAT
C   WE HAVE OMITTED THE NEARBY LATER CASE (EMISS(I,K,K)) AS WELL
C   AS ALL CASES WITH K=L OR LP1. BUT THESE CASES HAVE ALWAYS
C   BEEN HANDLED AS SPECIAL CASES, SO WE MAY AS WELL COMPUTE
C    THEIR FLUXES SEPARASTELY.
C
      CALL E290(EMISSB,EMISS,AVEPHI,KLEN,FXOE2,DTE2)
      DO 322 KK=1,LP1-K
      DO 322 I=MYIS,MYIE
      AVMO3(I,KK+K-1)=TOTO3(I,KK+K)-TOTO3(I,K)
      AVPHO3(I,KK+K-1)=TPHIO3(I,KK+K)-TPHIO3(I,K)
      AVVO2(I,KK+K-1)=TOTVO2(I,KK+K)-TOTVO2(I,K)
      CONT1D(I,KK+K-1)=CNTTAU(I,KK+K-1)*TOTEVV(I,K-1)
322   CONTINUE
C
      DO 3221 KK=1,LP1-K
      DO 3221 I=MYIS,MYIE
      FAC1(I,K+KK-1)=BO3RND(2)*AVPHO3(I,K+KK-1)/AVMO3(I,K+KK-1)
      VTMP3(I,K+KK-1)=HAF*(FAC1(I,K+KK-1)*
     1  (SQRT(ONE+(FOUR*AO3RND(2)*AVMO3(I,K+KK-1))/
     2   FAC1(I,K+KK-1))-ONE))
      TO31D(I,K+KK-1)=EXP(HM1EZ*(VTMP3(I,K+KK-1)
     1                   +SKO3R*AVVO2(I,K+KK-1)))
      OVER1D(I,K+KK-1)=EXP(HM1EZ*(SQRT(AB15WD*AVEPHI(I,K+KK-1))+
     1            SKC1R*AVVO2(I,K+KK-1)))
      CO21(I,K+KK,K)=OVER1D(I,K+KK-1)*CO21(I,K+KK,K)
3221  CONTINUE
      DO 3223 KP=K+1,LP1
      DO 3223 I=MYIS,MYIE
      CO21(I,K,KP)=OVER1D(I,KP-1)*CO21(I,K,KP)
3223  CONTINUE
C---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION
      DO 1804 I=MYIS,MYIE
      RLOG(I,K)=OVER1D(I,K)*CO2NBL(I,K)
1804  CONTINUE
C---THE KP TERMS FOR ARBIRRARY K..
      DO 3423 KP=K+1,LP1
      DO 3423 I=MYIS,MYIE
      FLX(I,K)=FLX(I,K)+(OSS(I,KP)*TO31D(I,KP-1)
     1                  +SS2(I,KP)*CONT1D(I,KP-1)
     2                  +CSS(I,KP)*CO21(I,KP,K)
     3                  +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,K)
3423  CONTINUE
      DO 3425 KP=K+1,LP1
      DO 3425 I=MYIS,MYIE
      FLX(I,KP)=FLX(I,KP)+(OSS(I,K)*TO31D(I,KP-1)
     1                   +SS2(I,K)*CONT1D(I,KP-1)
     2                   +CSS(I,K)*CO21(I,K,KP)
     3                   +DTC(I,K)*EMISSB(I,KP-1))*CLDFAC(I,K,KP)
3425  CONTINUE
321   CONTINUE
C
C   NOW DO K=L CASE. SINCE THE KP LOOP IS LENGTH 1, MANY SIMPLIFI-
C   CATIONS OCCUR. ALSO, THE CO2 QUANTITIES (AS WELL AS THE EMISS
C  QUANTITIES) ARE COMPUTED IN THE NBL SEDCTION; THEREFORE, WE WANT
C  ONLY OVER,TO3 AND CONT1D (OVER(I,L),TO31D(I,L) AND CONT1D(I,L)
C  ACCORDING TO THE NOTATION. THUS NO CALL IS MADE TO THE E290
C  SUBROUTINE.
C         THE THIRD SECTION CALCULATES BOUNDARY LAYER AND NEARBY LAYER
C     CORRECTIONS TO THE TRANSMISSION FUNCTIONS OBTAINED ABOVE. METHODS
C     ARE GIVEN IN REF. (4).
C          THE FOLLOWING RATIOS ARE USED IN VARIOUS NBL CALCULATIONS:
C
C   THE REMAINING CALCULATIONS ARE FOR :
C                        1) THE (K,K) TERMS, K=2,LM1;
C                        2) THE (L,L) TERM
C                        3) THE (L,LP1) TERM
C                        4) THE (LP1,L) TERM
C                        5) THE (LP1,LP1) TERM.
C     EACH IS UNIQUELY HANDLED; DIFFERENT FLUX TERMS ARE COMPUTED
C     DIFFERENTLY
C
C
C          FOURTH SECTION OBTAINS WATER TRANSMISSION FUNCTIONS
C     USED IN Q(APPROX) CALCULATIONS AND ALSO MAKES NBL CORRECTIONS:
C     1) EMISS (I,J) IS THE TRANSMISSION FUNCTION MATRIX OBTAINED
C     BY CALLING SUBROUTINE E1E288;
C     2) "NEARBY LAYER" CORRECTIONS (EMISS(I,I)) ARE OBTAINED
C     USING SUBROUTINE E3V88;
C     3) SPECIAL VALUES AT THE SURFACE (EMISS(L,LP1),EMISS(LP1,L),
C     EMISS(LP1,LP1)) ARE CALCULATED.
C
C
C      OBTAIN ARGUMENTS FOR E1E288 AND E3V88:
C
      DO 821 I=MYIS,MYIE
      TPL(I,1)=TEMP(I,L)
      TPL(I,LP1)=HAF*(T(I,LP1)+TEMP(I,L))
      TPL(I,LLP1)=HAF*(T(I,L)+TEMP(I,L))
821   CONTINUE
      DO 823 K=2,L
      DO 823 I=MYIS,MYIE
      TPL(I,K)=T(I,K)
      TPL(I,K+L)=T(I,K)
823   CONTINUE
C
C---E2 FUNCTIONS ARE REQUIRED IN THE NBL CALCULATIONS FOR 2 CASES,
C   DENOTED (IN OLD CODE) AS (L,LP1) AND (LP1,LP1)
      DO 833 I=MYIS,MYIE
      AVEPHI(I,1)=VAR2(I,L)
      AVEPHI(I,2)=VAR2(I,L)+EMPL(I,L)
833   CONTINUE
      CALL E2SPEC(EMISS,AVEPHI,FXOSP,DTSP)
C
C     CALL E3V88 FOR NBL H2O TRANSMISSIVITIES
           CALL E3V88(EMD,TPL,EMPL)
C
C   COMPUTE NEARBY LAYER AND SPECIAL-CASE TRANSMISSIVITIES FOR EMISS
C    USING METHODS FOR H2O GIVEN IN REF. (4)
      DO 851 K=2,L
      DO 851 I=MYIS,MYIE
      EMISDG(I,K)=EMD(I,K+L)+EMD(I,K)
851   CONTINUE
C
C   NOTE THAT EMX1/2 (PRESSURE SCALED PATHS) ARE NOW COMPUTED IN
C   LWR88
      DO 861 I=MYIS,MYIE
      EMSPEC(I,1)=(EMD(I,1)*EMPL(I,1)-EMD(I,LP1)*EMPL(I,LP1))/
     1 EMX1(I) + QUARTR*(EMISS(I,1)+EMISS(I,2))
      EMISDG(I,LP1)=TWO*EMD(I,LP1)
      EMSPEC(I,2)=TWO*(EMD(I,1)*EMPL(I,1)-EMD(I,LLP1)*EMPL(I,LLP1))/
     * EMX2(I)
861   CONTINUE
      DO 331 I=MYIS,MYIE
      FAC1(I,L)=BO3RND(2)*VAR4(I,L)/VAR3(I,L)
      VTMP3(I,L)=HAF*(FAC1(I,L)*
     1    (SQRT(ONE+(FOUR*AO3RND(2)*VAR3(I,L))/FAC1(I,L))-ONE))
      TO31D(I,L)=EXP(HM1EZ*(VTMP3(I,L)+SKO3R*CNTVAL(I,L)))
      OVER1D(I,L)=EXP(HM1EZ*(SQRT(AB15WD*VAR2(I,L))+
     1            SKC1R*CNTVAL(I,L)))
      CONT1D(I,L)=CNTTAU(I,L)*TOTEVV(I,LM1)
      RLOG(I,L)=OVER1D(I,L)*CO2NBL(I,L)
331   CONTINUE
      DO 618 K=1,L
      DO 618 I=MYIS,MYIE
      RLOG(I,K)=LOG(RLOG(I,K))
618   CONTINUE
      DO 601 K=1,LM1
      DO 601 I=MYIS,MYIE
      DELPR1(I,K+1)=DELP(I,K+1)*(PRESS(I,K+1)-P(I,K+1))
      ALP(I,LP1+K-1)=-SQRT(DELPR1(I,K+1))*RLOG(I,K+1)
601   CONTINUE
      DO 603 K=1,L
      DO 603 I=MYIS,MYIE
      DELPR2(I,K+1)=DELP(I,K)*(P(I,K+1)-PRESS(I,K))
      ALP(I,K)=-SQRT(DELPR2(I,K+1))*RLOG(I,K)
603   CONTINUE
      DO 625 I=MYIS,MYIE
      ALP(I,LL)=-RLOG(I,L)
      ALP(I,LLP1)=-RLOG(I,L)*SQRT(DELP(I,L)*(P(I,LP1)-PRESS(I,LM1)))
625   CONTINUE
C        THE FIRST COMPUTATION IS FOR THE 15 UM BAND,WITH THE
C     FOR THE COMBINED H2O AND CO2 TRANSMISSION FUNCTION.
C
C       PERFORM NBL COMPUTATIONS FOR THE 15 UM BAND
C***THE STATEMENT FUNCTION SF IN PREV. VERSIONS IS NOW EXPLICITLY
C   EVALUATED.
      DO 631 K=1,LLP1
      DO 631 I=MYIS,MYIE
      C(I,K)=ALP(I,K)*(HMP66667+ALP(I,K)*(QUARTR+ALP(I,K)*HM6666M2))
631   CONTINUE
      DO 641 I=MYIS,MYIE
      CO21(I,LP1,LP1)=ONE+C(I,L)
      CO21(I,LP1,L)=ONE+(DELP2(I,L)*C(I,LL)-(PRESS(I,L)-P(I,L))*
     1 C(I,LLM1))/(P(I,LP1)-PRESS(I,L))
      CO21(I,L,LP1)=ONE+((P(I,LP1)-PRESS(I,LM1))*C(I,LLP1)-
     1 (P(I,LP1)-PRESS(I,L))*C(I,L))/(PRESS(I,L)-PRESS(I,LM1))
641   CONTINUE
      DO 643 K=2,L
      DO 643 I=MYIS,MYIE
      CO21(I,K,K)=ONE+HAF*(C(I,LM1+K)+C(I,K-1))
643   CONTINUE
C
C    COMPUTE NEARBY-LAYER TRANSMISSIVITIES FOR THE O3 BAND AND FOR THE
C    ONE-BAND CONTINUUM BAND (TO3 AND EMISS2). THE SF2 FUNCTION IS
C    USED. THE METHOD IS THE SAME AS DESCRIBED FOR CO2 IN REF (4).
      DO 651 K=1,LM1
      DO 651 I=MYIS,MYIE
      CSUB(I,K+1)=CNTVAL(I,K+1)*DELPR1(I,K+1)
      CSUB(I,LP1+K-1)=CNTVAL(I,K)*DELPR2(I,K+1)
651   CONTINUE
C---THE SF2 FUNCTION IN PREV. VERSIONS IS NOW EXPLICITLY EVALUATED
      DO 655 K=1,LLM2
      DO 655 I=MYIS,MYIE
      CSUB2(I,K+1)=SKO3R*CSUB(I,K+1)
      C(I,K+1)=CSUB(I,K+1)*(HMP5+CSUB(I,K+1)*
     1          (HP166666-CSUB(I,K+1)*H41666M2))
      C2(I,K+1)=CSUB2(I,K+1)*(HMP5+CSUB2(I,K+1)*
     1           (HP166666-CSUB2(I,K+1)*H41666M2))
655   CONTINUE
      DO 661 I=MYIS,MYIE
      CONTDG(I,LP1)=1.+C(I,LLM1)
      TO3DG(I,LP1)=1.+C2(I,LLM1)
661   CONTINUE
      DO 663 K=2,L
      DO 663 I=MYIS,MYIE
      CONTDG(I,K)=ONE+HAF*(C(I,K)+C(I,LM1+K))
      TO3DG(I,K)=ONE+HAF*(C2(I,K)+C2(I,LM1+K))
663   CONTINUE
C---NOW OBTAIN FLUXES
C
C    FOR THE DIAGONAL TERMS...
      DO 871 K=2,LP1
      DO 871 I=MYIS,MYIE
      FLX(I,K)=FLX(I,K)+(DTC(I,K)*EMISDG(I,K)
     1                 +SS2(I,K)*CONTDG(I,K)
     2                 +OSS(I,K)*TO3DG(I,K)
     3                 +CSS(I,K)*CO21(I,K,K))*CLDFAC(I,K,K)
871   CONTINUE
C     FOR THE TWO OFF-DIAGONAL TERMS...
      DO 873 I=MYIS,MYIE
      FLX(I,L)=FLX(I,L)+(CSS(I,LP1)*CO21(I,LP1,L)
     1                  +DTC(I,LP1)*EMSPEC(I,2)
     2                  +OSS(I,LP1)*TO31D(I,L)
     3                  +SS2(I,LP1)*CONT1D(I,L))*CLDFAC(I,LP1,L)
      FLX(I,LP1)=FLX(I,LP1)+(CSS(I,L)*CO21(I,L,LP1)
     1                      +OSS(I,L)*TO31D(I,L)
     2                      +SS2(I,L)*CONT1D(I,L)
     3                      +DTC(I,L)*EMSPEC(I,1))*CLDFAC(I,L,LP1)
873   CONTINUE
C
C     FINAL SECTION OBTAINS EMISSIVITY HEATING RATES,
C     TOTAL HEATING RATES AND THE FLUX AT THE GROUND
C
C     .....CALCULATE THE EMISSIVITY HEATING RATES
      DO 1101 K=1,L
      DO 1101 I=MYIS,MYIE
      HEATEM(I,K)=RADCON*(FLX(I,K+1)-FLX(I,K))*DELP(I,K)
1101  CONTINUE
C     .....CALCULATE THE TOTAL HEATING RATES
      DO 1103 K=1,L
      DO 1103 I=MYIS,MYIE
      HEATRA(I,K)=HEATEM(I,K)-CTS(I,K)-CTSO3(I,K)+EXCTS(I,K)
1103  CONTINUE
C     .....CALCULATE THE FLUX AT EACH FLUX LEVEL USING THE FLUX AT THE
C    TOP (FLX1E1+GXCTS) AND THE INTEGRAL OF THE HEATING RATES (VSUM1)
      DO 1111 K=1,L
      DO 1111 I=MYIS,MYIE
      VSUM1(I,K)=HEATRA(I,K)*DELP2(I,K)*RADCON1
1111  CONTINUE
      DO 1115 I=MYIS,MYIE
      TOPFLX(I)=FLX1E1(I)+GXCTS(I)
      FLXNET(I,1)=TOPFLX(I)
1115  CONTINUE
C---ONLY THE SURFACE VALUE OF FLUX (GRNFLX) IS NEEDED UNLESS
C    THE THICK CLOUD SECTION IS INVOKED.
      DO 1123 K=2,LP1
      DO 1123 I=MYIS,MYIE
      FLXNET(I,K)=FLXNET(I,K-1)+VSUM1(I,K-1)
1123  CONTINUE
      DO 1125 I=MYIS,MYIE
      GRNFLX(I)=FLXNET(I,LP1)
1125  CONTINUE
C
C     THIS IS THE THICK CLOUD SECTION.OPTIONALLY,IF THICK CLOUD
C     FLUXES ARE TO BE "CONVECTIVELY ADJUSTED",IE,DF/DP IS CONSTANT,
C     FOR CLOUDY PART OF GRID POINT, THE FOLLOWING CODE IS EXECUTED.
C***FIRST,COUNT THE NUMBER OF CLOUDS ALONG THE LAT. ROW. SKIP THE
C   ENTIRE THICK CLOUD COMPUTATION OF THERE ARE NO CLOUDS.
      ICNT=0
      DO 1301 I=MYIS,MYIE
      ICNT=ICNT+NCLDS(I)
1301  CONTINUE
      IF (ICNT.EQ.0) GO TO 6999
C---FIND THE MAXIMUM NUMBER OF CLOUDS IN THE LATITUDE ROW
      KCLDS=NCLDS(1)
      DO 2106 I=MYIS,MYIE
      KCLDS=MAX(NCLDS(I),KCLDS)
2106  CONTINUE
C
C
C***OBTAIN THE PRESSURES AND FLUXES OF THE TOP AND BOTTOM OF
C   THE NC'TH CLOUD (IT IS ASSUMED THAT ALL KTOP AND KBTM'S HAVE
C   BEEN DEFINED!).
      DO 1361 KK=1,KCLDS
      KMIN=LP1
      KMAX=0
      DO 1362 I=MYIS,MYIE
        J1=KTOP(I,KK+1)
C       IF (J1.EQ.1) GO TO 1362
        J3=KBTM(I,KK+1)
        IF (J3.GT.J1) THEN
          PTOP(I)=P(I,J1)
          PBOT(I)=P(I,J3+1)
          FTOP(I)=FLXNET(I,J1)
          FBOT(I)=FLXNET(I,J3+1)
C***OBTAIN THE "FLUX DERIVATIVE" DF/DP (DELPTC)
          DELPTC(I)=(FTOP(I)-FBOT(I))/(PTOP(I)-PBOT(I))
          KMIN=MIN(KMIN,J1)
          KMAX=MAX(KMAX,J3)
        ENDIF
1362  CONTINUE
      KMIN=KMIN+1
C***CALCULATE THE TOT. FLUX CHG. FROM THE TOP OF THE CLOUD, FOR
C   ALL LEVELS.
      DO 1365 K=KMIN,KMAX
      DO 1363 I=MYIS,MYIE
C       IF (KTOP(I,KK+1).EQ.1) GO TO 1363
        IF(KTOP(I,KK+1).LT.K .AND. K.LE.KBTM(I,KK+1)) THEN
          Z1(I,K)=(P(I,K)-PTOP(I))*DELPTC(I)+FTOP(I)
CORIGINAL FLXNET(I,K)=FLXNET(I,K)*(ONE-CAMT(I,KK+1)) +
CORIGINAL1            Z1(I,K)*CAMT(I,KK+1)
          FLXNET(I,K)=Z1(I,K)
        ENDIF
1363  CONTINUE
1365  CONTINUE
1361  CONTINUE
C***USING THIS FLUX CHG. IN THE CLOUDY PART OF THE GRID BOX, OBTAIN
C   THE NEW FLUXES, WEIGHTING THE CLEAR AND CLOUDY FLUXES:AGAIN, ONLY
C    THE FLUXES IN THICK-CLOUD LEVELS WILL EVENTUALLY BE USED.
C     DO 6051 K=1,LP1
C     DO 6051 I=MYIS,MYIE
C     FLXNET(I,K)=FLXNET(I,K)*(ONE-CAMT(I,NC)) +
C    1            Z1(I,K)*CAMT(I,NC)
C051  CONTINUE
C***MERGE FLXTHK INTO FLXNET FOR APPROPRIATE LEVELS.
C     DO 1401 K=1,LP1
C     DO 1401 I=MYIS,MYIE
C     IF (K.GT.ITOP(I) .AND. K.LE.IBOT(I)
C    1  .AND.  (NC-1).LE.NCLDS(I))  THEN
C          FLXNET(I,K)=FLXTHK(I,K)
C     ENDIF
C401  CONTINUE
C
C******END OF CLOUD LOOP*****
6001  CONTINUE
6999  CONTINUE
C***THE FINAL STEP IS TO RECOMPUTE THE HEATING RATES BASED ON THE
C   REVISED FLUXES:
      DO 6101 K=1,L
      DO 6101 I=MYIS,MYIE
      HEATRA(I,K)=RADCON*(FLXNET(I,K+1)-FLXNET(I,K))*DELP(I,K)
6101  CONTINUE
C     THE THICK CLOUD SECTION ENDS HERE.
      RETURN
      END
