      SUBROUTINE LONGW2O(C MTX, EMD, EMM, IL1, IL2, LON, LAY, LCLW,
     1                  DSIG, PSOL, WV, TAVE, TAUAE, CLT, CLB,
     2                  T, O, RESULT, F, BINT, BSOLIN, DT0, EM0, DZ,
     3                  SSIG, XADJD, XADJM, XDESC, XMONT, 
     4                  U, V, W,
     5                  FDN, FUP, 
     6                  TRTOP, CNTRB, FDES, FMON, 
     7                  BSOL, XN, XD, XT, 
     8                  UU, TT1, TT2, TT, 
     9                  BTOP, XQ, XO, TAVIC,
     A                  TI, RES, XDG, XMG,
     B                  BGND, FM, FD, ZZ, 
     C                  DUC, UAER, FACTC, EPSH1,
     D                  EPSH2, EPSH3, EPSH4, EPSH5, 
     E                  EPSC2, EPSC3, EPSC4, EPSIO, 
     F                  EPHIO, DSTSOL, DSTTOP, DST, 
     G                  DTS, IXTSOL, IXTTOP, IXT, 
     H                  ITX)
  
C     * MAY  17/95 - M.LAZARE. - CALCULATE AND PASS OUT ADDITIONAL 
C     *                          ARRAYS "CLT" AND "CLB" TO BE USED
C     *                          IN ACCUMULATING CLEAR-SKY L/W FLUXES.
C     * SEPT 23/88 - M.LAZARE. - PREVIOUS VERSION LONGW2.
C 
C***********************************************************************
C                                                                      C
C   DEVELOPED, 'AEROSOLIZED' AND VECTORIZED  (ORDER ILX) BY            C
C                                                                      C
C    JEAN-JACQUES MORCRETTE                                            C
C    LABORATOIRE D'OPTIQUE ATMOSPHERIQUE                               C
C    UNIVERSITE DES SCIENCES ET TECHNIQUES DE LILLE                    C
C    59655 - VILLENEUVE D'ASCQ CEDEX, FRANCE                           C
C                                                                      C
C    THESE ETAT, 3E CHAPITRE, JUIN 1984                                C
C    SUBMITTED TO QUART.J.ROY.METEOR.SOC. (PAPER 84/64)                C
C                                                                      C
C-----------------------------------------------------------------------
  
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      PARAMETER ( NG1=2, NINT=5, NTR=11,
     A            NG1P1=NG1+1, NUA=3*NINT+3, MP=NINT+3, 
     B            MP2=2*MP, NTRA=NUA+1, LN6=MP2*NTR*6 ) 
  
      COMMON /RADSIZE/ KMX, LMX, KMXP, NGLP1, NGL 
  
      COMMON /AERO/ TAUA(2,5),SSC ALB(2,5),CGA(2,5),CAER(5,5) 
      COMMON /ABSOR/  IABS,IAER 
      COMMON /PARAD1/ RT1(NG1),WG1(NG1),MPAD,AT(MP,5),BT(MP,5),CNTNU(5),
     1                O1H,O2H,OCT(4),TREF,XP(2,6,NINT),TSTAND,
     2                GA(MP2,NTR,6),GB(MP2,NTR,6),NPAD, 
     3                TINTP(NTR),TSTP,MXIXT 
      COMMON /CONSTN/ DAY, C DYN MB, ERG W M2, PI,
     2                STFN, DIFFU
      COMMON /LWCONS/ MAXC, CO2, SIG1 
      COMMON /SWCONS/ CSEC, CH2O, CCAR, CCO2, FCST, RAYL
  
C........ ARRAYS SHARED WITH RADIATN AND SHORTW 
  
      DIMENSION DSIG(LMX,KMX),PSOL(LMX) 
      DIMENSION WV(LMX,KMX),TAVE(LMX,KMX),TAUAE(LMX,KMX,5)
      DIMENSION CLT(LMX),CLB(LMX)
  
C........ ARRAYS SHARED WITH RADIATN ROUTINE ONLY 
  
      DIMENSION T(LMX,KMXP),O(LMX,KMXP),RESULT(LMX,KMX),F(LMX,2,KMXP),
     A BINT(LMX,KMXP),BSOLIN(LMX),DT0(LMX),EM0(LMX),DZ(LMX,KMX),
     B SSIG(LMX,NGLP1),XADJD(LMX,KMXP),XADJM(LMX,KMXP),XDESC(LMX,KMXP), 
     C XMONT(LMX,KMXP)
  
C........ VARIABLES USED IN LONGWAVE ROUTINE ONLY 
  
      DIMENSION U(LMX,8,NGLP1), V(LMX,NUA,NGLP1), 
     A W(LMX,NINT,NGLP1), FDN(LMX,KMXP,KMXP), FUP(LMX,KMXP,KMXP), 
     B TRTOP(LMX,KMXP), 
     C CNTRB(LMX,KMXP,KMXP),FDES(LMX,KMXP),FMON(LMX,KMXP),
     D BSOL(LMX,NINT),XN(LMX),XD(LMX),XT(LMX,NGLP1),
     E UU(LMX,NUA),TT1(LMX,NTRA),TT2(LMX,NTRA),TT(LMX,NTRA),
     F BTOP(LMX,NINT),XQ(LMX),XO(LMX),TAVIC(LMX),TI(LMX), 
     G RES(LMX),XDG(LMX),XMG(LMX),BGND(LMX),FM(LMX),FD(LMX),
     H ZZ(LMX),DUC(LMX,NGLP1),UAER(LMX,NINT),FACTC(LMX),EPSH1(LMX), 
     I EPSH2(LMX),EPSH3(LMX),EPSH4(LMX),EPSH5(LMX),EPSC2(LMX),
     J EPSC3(LMX),EPSC4(LMX),EPSIO(LMX),EPHIO(LMX), 
     K DSTSOL(LMX),DSTTOP(LMX),DST(LMX,KMXP),DTS(LMX),IXTSOL(LMX),
     L IXTTOP(LMX),IXT(LMX,KMXP),ITX(LMX) 
  
C........ I/O ARRAYS
  
      DIMENSION C MTX(LMX,KMXP,KMXP), EMD(LMX,KMX), EMM(LMX,KMX)
  
C........ THE FOLLOWING ARE EQUIVALENCED IN THE CALL TO LONGW.
C         NOTE THAT ARRAY "W" WILL ALWAYS BE HIDDEN INSIDE ARRAY "U", 
C         WHILE ARRAY "FDN" WILL BE HIDDEN INSIDE ARRAY "V" AS LONG AS
C         KMX.LE.52 (I.E. LESS THAN 51 FULL SIGMA LEVELS) 
C         SIMILARILY, ARRAY "FUP" WILL ALWAYS BE HIDDEN INSIDE ARRAY "U", 
C         AS LONG AS KMX.LE.23 (I.E. LESS THAN 22 FULL SIGMA LEVELS)
  
C     EQUIVALENCE (U(1,1,1), W(1,1,1), FUP(1,1,1)), 
C    1            (V(1,1,1), FDN(1,1,1)), 
C    2            (CNTRB(1,1,1), DUC(1,1)), 
C    3            (XT(1,1), FDES(1,1)), 
C    4            (FMON(1,1), IXT(1,1)),
C    5            (BTOP(1,1), UAER(1,1))
  
      LOGICAL LO1, LCLW 
  
C-----------------------------------------------------------------------
  
      IFTT=0
      IAER=0
      IABS=5
      KMXM1=KMX-1 
  
      IF (LCLW) THEN
C 
C     * INITIALIZE TRANSMISSIVITY ARRAYS FOR HORNEL.
C 
      DO 110 IA=1,NTRA
          DO 100 I=IL1,IL2
              TT(I,IA)=1. 
  100     CONTINUE
  110 CONTINUE
C 
C     THE PARAMETER WXT CONTROLS HOW THE TEMPERATURE GRADIENTS USED IN
C     THE FULL RADIATION CALCULATION ARE DEFINED. WXT=0 IMPLIES THAT
C     ONLY THE INFORMATION ON THE FULL TEMPERATURE LEVELS ARE USED. 
C     LARGER VALUES OF WXT IMPLY THAT PROGRESSIVELY MORE OF THE 
C     INFORMATION AT HALF LEVELS ARE USED AS WELL. VARIOUS TESTS OVER 
C     PARAMETER SPACE HAVE SHOWN THAT VALUES LARGER THAN 0.05 RESULT IN 
C     PROGRESSIVELY COLDER STRATOSPHERIC TEMPERATURES AND GREATER 
C     INSTABILITY AS WXT GETS LARGER. THERE SEEMS TO BE LITTLE DIFFERENCE 
C     BETWEEN VALUES OF 0. AND 0.05; HENCE THE LARGER VALUE IS USED TO
C     CONSIDER AS MUCH OF THE ACTUAL TEMPERATURE DATA (ON HALF LEVELS) AS 
C     POSSIBLE. 
C 
      WXT  = 0.05 
      WXTI = 1./(1.+WXT)
  
C........ CALCULATES INDEX OF THE TRANSMISSIVITY TO BE USED AND THE 
C         CORRESPONDING LINEAR WEIGHTING FACTOR 
  
      DO 120 I=IL1,IL2
          IXTS=MIN(MXIXT, INT((T(I,1)+DT0(I)-TINTP(1))/TSTP)+1)
          IXTSOL(I)=MAX( 1, IXTS ) 
          IXTT=MIN(MXIXT, INT((T(I,KMXP)    -TINTP(1))/TSTP)+1)
          IXTTOP(I)=MAX( 1, IXTT ) 
          DSTSOL(I)=(T(I,1)+DT0(I)-TINTP(IXTSOL(I)))/TSTP 
          DSTTOP(I)=(T(I,KMXP)    -TINTP(IXTTOP(I)))/TSTP 
  120 CONTINUE
  
      DO 122 K=1,KMX
          DO 121 I=IL1,IL2
C            IXTX=MIN(MXIXT, INT((TAVE(I,K)-TINTP(1))/TSTP)+1) 
C            IXT(I,K)=MAX( 1, IXTX ) 
C            DST(I,K)=(TAVE(I,K)-TINTP(IXT(I,K))))/TSTP 
             IXTX=MAX(1,MIN(MXIXT,INT((TAVE(I,K)-TINTP(1))/TSTP)+1))
             DST(I,K)=(TAVE(I,K)-TINTP(IXTX))/TSTP
C**          IXT(I,K) = CVMGT (IXTX, IXTX + 1, DST(I,K) .LE. 0.5) 
             IF(DST(I,K).LE.0.5)      THEN
                IXT(I,K)=IXTX
             ELSE
                IXT(I,K)=IXTX+1
             ENDIF 
             DST(I,K) = 0.
  121     CONTINUE
  122 CONTINUE
  
C........ CALCULATES THE TEMPERATURE GRADIENTS, MIXING RATIOS ON
C         SUB-LEVELS. INTERPOLATES WITH RESPECT TO SCALE HEIGHT 
C         Z=-LN(SIGMA)
  
C........ TEMPERATURE GRADIENTS 
  
  
      DO 132  K=1,KMX 
          KP1=K+1 
          KJ=(K-1)*NG1P1+1
          KJ1=KJ+1
          KJ2=KJ1+1 
          DO 130 I=IL1,IL2
            XTLO  = (T(I,KP1)-T(I,K))/(2.*DZ(I,K))
            XTHI1 = (TAVE(I,K)-T(I,K))/DZ(I,K)
            XTHI2 = (T(I,KP1)-TAVE(I,K))/DZ(I,K)
            XT(I,KJ1)=(XTLO+WXT*XTHI1)*WXTI 
            XT(I,KJ2)=(XTLO+WXT*XTHI2)*WXTI 
  130     CONTINUE
  132 CONTINUE
  
C........ ABSORBING AMOUNTS IN THE SUB-LAYERS 
C         COEF=PSOL(PA)/9.80665/10  (10= KG.M-2 / G.CM-2) 
C         SUB-LAYERS SIGMA THICKNESS AND MEAN SIGMA PRESSURE
  
      DO 134 KI=1,NGL 
          KIP1=KI+1 
          DO 133 I=IL1,IL2
              U(I,5,KI)=(SSIG(I,KI)+SSIG(I,KIP1))*0.5 
              U(I,3,KI)=(SSIG(I,KI)-SSIG(I,KIP1))*PSOL(I)/98.0665 
  133     CONTINUE
  134 CONTINUE
C 
      DO 138 K=1,KMX
          KP1=K+1 
          KJ=(K-1)*NG1P1+1
          KJPN=KJ+NG1 
          DO 135 I=IL1,IL2
              XQ(I)= WV(I,K)
              XO(I)= O(I,K) 
  135     CONTINUE
          DO 137 KK=KJ,KJPN 
              DO 136 I=IL1,IL2
                  UPM=U(I,5,KK)*PSOL(I)/101325. 
                  DPM=U(I,3,KK) 
                  DUC(I,KK)=DPM 
                  U(I,1,KK)=XO(I)*DPM 
                  U(I,2,KK)=XO(I)*DPM*UPM 
                  XXU=XQ(I)*DPM 
                  U6=XXU*UPM
C . . . . . .     28.9647 / 18.0153 
                  FPPW= 1.6078 *XQ(I)/(1.+0.608*XQ(I))
                  U(I,6,KK)=U6
                  U(I,7,KK)=U6*FPPW 
                  U(I,8,KK)=U6*(1.-FPPW)
                  XXC=CO2*DPM 
                  XXCP=XXC*UPM
                  U(I,3,KK)=XXC 
                  U(I,4,KK)=XXCP
  136         CONTINUE
  137     CONTINUE
  138 CONTINUE
  
C........ ARRAY V IS FILLED WITH THE ABSORBING QUANTITIES CUMULATED 
C         FROM THE TOP OF THE ATMOSPHERE (SET TO ZERO AT THE TOP).
  
      DO 141 KI=1,NUA 
          DO 140 I=IL1,IL2
              V(I,KI,NGLP1)=0.
  140     CONTINUE
  141 CONTINUE
  
      DO 148 K=1,KMX
          KJ=(K-1)*NG1P1+1
          KJPN=KJ+NG1 
          L=KMXP-K
  
C........ CALCULATES THE AMOUNTS OF AEROSOLS: IN THE FIVE SPECTRAL
C         INTERVALS CONSIDERED IN THE LONGWAVE SCHEME, CALCULATES 
C         EXTINCTION COEFFICIENTS FOR A MIXING OF THE 5 MODELS OF 
C         STANDARD AEROSOLS OF THE IAMAP-INT.RADIATION COMMISSION.
C             TAUAE = OPTICAL THICKNESS 
C             CAER = ABSORPTION COEFFICIENT NORMALIZED AT 0.55 MICRONS
  
          JAE1=NGLP1-KJ 
          JAE2=NGLP1-(KJ+1) 
          JAE3=NGLP1-KJPN 
          DO 144 IAE=1,5
C........ MOD (JPB) DEC 27, 85: FOR INITIALISATION. 
              CAER(IAE,5) = 0.
              DO 143 I=IL1,IL2
                  UAER(I,IAE) = (CAER(IAE,1)*TAUAE(I,L,1)+CAER(IAE,2)*
     1                          TAUAE(I,L,2)+CAER(IAE,3)*TAUAE(I,L,3)+
     2                          CAER(IAE,4)*TAUAE(I,L,4)+CAER(IAE,5)* 
     3                          TAUAE(I,L,5))/(DUC(I,JAE1)+DUC(I,JAE2)+ 
     4                          DUC(I,JAE3))
  143         CONTINUE
  144     CONTINUE
  
C........ INTRODUCES THE TEMPERATURE EFFECT ON THE ABSORBING AMOUNTS AND
C         CALCULATES THE ABSORBING QUANTITIES CUMULATED FROM THE TOP OF 
C         THE ATMOSPHERE
  
      IF (IFTT .NE. 0) THEN 
          DO 1441 I=IL1,IL2 
              EPSH1(I)=1. 
              EPSH2(I)=1. 
              EPSH3(I)=1. 
              EPSH4(I)=1. 
              EPSH5(I)=1. 
              EPSC2(I)=1. 
              EPSC3(I)=1. 
              EPSC4(I)=1. 
              FACTC(I)=1. 
              EPHIO(I)=1. 
              EPSIO(I)=1. 
 1441     CONTINUE
      ELSE
        DO 145 I=IL1,IL2
          TAVIC(I)=TAVE(I,L)
C . . . . TEMPERATURE DEPENDENCE OF THE H2O CONTINUUM ABSORPTION IN 
C         THE WINDOW
          FACTC(I)=EXP(6.08*(296./TAVIC(I)-1.)) 
          TX=TAVIC(I)-TREF
          TX2 = TX * TX 
C . . . . EPSH1  H2O  MP  0-500 CM-1
          UP = MIN (MAX (0.5 * LOG10(U(I,6,L)) + 5., 0.), 6.)
          CAH1=AT(1,1)+UP*(AT(1,2)+UP*(AT(1,3)+UP*(AT(1,4)+UP*AT(1,5))))
          CBH1=BT(1,1)+UP*(BT(1,2)+UP*(BT(1,3)+UP*(BT(1,4)+UP*BT(1,5))))
          PSI1=CAH1*TX+CBH1*TX2 
          EPSH1(I)=EXP(PSI1)
C . . . . EPSH2  H2O  MP  500-800 CM-1
          CAH2=AT(2,1)+UP*(AT(2,2)+UP*(AT(2,3)+UP*(AT(2,4)+UP*AT(2,5))))
          CBH2=BT(2,1)+UP*(BT(2,2)+UP*(BT(2,3)+UP*(BT(2,4)+UP*BT(2,5))))
          PSI2=CAH2*TX+CBH2*TX2 
          EPSH2(I)=EXP(PSI2)
C . . . . EPSH3  H2O  MP  800-970 + 1110-1250 CM-1
          CAH3=AT(3,1)+UP*(AT(3,2)+UP*(AT(3,3)+UP*(AT(3,4)+UP*AT(3,5))))
          CBH3=BT(3,1)+UP*(BT(3,2)+UP*(BT(3,3)+UP*(BT(3,4)+UP*BT(3,5))))
          PSI3=CAH3*TX+CBH3*TX2 
          EPSH3(I)=EXP(PSI3)
C . . . . EPSH4  H2O  MP  970-1110 CM-1 
          CAH4=AT(4,1)+UP*(AT(4,2)+UP*(AT(4,3)+UP*(AT(4,4)+UP*AT(4,5))))
          CBH4=BT(4,1)+UP*(BT(4,2)+UP*(BT(4,3)+UP*(BT(4,4)+UP*BT(4,5))))
          PSI4=CAH4*TX+CBH4*TX2 
          EPSH4(I)=EXP(PSI4)
C . . . . EPSH5  H2O  MP  1250-2820 CM-1
          CAH5=AT(5,1)+UP*(AT(5,2)+UP*(AT(5,3)+UP*(AT(5,4)+UP*AT(5,5))))
          CBH5=BT(5,1)+UP*(BT(5,2)+UP*(BT(5,3)+UP*(BT(5,4)+UP*BT(5,5))))
          PSI5=CAH5*TX+CBH5*TX2 
          EPSH5(I)=EXP(PSI5)
  145   CONTINUE
        DO 1451 I=IL1,IL2 
          TX = TAVE(I,L) - TREF 
          TX2 = TX * TX 
          UP = 0.5 * LOG10(U(I,4,L)) 
C . . . . EPSC2  CO2  MP  500-800 CM-1
          CAC6=AT(6,1)+UP*(AT(6,2)+UP*(AT(6,3)+UP*(AT(6,4)+UP*AT(6,5))))
          CBC6=BT(6,1)+UP*(BT(6,2)+UP*(BT(6,3)+UP*(BT(6,4)+UP*BT(6,5))))
          PSC2=CAC6*TX+CBC6*TX2 
          EPSC2(I)=EXP(PSC2)
C . . . . EPSC3  CO2  MP  800-970 + 1110-1250 CM-1
          CAC7=AT(7,1)+UP*(AT(7,2)+UP*(AT(7,3)+UP*(AT(7,4)+UP*AT(7,5))))
          CBC7=BT(7,1)+UP*(BT(7,2)+UP*(BT(7,3)+UP*(BT(7,4)+UP*BT(7,5))))
          PSC3=CAC7*TX+CBC7*TX2 
          EPSC3(I)=EXP(PSC3)
C . . . . EPSC4  CO2  MP  970-1110 CM-1 
          CAC8=AT(8,1)+UP*(AT(8,2)+UP*(AT(8,3)+UP*(AT(8,4)+UP*AT(8,5))))
          CBC8=BT(8,1)+UP*(BT(8,2)+UP*(BT(8,3)+UP*(BT(8,4)+UP*BT(8,5))))
          PSC4=CAC8*TX+CBC8*TX2 
          EPSC4(I)=EXP(PSC4)
C . . . . EPHIO   O3  M  970-1110 CM-1
          PHIO=OCT(1)*TX+OCT(2)*TX2 
          PSIO=OCT(3)*TX+OCT(4)*TX2 
          EPHIO(I)=EXP(PHIO)
C . . . . EPSIO   O3  MP 970-1110 CM-1
          EPSIO(I)=EXP(2.*PSIO) 
 1451   CONTINUE
      END IF
  
          DO 147 KK=KJ,KJPN 
              KC=NGLP1-KK 
              KCP1=KC+1 
              DO 146 I=IL1,IL2
                  V(I,1,KC)=V(I,1,KCP1)+U(I,6,KC)*EPSH1(I)*DIFFU
                  V(I,2,KC)=V(I,2,KCP1)+U(I,6,KC)*EPSH2(I)*DIFFU
                  V(I,3,KC)=V(I,3,KCP1)+U(I,6,KC)*EPSH3(I)*DIFFU
                  V(I,4,KC)=V(I,4,KCP1)+U(I,6,KC)*EPSH4(I)*DIFFU
                  V(I,5,KC)=V(I,5,KCP1)+U(I,6,KC)*EPSH5(I)*DIFFU
  
                  V(I,6,KC)=V(I,6,KCP1)+U(I,4,KC)*EPSC2(I)*DIFFU
                  V(I,7,KC)=V(I,7,KCP1)+U(I,4,KC)*EPSC3(I)*DIFFU
                  V(I,8,KC)=V(I,8,KCP1)+U(I,4,KC)*EPSC4(I)*DIFFU
  
                  V(I,9,KC)=V(I,9,KCP1)+U(I,8,KC)           *DIFFU
                  V(I,10,KC)=V(I,10,KCP1)+U(I,7,KC)         *DIFFU
                  V(I,11,KC)=V(I,11,KCP1)+U(I,7,KC)*FACTC(I)*DIFFU
  
                  V(I,12,KC)=V(I,12,KCP1)+U(I,1,KC)*EPHIO(I)*DIFFU
                  V(I,13,KC)=V(I,13,KCP1)+U(I,2,KC)*EPSIO(I)*DIFFU
  
                  V(I,14,KC)=V(I,14,KCP1)+UAER(I,1)*DUC(I,KC)*DIFFU 
                  V(I,15,KC)=V(I,15,KCP1)+UAER(I,2)*DUC(I,KC)*DIFFU 
                  V(I,16,KC)=V(I,16,KCP1)+UAER(I,3)*DUC(I,KC)*DIFFU 
                  V(I,17,KC)=V(I,17,KCP1)+UAER(I,4)*DUC(I,KC)*DIFFU 
                  V(I,18,KC)=V(I,18,KCP1)+UAER(I,5)*DUC(I,KC)*DIFFU 
  146         CONTINUE
  147     CONTINUE
  148 CONTINUE
  
C........ ARRAY CONTAINS PRESSURE AND TEMP SCALED AMOUNTS OF ABSORBERS
C             V(1,   H2O    FOR 1 ST SPECTRAL INTERVAL
C             V(2,   H2O    FOR 2 ND SPECTRAL INTERVAL
C             V(3,   H2O    FOR 3 RD SPECTRAL INTERVAL
C             V(4,   H2O    FOR 4 TH SPECTRAL INTERVAL
C             V(5,   H2O    FOR 5 TH SPECTRAL INTERVAL
C             V(6,   CO2    FOR 2 ND SPECTRAL INTERVAL
C             V(7,   CO2    FOR 3 RD SPECTRAL INTERVAL
C             V(8,   CO2    FOR 4 TH SPECTRAL INTERVAL
C             V(9,   P-CONT FOR 2 ND, 3 RD, 4 TH SPECTRAL INTERVALS 
C             V(10,  E-CONT FOR 1 ST AND 5 TH SPECTRAL INTERVALS
C             V(11,  E-CONT FOR 2 ND, 3 RD, 4 TH SPECTRAL INTERVALS 
C             V(12   O3-M   FOR 4 TH SPECTRAL INTERVAL
C             V(13   O3-MP  FOR 4 TH SPECTRAL INTERVAL
C             V(14,  AEROSOLS IN THE 1 ST INTERVAL
C             V(15,  AEROSOLS IN THE 2 ND INTERVAL
C             V(16,  AEROSOLS IN THE 3 RD INTERVAL
C             V(17,  AEROSOLS IN THE 4 TH INTERVAL
C             V(18,  AEROSOLS IN THE 5 TH INTERVAL
  
C     WRITE(6,888) ((V(1,KI,NK),KI=1,13),NK=1,NGLP1)
C     DO 149 I=IL1,IL2
C         WRITE(6,888) (V(I,KI,1),KI=1,13)
C 149 CONTINUE
  
C.......... CALCULATES PLANCK FUNCTIONS 
C           B(T) AND DB/DT CALCULATED FROM POLYNOMIALS OF 5 TH ORDER
C           THE LONGWAVE SCHEME CONSIDERS 5 SPECTRAL INTERVALS: 
C                    1-   0 - 500 CM-1
C                    2- 500 - 800 CM-1
C                    3- 800 - 970 CM-1  +  1110 - 1250 CM-1 
C                    4- 970 - 1110 CM-1 
C                    5- 1250- 2820 CM-1 (NORMALIZED FOR INFINITY
C                                             TO SIGMA*T**4)
  
      DO 187 INU=1,NINT 
  
C . . . . UPWARD: LEVELS FROM SURFACE TO KMX (TOP)
  
          DO 154 K=1,KMX
              DO 150 I=IL1,IL2
                  TI(I)=(T(I,K)-TSTAND)/TSTAND
                  RES(I)=XP(1,6,INU)
  150         CONTINUE
              DO 152 LL=1,5 
                  L=6-LL
                  DO 151 I=IL1,IL2
                      RES(I)=RES(I)*TI(I)+XP(1,L,INU) 
  151             CONTINUE
  152         CONTINUE
              DO 153 I=IL1,IL2
                  BINT(I,K)=BINT(I,K)+RES(I)
  153         CONTINUE
  154     CONTINUE
  
C . . . . TOP OF THE ATMOSPHERE (KMXP)
  
          DO 160 I=IL1,IL2
              TI(I)=(T(I,KMXP)-TSTAND)/TSTAND 
              RES(I)=XP(1,6,INU)
  160     CONTINUE
          DO 162 LL=1,5 
              L=6-LL
              DO 161 I=IL1,IL2
                  RES(I)=RES(I)*TI(I)+XP(1,L,INU) 
  161         CONTINUE
  162     CONTINUE
          DO 163 I=IL1,IL2
              BINT(I,KMXP)=BINT(I,KMXP)+RES(I)
              BTOP(I,INU)=RES(I)
  163     CONTINUE
  
C . . . . SURFACE: TEMPERATURE DISCONTINUITY DT0
  
          DO 170 I=IL1,IL2
              TI(I)=(T(I,1)+DT0(I)-TSTAND)/TSTAND 
              RES(I)=XP(1,6,INU)
  170     CONTINUE
          DO 172 LL=1,5 
              L=6-LL
              DO 171 I=IL1,IL2
                  RES(I)=RES(I)*TI(I)+XP(1,L,INU) 
  171         CONTINUE
  172     CONTINUE
          DO 173 I=IL1,IL2
              BSOL(I,INU)=RES(I)
              BSOLIN(I)=BSOLIN(I)+BSOL(I,INU) 
  173     CONTINUE
  
C . . . .  DERIVATIVE OF PLANCK FUNCTION WITHIN THE LAYERS
  
          DO 185 K=1,KMX
              KJ=(K-1)*NG1P1+1
              DO 180 I=IL1,IL2
                  TI(I)=(TAVE(I,K)-TSTAND)/TSTAND 
                  RES(I)=XP(2,6,INU)
  180         CONTINUE
              DO 182 LL=1,5 
              L=6-LL
                  DO 181 I=IL1,IL2
                      RES(I)=RES(I)*TI(I)+XP(2,L,INU) 
  181             CONTINUE
  182         CONTINUE
              DO 184 IG=1,NG1 
                  KJPIG=KJ+IG 
                  DO 183 I=IL1,IL2
                      W(I,INU,KJPIG)=RES(I) 
  183             CONTINUE
  184         CONTINUE
  185     CONTINUE
  187 CONTINUE
  
C........ VERTICAL INTEGRATION: 
C         TRAPEZOIDAL RULE EXCEPT FOR THE LAYERS ADJACENT TO THE LEVEL
C         OF COMPUTATION WHERE A NG1-POINTS GAUSS QUADRATURE IS USED
  
      IND1=0
      IND3=0
      IND4=1
      IND2=1
  
C . . . . CALCULATES CONTRIBUTION FROM ADJACENT LAYERS (XADJ#)
  
C . . . . DOWNWARD TRANSMITTANCE (XADJD): 
  
      DO 206 K=1,KMX
          KN=(K-1)*NG1P1+1
          KJ=K
          KJP1=KJ+1 
          KX=(KJ-1)*NG1P1+1 
          DO 201 I=IL1,IL2
              XDG(I)=0. 
              ITX(I)=IXT(I,K) 
              DTS(I)=DST(I,K) 
  201     CONTINUE
          DO 204 IG=1,NG1 
              KD=KX+IG
              DO 2021 IA=1,NUA
                  DO 202 I=IL1,IL2
                      UU(I,IA)=V(I,IA,KN)-V(I,IA,KD)
  202             CONTINUE
 2021         CONTINUE
  
              CALL HORNEL(UU,TT,XN,XD,IND1,ITX,IL1,IL2) 
  
              DO 203 I=IL1,IL2
               WTR=W(I,1,KD)*TT(I,1)        *TT(I,9)          *TT(I,15) 
     1           + W(I,2,KD)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16) 
     2           + W(I,3,KD)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17) 
     3           + W(I,4,KD)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18) 
     4           + W(I,5,KD)*TT(I,5)        *TT(I,13)         *TT(I,19) 
                  XDG(I) = XDG(I)+WTR*XT(I,KD)*WG1(IG)
  203         CONTINUE
  204     CONTINUE
          DO 205 I=IL1,IL2
              DZXDG=DZ(I,KJ)*XDG(I) 
              XADJD(I,K)=DZXDG
              CNTRB(I,K,KJP1)=DZXDG 
  205     CONTINUE
  206 CONTINUE
  
C . . . . UPWARD TRANSMITTANCE (XADJM): 
  
      DO 212 K=2,KMXP 
          KM1=K-1 
          KN=(K-1)*NG1P1+1
          KJ=K-1
          KX=(KJ-1)*NG1P1+1 
          DO 207 I=IL1,IL2
              XMG(I)=0. 
              ITX(I)=IXT(I,KM1) 
              DTS(I)=DST(I,KM1) 
  207     CONTINUE
          DO 210 IG=1,NG1 
              KM=KX+IG
              DO 2081 IA=1,NUA
                  DO 208 I=IL1,IL2
                      UU(I,IA)=V(I,IA,KM)-V(I,IA,KN)
  208             CONTINUE
 2081         CONTINUE
  
              CALL HORNEL(UU,TT,XN,XD,IND3,ITX,IL1,IL2) 
  
              DO 209 I=IL1,IL2
               WTR=W(I,1,KM)*TT(I,1)        *TT(I,9)          *TT(I,15) 
     1           + W(I,2,KM)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16) 
     2           + W(I,3,KM)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17) 
     3           + W(I,4,KM)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18) 
     4           + W(I,5,KM)*TT(I,5)        *TT(I,13)         *TT(I,19) 
                  XMG(I) = XMG(I)+WTR*XT(I,KM)*WG1(IG)
  209         CONTINUE
  210     CONTINUE
          DO 211 I=IL1,IL2
              DZXMG=DZ(I,KJ)*XMG(I) 
              XADJM(I,K)=DZXMG
              CNTRB(I,K,KJ)=DZXMG 
  211     CONTINUE
  212 CONTINUE
  
C . . . . CALCULATES CONTRIBUTION BY DISTANT LAYERS 
  
C . . . . DOWNWARD TRANSMITTANCE (XDESC): 
  
      DO 219 K=1,KMXM1
          KP1=K+1 
          KN=(K-1)*NG1P1+1
          KX=(KP1-1)*NG1P1+1
          KD1=KX
          DO 2130 I=IL1,IL2 
              ITX(I)=IXT(I,KP1) 
              DTS(I)=DST(I,KP1) 
 2130     CONTINUE
          DO 2131 IA=1,NUA
              DO 213 I=IL1,IL2
                  UU(I,IA)=V(I,IA,KN)-V(I,IA,KD1) 
  213         CONTINUE
 2131     CONTINUE
  
          CALL HORNEL (UU,TT,XN,XD,IND1,ITX,IL1,IL2)
  
          DO 2135 IA=1,NTRA 
              DO 2134 I=IL1,IL2 
                  TT1(I,IA)=TT(I,IA)
 2134         CONTINUE
 2135     CONTINUE
          DO 218 KJ=KP1,KMX 
              KJP1=KJ+1 
              KX=(KJ-1)*NG1P1+1 
              KD2=KX+NG1P1
              K2=KD2-1
              DO 2140 I=IL1,IL2 
                  ITX(I)=IXT(I,KJ)
                  DTS(I)=DST(I,KJ)
 2140         CONTINUE
              DO 2141 IA=1,NUA
                  DO 214 I=IL1,IL2
                      UU(I,IA)=V(I,IA,KN)-V(I,IA,KD2) 
  214             CONTINUE
 2141         CONTINUE
  
              CALL HORNEL (UU,TT,XN,XD,IND1,ITX,IL1,IL2)
  
              DO 2145 IA=1,NTRA 
                  DO 2144 I=IL1,IL2 
                      TT2(I,IA)=TT(I,IA)
                      TT(I,IA)=(TT1(I,IA)+TT2(I,IA))*0.5
 2144             CONTINUE
 2145         CONTINUE
              DO 215 I=IL1,IL2
               WW=W(I,1,K2)*TT(I,1)        *TT(I,9)          *TT(I,15)
     1          + W(I,2,K2)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16)
     2          + W(I,3,K2)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17)
     3          + W(I,4,K2)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18)
     4          + W(I,5,K2)*TT(I,5)        *TT(I,13)         *TT(I,19)
                  XDG(I)=WW*XT(I,K2)*2. 
                  DZXDG=DZ(I,KJ)*XDG(I) 
                  XDESC(I,K)=XDESC(I,K)+DZXDG 
                  CNTRB(I,K,KJP1)=DZXDG 
  215         CONTINUE
              DO 217 IA=1,NTRA
                  DO 216 I=IL1,IL2
                      TT1(I,IA)=TT2(I,IA) 
  216             CONTINUE
  217         CONTINUE
  218     CONTINUE
  219 CONTINUE
  
C . . . . UPWARD TRANSMITTANCE (XMONT): 
  
      DO 226 K=3,KMXP 
          KN=(K-1)*NG1P1+1
          KM1=K-1 
          KM2=K-2 
          KJ=KM2
          KX=(KJ-1)*NG1P1+1 
          KU1=KX+NG1P1
          DO 2200 I=IL1,IL2 
              ITX(I)=IXT(I,KM2) 
              DTS(I)=DST(I,KM2) 
 2200     CONTINUE
          DO 2201 IA=1,NUA
              DO 220 I=IL1,IL2
                  UU(I,IA)=V(I,IA,KU1)-V(I,IA,KN) 
  220         CONTINUE
 2201     CONTINUE
  
          CALL HORNEL (UU,TT,XN,XD,IND3,ITX,IL1,IL2)
  
          DO 2205 IA=1,NTRA 
              DO 2204 I=IL1,IL2 
                  TT1(I,IA)=TT(I,IA)
 2204         CONTINUE
 2205     CONTINUE
          DO 225 L=1,KM2
              KJ=KM1-L
              KX=(KJ-1)*NG1P1+1 
              KU2=KX
              K2=KX+1 
              DO 2210 I=IL1,IL2 
                  ITX(I)=IXT(I,KJ)
                  DTS(I)=DST(I,KJ)
 2210         CONTINUE
              DO 2211 IA=1,NUA
                  DO 221 I=IL1,IL2
                      UU(I,IA)=V(I,IA,KU2)-V(I,IA,KN) 
  221             CONTINUE
 2211         CONTINUE
  
              CALL HORNEL (UU,TT,XN,XD,IND3,ITX,IL1,IL2)
  
              DO 2215 IA=1,NTRA 
                  DO 2214 I=IL1,IL2 
                      TT2(I,IA)=TT(I,IA)
                      TT(I,IA)=(TT1(I,IA)+TT2(I,IA))*0.5
 2214             CONTINUE
 2215         CONTINUE
              DO 222 I=IL1,IL2
               WW=W(I,1,K2)*TT(I,1)        *TT(I,9)          *TT(I,15)
     1          + W(I,2,K2)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16)
     2          + W(I,3,K2)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17)
     3          + W(I,4,K2)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18)
     4          + W(I,5,K2)*TT(I,5)        *TT(I,13)         *TT(I,19)
                  XMG(I)=WW*XT(I,K2)*2. 
                  DZXMG=DZ(I,KJ)*XMG(I) 
                  XMONT(I,K)=XMONT(I,K)+DZXMG 
                  CNTRB(I,K,KJ)=DZXMG 
 222          CONTINUE
              DO 224 IA=1,NTRA
                  DO 223 I=IL1,IL2
                      TT1(I,IA)=TT2(I,IA) 
  223             CONTINUE
  224         CONTINUE
  225     CONTINUE
  226 CONTINUE
  
C . . . . EXCHANGE WITH THE UPPER LIMIT OF THE ATMOSPHERE 
  
      DO 230 K=1,KMXP 
          KN=(K-1)*NG1P1+1
          DO 228 IA=1,NUA 
              DO 227 I=IL1,IL2
                  UU(I,IA)=V(I,IA,KN) 
  227         CONTINUE
  228     CONTINUE
  
          CALL HORNEL (UU,TT,XN,XD,IND2,IXTTOP,IL1,IL2) 
  
          DO 229 I=IL1,IL2
            CNTOP=BTOP(I,1)*TT(I,1)        *TT(I,9)          *TT(I,15)
     1          + BTOP(I,2)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16)
     2          + BTOP(I,3)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17)
     3          + BTOP(I,4)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18)
     4          + BTOP(I,5)*TT(I,5)        *TT(I,13)         *TT(I,19)
              TR TOP(I,K) = CNTOP/BINT(I,KMXP)
              FD(I)=CNTOP-BINT(I,K)-XDESC(I,K)-XADJD(I,K) 
              FDES(I,K)=FD(I) 
              F(I,2,K)=FD(I)
  229     CONTINUE
  230 CONTINUE
  
C . . . . EXCHANGE WITH THE LOWER LIMIT OF THE ATMOSPHERE:  
  
      DO 235 K=1,KMXP 
          KN=(K-1)*NG1P1+1
          DO 231 I=IL1,IL2
              BGND(I)=BSOLIN(I)*EM0(I)-(1.-EM0(I))*F(I,2,1)-BINT(I,1) 
  231     CONTINUE
          DO 233 IA=1,NUA 
              DO 232 I=IL1,IL2
                  UU(I,IA)=V(I,IA,1)-V(I,IA,KN) 
  232         CONTINUE
  233     CONTINUE
  
          CALL HORNEL (UU,TT,XN,XD,IND4,IXTSOL,IL1,IL2) 
  
          DO 234 I=IL1,IL2
             CNSOL=BSOL(I,1)*TT(I,1)        *TT(I,9)          *TT(I,15) 
     1           + BSOL(I,2)*TT(I,2)*TT(I,6)*TT(I,10)         *TT(I,16) 
     2           + BSOL(I,3)*TT(I,3)*TT(I,7)*TT(I,11)         *TT(I,17) 
     3           + BSOL(I,4)*TT(I,4)*TT(I,8)*TT(I,12)*TT(I,14)*TT(I,18) 
     4           + BSOL(I,5)*TT(I,5)        *TT(I,13)         *TT(I,19) 
              CNSOL=CNSOL*BGND(I)/BSOLIN(I) 
              FM(I)=CNSOL+BINT(I,K)-XMONT(I,K)-XADJM(I,K) 
              FMON(I,K)=FM(I) 
              F(I,1,K)=FM(I)
  234     CONTINUE
  235 CONTINUE

C
C........ SAVE CLEAR-SKY FLUXES AT THE TOA AND SURFACE (FULL CALC.)
C
      DO 236 I=IL1,IL2
         CLT(I) = -1.*F(I,1,KMXP)
         CLB(I) = -1.*(F(I,1,1) + F(I,2,1))
  236 CONTINUE  

C........ END OF THE CLEAR-SKY COMPUTATIONS ............................
  
      DO 861 K = 1, KMXP
  861 CONTINUE
  
C........ CALCULATES EFFECTIVE UP AND DOWN EMISSIVITIES TO BE USED IN 
C         'PARTIAL' RADIATIVE  COMPUTATIONS FOR CLEAR SKY.
  
      ZDENO = 0.01
      DO 302 K=1,KMX
          KP1 = K+1 
          DO 301 I=IL1,IL2
              DENO 1 = STFN * T(I,K)**4 + F(I,2,KP1)
              DD = MERGE(DENO 1, ZDENO, ABS(DENO 1) .GT. ZDENO) 
              EMD(I,K)=MAX(MIN((-F(I,2,K)+F(I,2,KP1))/DD,1.),-1.) 
C             IF (K .EQ. 1) EMD(I,K) = MAX(EMD(I,K), 0.3) 
  
              DENO 2 = STFN * T(I,KP1)**4 - F(I,1,K)
              DM = MERGE(DENO 2, ZDENO, ABS(DENO 2) .GT. ZDENO) 
              EMM(I,K)=MAX(MIN((-F(I,1,K)+F(I,1,KP1))/DM,1.),-1.) 
  301     CONTINUE
  302 CONTINUE
  
C........ DO PARTIAL CACULATION OF FLUXES AND HR BASED ON SAVED 
C         EMISSIVITIES IN EACH LAYER (EMM AND EMD). 
  
      ELSE
  
      DO 312 K = 1, KMXP
          DO 310 I = IL1, IL2 
              B INT(I,K) = STFN * T(I,K) ** 4 
  310     CONTINUE
  312 CONTINUE
  
      DO 315 I = IL1, IL2 
          F(I,2,KMXP) = 0.
  315 CONTINUE
  
      DO 321 K = KMX, 1, -1 
          KP1 = K + 1 
          DO 320 I = IL1, IL2 
              F(I,2,K) = F(I,2,KP1) * (1. - EMD(I,K)) - 
     1                    EMD(I,K) * B INT(I,K) 
  320     CONTINUE
  321 CONTINUE
  
      DO 325 I = IL1, IL2 
          F(I,1,1) = EM0(I) * STFN * (T(I,1) + DT0(I))** 4 +
     1               (1. - EM0(I)) * F(I,2,1) 
          B SOL IN(I) = STFN * (T(I,1) + DT0(I)) ** 4 
  325 CONTINUE
  
      DO 331 K = 1, KMX 
          KP1 = K + 1 
          DO 330 I = IL1, IL2 
              F(I,1,KP1) = F(I,1,K) * (1. - EMM(I,K)) + 
     1                     EMM(I,K) * B INT(I,KP1)
  330     CONTINUE
  331 CONTINUE
  
C
C........ SAVE CLEAR-SKY FLUXES AT THE TOA AND SURFACE (PARTIAL CALC.)
C
      DO 333 I=IL1,IL2
         CLT(I) = -1.*F(I,1,KMXP)
         CLB(I) = -1.*(F(I,1,1) + F(I,2,1))
  333 CONTINUE

C........ PARTIAL CALCULATIONS FOR CLOUDY CASES (SET F UP AND F DN) 
  
      IF (MAX C .GT. 0) THEN
  
          DO 746 K1 = 1, KMXP 
              DO 744 K2 = 1, KMXP 
                  DO 740 I = IL1, IL2 
                       F UP(I,K2,K1) = F(I,1,K1)
                       F DN(I,K2,K1) = F(I,2,K1)
  740             CONTINUE
  744         CONTINUE
  746     CONTINUE
  
      DO 360 KC = 1, MAXC 
          KC P 1 = KC + 1 
          DO 350 I = IL1, IL2 
               F UP(I,KCP1,KCP1) = B INT(I,KCP1)
               F DN(I,KCP1,KC) = - B INT(I,KC)
  350    CONTINUE 
  
      DO 352 K = KCP1, KMX
          K P 1 = K + 1 
          DO 351 I = IL1, IL2 
              F UP(I,KCP1,KP1) = F UP(I,KCP1,K) * (1.-EMM(I,K)) + 
     1                         EMM(I,K) * B INT(I,KP1)
  351     CONTINUE
  352 CONTINUE
  
C . . . . CORRECTION TERM FOR THE UPWARD COMPONENTS IN CLOUDY SKIES.
      DO 359 K = KCP1, KMX
          K P 1 = K + 1 
          DO 358 I = IL1, IL2 
              X = (F(I,1,KP1) - B INT(I,KP1)) * 0.05
              F UP(I,KCP1,KP1) = F UP(I,KCP1,KP1) - X 
  358     CONTINUE
  359 CONTINUE
  
      IF (KC .GT. 1) THEN 
      DO 354 K = KC, 2, -1
         K M 1 = K - 1
          DO 353 I = IL1, IL2 
              ED = EMD(I,KM1) 
              IF (KM1 .EQ. 1) ED = MAX (EMD(I,KM1), 0.2)
              F DN(I,KCP1,KM1) = F DN(I,KCP1,K) * (1.-ED) - 
     1                         B INT(I,KM1) * ED
  353     CONTINUE
  354 CONTINUE
  
C . . . . CORRECTION TERM FOR THE DOWNWARD COMPONENT IN CLOUDY SKIES. 
      DO 357 K = KC, 2, -1
          K M 1 = K - 1 
          DO 356 I = IL1, IL2 
              X = (F(I,2,KM1) + B INT(I,KM1)) * 0.075 
              F DN(I,KCP1,KM1) = F DN(I,KCP1,KM1) - X 
  356     CONTINUE
  357 CONTINUE
      END IF
  360 CONTINUE
  
      END IF
      END IF
C . . . . END OF BLOCK SPECIFIC TO PARTIAL CALCULATIONS 
  
  
C........ TAKING THE CLOUD COVER INTO ACCOUNT 
  
      IF (MAXC .GT. 0) THEN 
      MXP1=MAXC+1 
      MXM1=MAXC-1 
      IF (LCLW) THEN
  
C . . . . INITIALIZING TO CLEAR-SKY FLUXES
  
          DO 703 K1=1,KMXP
              DO 702 K2=1,KMXP
                  DO 701 I=IL1,IL2
                      FUP(I,K2,K1)=FMON(I,K1) 
                      FDN(I,K2,K1)=FDES(I,K1) 
  701             CONTINUE
  702         CONTINUE
  703     CONTINUE
  
C . . . . CALCULATES FLUX DISTRIBUTIONS FOR A TOTAL CLOUD COVER OF
C         EMISSIVITY UNITY FILLING THE DIFFERENT LAYERS BETWEEN THE 
C         SURFACE AND 'MAXC' (INDEX OF THE HIGHEST CLOUDY LAYER OVER
C         THE 'ILX' ATMOSPHERIC COLUMNS  TAKEN INTO ACCOUNT IN ONE
C         CALCULATION). 
  
          DO 714 KC=1,MAXC
              KCLOUD=KC 
              KCP1=KCLOUD+1 
              DO 708 K=KCP1,KMXP
                  KM1=K-1 
                  DO 704 I=IL1,IL2
                      FM(I)=0.
  704             CONTINUE
                  LO1 = KCP1 .GT. KM1 
                  DO 706 KJ=KCP1,KM1
                      DO 705 I=IL1,IL2
                          CNTRI=MERGE(0., CNTRB(I,K,KJ), LO1) 
                          FM(I)=FM(I)+CNTRI 
  705                 CONTINUE
  706             CONTINUE
                  DO 707 I=IL1,IL2
                      FUP(I,KCP1,K)=BINT(I,K)-FM(I) 
  707             CONTINUE
  708         CONTINUE
              DO 713 K=1,KCLOUD 
                  KP1=K+1 
                  DO 709 I=IL1,IL2
                      FD(I)=0.
  709             CONTINUE
                  DO 711 KJ=KP1,KCLOUD
                      DO 710 I=IL1,IL2
                          CNTRI=MERGE(0., CNTRB(I,K,KJ), LO1) 
                          FD(I)=FD(I)+CNTRI 
  710                 CONTINUE
  711             CONTINUE
                  DO 712 I=IL1,IL2
                      FDN(I,KCP1,K)=-BINT(I,K)-FD(I)
  712             CONTINUE
  713         CONTINUE
  714     CONTINUE
  
      END IF
  
C . . . . CALCULATES FLUX DISTRIBUTIONS FOR PARTIAL AND/OR MULTIPLE 
C         LAYERED CLOUDINESS USING LINEAR COMBINATIONS OF FLUXES
C         PREVIOUSLY DEFINED ASSUMING RANDOM OVERLAPPING OF THE CLOUD 
C         LAYERS. 
  
C . . . . DOWNWARD FLUXES ABOVE HIGHEST CLOUD ARE CLEAR-SKY FLUXES
  
          DO 716 K=MXP1,KMXP
              DO 715 I=IL1,IL2
                  F(I,2,K)=FDN(I,1,K) 
  715         CONTINUE
  716     CONTINUE
          IF(MXM1 .NE. 0) THEN
              DO 721 K=1,MXM1 
                  KP1 = K + 1 
                  KR = KMXP - K + 1 
                  KP1R = KMXP - KP1 + 1 
                  DO 717 I=IL1,IL2
                      FD(I)=FDN(I,KP1,K)*C MTX(I,KR,KP1R) 
  717             CONTINUE
                  DO 719 KJ=K,MXM1
                      KJ1=KJ+1
                      KJ2=KJ+2
                      KJR = KMXP - KJ + 1 
                      KJ1R= KMXP - KJ1 + 1
                      KJ2R= KMXP - KJ2 + 1
                      DO 718 I=IL1,IL2
                          CCLD = C MTX(I,KR,KJ2R) - C MTX(I,KR,KJ1R)
                          FD(I)=FD(I)+FDN(I,KJ2,K)*CCLD 
  718                 CONTINUE
  719             CONTINUE
                  MXP1R = KMXP - MXP1 + 1 
                  DO 720 I=IL1,IL2
                      FD(I) = FD(I) + FDN(I,1,K)*(1.-CMTX(I,KR,MXP1R))
                      F(I,2,K)=FD(I)
  720             CONTINUE
  721         CONTINUE
          END IF
          K=MAXC
          KP1=K+1 
          KR = KMXP - K + 1 
          KP1R = KMXP - KP1 + 1 
  
C . . . . DOWNWARD FLUX BELOW THE HIGHEST CLOUD  LAYER AND UPWARD FLUX
C         AT THE SURFACE
  
          DO 722 I=IL1,IL2
              F(I,2,K) = FDN(I,KP1,K) * C MTX(I,KR,KP1R) +
     1                   FDN(I,1,K) * (1. - C MTX(I,KR,KP1R)) 
              F(I,1,1)=EM0(I)*BSOLIN(I)-(1.-EM0(I))*F(I,2,1)
  722     CONTINUE
  
C . . . . UPWARD FLUXES 
  
          DO 727 K=2,KMXP 
              K1=MIN(K,MXP1) 
              K2=K1-1 
              K2M1=K2-1 
              K1R = KMXP - K1 + 1 
              K2R = KMXP - K2 + 1 
              DO 723 I=IL1,IL2
                  FM(I) = FUP(I,K1,K) * C MTX(I,K2R,K1R)
  723         CONTINUE
              IF (K2M1 .GE. 1) THEN 
                  DO 725 KJ=1,K2M1
                      KIJ=K1-KJ 
                      KIJM1 = KIJ - 1 
                      KIJM1R = KMXP - KIJM1 + 1 
                      KIJR = KMXP - KIJ + 1 
                      DO 724 I=IL1,IL2
                          CCLD = C MTX(I,KIJM1R,K1R) - C MTX(I,KIJR,K1R)
                          FM(I)=FM(I)+FUP(I,KIJ,K)*CCLD 
  724                 CONTINUE
  725             CONTINUE
              END IF
              DO 726 I=IL1,IL2
                  FM(I) = FM(I) + FUP(I,1,K)*(1.-CMTX(I,KMXP,K1R))
                  F(I,1,K)=FM(I)
  726         CONTINUE
  727     CONTINUE
      END IF
  
C     WRITE(6,887) (F(I,2,K),K=1,KMXP)
C     WRITE(6,887) (F(1,2,K),K=1,KMXP)
  
      DO 802 K=1,KMXM1
          KP1=K+1 
          DO 801 I=IL1,IL2
              DFNET=F(I,1,KP1)+F(I,2,KP1)-F(I,1,K)-F(I,2,K) 
C . . . . . . 86400 * 9.81 / 1004.64   (PSOL IN PA , DFNET IN W.M-2)
              RESULT(I,K)= CSEC  *DFNET/(PSOL(I)*DSIG(I,K)) 
  801     CONTINUE
  802 CONTINUE
  
C........ COOLING TO SPACE TERM (ONLY REALISTIC ABOVE 10 MBS).
  
      IF(SIG1.LE.0.001 .AND. LCLW) THEN 
          DO 805 I=IL1,IL2
              RESULT(I,KMXM1) = CSEC * STFN * TAVE(I,KMXM1)**4 *
     1                          (TRTOP(I,KMX) - TRTOP(I,KMXM1))/
     2                          (PSOL(I) * DSIG(I,KMXM1)) 
  805     CONTINUE
      ENDIF 
  
C........ CALCULATES EFFECTIVE DOWNWARD EMISSIVITIES TO BE USED IN
C         PARTIAL RADIATIVE COMPUTATIONS TO EVALUATE THE INCIDENT 
C         LONG WAVE FLUX AT THE SURFACE FOR CLOUDY ATMOSPHERE.
  
C     Z DENO = 0.01 
C     DO 820 K = 1, LAY 
C         K P 1 = K + 1 
C         K R   = KMX - K 
C         DO 815 I = IL1, IL2 
C              DENO 1 = STFN * T(I,K) ** 4 + F(I,2,KP1) 
C              DENO = CVMGT (DENO 1, Z DENO, ABS(DENO 1) .GT. Z DENO) 
C              EMD(I,KL) = MAX(MIN((-F(I,2,K)+F(I,2,KP1)) / DENO, 
C    *                           1.), -1.)
C 815     CONTINUE
C 820 CONTINUE
  
C........ END OF LONG WAVE
  
      RETURN
  
C-----------------------------------------------------------------------
  
 885  FORMAT (1X,14E9.2)
 886  FORMAT(1X,16I5) 
 887  FORMAT(1X,10E12.6)
 888  FORMAT (1X,13E9.2)
 889  FORMAT(6X,11E10.4)
  
      END 
