       SUBROUTINE SHORTW2(C MTX, RMU0, IL1, IL2, LON, IFIXCLD,
     1                    DSIG, PSOL, WV, TAVE, TAUAE,
     2                    DSH, DSC, ALBSUR, FRACT, FSDOWN, FSUP,
     3                    FSNET, HEAT, ALBPLA, FSSUR, QOF, TAU, 
     4                    OMEGA, CG,
     5                    UD, UM, REFZ, RJ, RK, TR, 
     6                    C1I, RL, RMUE, RAY1, RUEF,
     7                    RE1, RE2, TR1, TE2, S, G, 
     8                    R23, REF, RMU, R1, W, TO1,
     9                    GG, RNEB, RMUZ, KK, CNEB, TRA1, 
     A                    TAUAZ, SSCAZ, CGAZ, FU, FD, 
     B                    TOT, U1D, U2D, AKI, SS1, RAYL,
     C                    RAY2, TRA2, FACT, SEC, ALB1,
     D                    ALB2, FSU1, FSU2, FAB1, FAB2, TAUCUM, 
     E                    CORR, TAULST, TAUOLD, GLL0) 
C 
C     * DEC  05/88. - M.LAZARE. EFFECTIVE ABSORBER AMOUNTS (NEAR END) 
C     *                         BOUNDED TO ENSURE POSITIVE, IN ORDER
C     *                         THAT INFREQUENT BLOWUPS DUE TO USING
C     *                         VERY SMALL NEGATIVE RESULTS IN TRFCTN 
C     *                         WILL NO LONGER OCCUR. 
C     * JULY 25/88. - M.LAZARE. BASED ON SHORTW EXCEPT GENERALIZED FOR
C     *                         HYBRID MODEL (ARRAYS DSH,DSC,DSIG). 
C 
C  *** SUNRAY 2 
C 
C    THE 'IMX' VERTICAL COLUMNS WHICH ARE PROCESSED TOGETHER ARE ASSUMED
C   TO RECEIVE THE SAME EXTRATERRESTIRAL SOLAR FLUX 'F' WITH DIFFERENT
C   ZENITH ANGLE 'RMU0' 
C 
C      Y. FOUQUART AND B. BONNEL, 1980, 
C      BEITR.PHYS.ATMOSPH. 53,1, 35-62
C 
C      MODIFIED FOR INCLUSION OF AEROSOL EFFECT BY D. TANRE, AUG. 1982
C 
C      VECTORIZED AND OPTIMIZED BY J.J. MORCRETTE, MARCH 1984 
C      ADAPTED TO 2 SPECTRAL INTERVALS BY 
C       Y. FOUQUART, B.BONNEL AND J.J. MORCRETTE, MAY 1985
C 
C-----------------------------------------------------------------------
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      COMMON /RADSIZE/ KMX, LMX, KMXP, NGLP1, NGL 
  
      COMMON /AERO/ TAUA(2,5),SSC ALB(2,5),CGA(2,5),CAER(5,5) 
C 
      COMMON /SWCONS/ CSEC, CH2O, CCAR, CONCO2, F, RAYLG
C 
C 
      COMMON/PADE/APAD(2,3,7),BPAD(2,3,7),D(2,3),CRAY(2,6),SUN(2) 
  
C........ ARRAYS SHARED WITH LONGW AND RADIATN
  
      DIMENSION DSIG(LMX,KMX),PSOL(LMX) 
      DIMENSION WV(LMX,KMX),TAVE(LMX,KMX),TAUAE(LMX,KMX,5)
  
C........ ARRAYS SHARED WITH RADIATN ONLY 
  
      DIMENSION DSH(LMX,KMX),DSC(LMX,KMX),ALBSUR(LMX,2),FRACT(LMX)
      DIMENSION FSDOWN(LMX,KMXP),FSUP(LMX,KMXP),FSNET(LMX,KMXP),
     A          HEAT(LMX,KMX),ALBPLA(LMX),FSSUR(LMX),QOF(LMX,KMX),
     B          TAU(LMX,2,KMX),OMEGA(LMX,2,KMX),CG(LMX,2,KMX) 
  
C........ VARIABLES COMMON TO SHORTWAVE ROUTINE ONLY (GLL0 IS WORK SPACE
C         FOR THE SUBROUTINE TR FCTN) 
  
      DIMENSION UD(LMX,3,KMXP),UM(LMX,KMXP),REFZ(LMX,2,KMXP), 
     A RJ(LMX,6,KMXP),RK(LMX,6,KMXP),TR(LMX,2,KMXP),
     B C1I(LMX,KMXP),RL(LMX,8),RMUE(LMX,KMXP),RAY1(LMX,KMXP), 
     C RUEF(LMX,8),RE1(LMX),RE2(LMX),TR1(LMX),TE2(LMX),S(LMX),G(LMX), 
     D R23(LMX),REF(LMX),RMU(LMX),R1(LMX),W(LMX),TO1(LMX),
     E GG(LMX),RNEB(LMX),RMUZ(LMX),KK(LMX),CNEB(LMX,KMX),TRA1(LMX,KMXP),
     F TAUAZ(LMX,KMX),SSCAZ(LMX,KMX),CGAZ(LMX,KMX), 
     G FU(LMX,KMXP),FD(LMX,KMXP),TOT(LMX,KMXP), 
     H U1D(LMX),U2D(LMX),AKI(LMX,2),SS1(LMX),RAYL(LMX), 
     I RAY2(LMX,KMXP),TRA2(LMX,KMXP),FACT(LMX),SEC(LMX),
     J ALB1(LMX),ALB2(LMX),FSU1(LMX),FSU2(LMX),FAB1(LMX),FAB2(LMX), 
     K TAUCUM(LMX),CORR(LMX),TAULST(LMX,KMX),TAUOLD(LMX,KMX),GLL0(LMX)
  
C........ I/O ARRAYS
  
      DIMENSION C MTX(LMX,KMXP,KMXP), RMU0(LMX) 
  
      LOGICAL LO1 
C 
C-----------------------------------------------------------------------
C 
C-- OBTAIN COMMON BLOCK DATA. 
  
      CALL PADA 
  
C-- THRESHOLDS FOR CLOUDINESS AND CLOUD OPTICAL THICKNESS.
C 
      ZEPSC =1.E-02 
      ZEPSCT=1.E-03 
  
C 
C 
C         Y  : MAGNIFICATION FACTOR 
C         UD : DOWNWARD RADIATION 
C         UM : UPWARD RADIATION 
C 
C     INTERACTIONS BETWEEN OZONE ABSORPTION AND SCATTERING ARE NEGLECTED
C 
C     ABSORBER AMOUNT IN THE K TH LAYER 
C     H2O : UD(1,K) 
C     CO2 : UD(2,K) 
C 
      DO 1 I=IL1,IL2
         C1I(I,KMXP)=0. 
         UD(I,1,KMXP)=0.
         UD(I,2,KMXP)=0.
         UD(I,3,KMXP)=0.
         FACT(I)=RMU0(I)*F*FRACT(I) 
         RMU(I)=SQRT(1224.*RMU0(I)*RMU0(I)+1.)/35.
         SEC(I)=1./RMU(I) 
 1    CONTINUE
C 
C-----------------------------------------------------------------------
C-- CALCULATES OZONE AMOUNTS FOR DOWNWARD LOOKING PATHS (SEC=1./RMU)
C-----------------------------------------------------------------------
      DO 3 K=1,KMX
         KP1=K+1
         L=KMXP-K 
         LP1=L+1
         DO 2 I=IL1,IL2 
            UD(I,3,L)=UD(I,3,LP1)+QOF(I,L)*SEC(I) 
 2       CONTINUE 
 3    CONTINUE
C 
C-----------------------------------------------------------------------
C  CALCULATES OZONE AMOUNTS FOR UPWARD LOOKING PATHS (SEC=1.66) 
C-----------------------------------------------------------------------
      Y=.6024 
      DO 4 I=IL1,IL2
         UM(I,1)=UD(I,3,1)
         U1D(I)=0.
         U2D(I)=0.
 4    CONTINUE
C 
      DO 6 K=2,KMXP 
         KM1=K-1
         DO 5 I=IL1,IL2 
            UM(I,K)=UM(I,KM1)+QOF(I,KM1)*1.66 
 5       CONTINUE 
 6    CONTINUE
C 
C-----------------------------------------------------------------------
C   CALCULATES AMOUNTS OF WATER VAPOR AND UNIFORMLY MIXED GASES (O2+CO2)
C   FOR A VERTICAL PATH 
C-----------------------------------------------------------------------
      DO 8 K=1,KMX
         KP1=K+1
         DO 7 I=IL1,IL2 
            TZ=TAVE(I,K)
            Z=273.15/TZ 
            WH2O=WV(I,K)
            DSHH= DSH(I,K) * CH2O*WH2O
            DSCC= DSC(I,K) * CCAR 
            UD(I,1,K)= PSOL(I) * DSHH * Z**0.45 
            UD(I,2,K)= PSOL(I) * DSCC * Z**1.375
            U1D(I)=U1D(I)+UD(I,1,K) 
            U2D(I)=U2D(I)+UD(I,2,K) 
 7       CONTINUE 
 8    CONTINUE
C 
C-----------------------------------------------------------------------
C-- CALCULATES THE GREY ABSORPTION COEFFICIENTS THAT WOULD GIVE THE SAME
C   CLEAR-SKY DOWNWARD FLUXES AT THE SURFACE (IN SPECTRAL INTERVAL 2) 
C-----------------------------------------------------------------------
      DO 9 I=IL1,IL2
         U1D(I)=U1D(I)*SEC(I) 
         U2D(I)=U2D(I)*SEC(I) 
 9    CONTINUE
C 
      CALL TR FCTN (2, 1, U1D, R1, IL1, IL2, LON, GLL0) 
      CALL TR FCTN (2, 2, U2D, S ,IL1, IL2, LON, GLL0)
C 
      DO 10 I=IL1,IL2 
         AKI(I,1)=-LOG(R1(I))/U1D(I) 
         AKI(I,2)=-LOG(S(I) )/U2D(I) 
 10   CONTINUE
C 
C***********************************************************************
C                  1.0. FIRST SPECTRAL INTERVAL (0.25-0.68 MICRON)
C                       ------------------------------------------
C 
      INU=1 
C 
C-----------------------------------------------------------------------
C   1.1. C1I(*) TOTAL CLOUDINESS ABOVE LEVEL K ASSUMING A RANDOM
C       OVERLAPPING OF THE CLOUD LAYERS, TAKING INTO ACCOUNT THE ACTUAL 
C       OPTICAL THICKNESS OF THE CLOUD (VIA THE FORWARD SCATTERING PEAK)
C-----------------------------------------------------------------------
C 
      DO 110 I=IL1,IL2
         R23(I)=0.
         C1I(I,KMXP)=0. 
         TAU CUM(I)=0.
         CORR(I)=0. 
         DO 109 K = 1, KMX
             TAU LST(I,K) = 0.
 109     CONTINUE 
 110  CONTINUE
      DO 112 K=1,KMX
         L=KMXP-K 
         K P 1 = K + 1
         DO 111 I=IL1,IL2 
            C NEB 1 = C MTX(I,K,KP1)
C . . . . COMPUTES A DOWNWARD REDUCED CLOUDINESS TO CANCELL RANDOM OVERLAP. 
            C NEB 2 = MERGE(0., (C MTX(I,1,KP1) - C MTX(I,1,K)) /
     1                      MAX(1. - C MTX(I,1,K), 0.001)    ,
     2                         C MTX(I,1,K) .GE. 0.999          ) 
            C NEB(I,L) = (1. - CORR(I)) * C NEB 1 + CORR(I) * C NEB 2 
            C NEB(I,L) = MIN(C NEB(I,L),CNEB1)
            TAU CUM(I) = TAU CUM(I) + TAU(I,INU,L)
            TAULST(I,L) = TAULST(I,L) + TAU(I,INU,L)*(CNEB1-CNEB(I,L))
            FACOR = 1.-OMEGA(I,INU,L)*CG(I,INU,L)*CG(I,INU,L) 
            CORR(I) = 1. - EXP(-FACOR * TAU CUM(I) * SEC(I))
            C1I(I,L) = C MTX (I,1,KP1) * CORR(I)
 111     CONTINUE 
 112  CONTINUE
  
      KMX M 1 = KMX - 1 
      DO 116 K = 1, KMX M 1 
          L = KMXP - K
          DO 115 I = IL1, IL2 
            TAUOLD(I,L) = TAU (I,2,L) 
            IF (CNEB(I,L) .GT. 0.) THEN 
              TAU (I,1,L) = (TAU(I,1,L) * CNEB(I,L) + TAULST(I,L-1)) /
     1                                 CNEB(I,L)
              TAU (I,2,L) = (TAU(I,2,L) * CNEB(I,L) + TAULST(I,L-1)) /
     1                                 CNEB(I,L)
              OMEGA(I,2,L)=1.-(1.-OMEGA(I,2,L))*TAUOLD(I,L)/TAU(I,2,L)
              OMEGA(I,2,L) = MIN (OMEGA(I,2,L), OMEGA(I,1,L)) 
            END IF
  115     CONTINUE
  116 CONTINUE
C 
C-----------------------------------------------------------------------
C--  1.2. OPTICAL PARAMETERS FOR THE MIXING OF AEROSOLS 
C-----------------------------------------------------------------------
      DO 124 K=1,KMX
         DO 120 I=IL1,IL2 
            CGAZ(I,K)=0.
            SSCA Z(I,K)=0.
            TAUAZ(I,K)=0. 
 120     CONTINUE 
         DO 122 IAE=1,5 
            DO 121 I=IL1,IL2
               TAUAZ(I,K)=TAUAZ(I,K)+TAUAE(I,K,IAE)*TAUA(INU,IAE) 
               SSCA Z(I,K)=SSCA Z(I,K)+TAUAE(I,K,IAE) 
     S                   * TAUA(INU,IAE)*SSC ALB(INU,IAE) 
               CGAZ(I,K) = CGAZ(I,K) +TAUAE(I,K,IAE)
     S                   * TAUA(INU,IAE)*SSC ALB(INU,IAE)*CGA(INU,IAE)
 121        CONTINUE
 122     CONTINUE 
         DO 123 I=IL1,IL2 
            CGAZ(I,K)=CGAZ(I,K)/SSCA Z(I,K) 
            SSCA Z(I,K)=SSCA Z(I,K)/TAUAZ(I,K)
 123     CONTINUE 
 124  CONTINUE
C 
C-----------------------------------------------------------------------
C   1.3. TRANSMISSION IN THE FORWARD SCATTERING PEAK
C-----------------------------------------------------------------------
      DO 131 I=IL1,IL2
         SS1(I)=0.
         TOT(I,KMXP)=1. 
 131  CONTINUE
      DO 133 K=1,KMX
         L=KMXP-K 
         DO 132 I=IL1,IL2 
            SS1(I)=SS1(I)+TAU(I,INU,L)*(1.-OMEGA(I,INU,L)*CG(I,INU,L))
            TOT(I,L)=EXP(-SS1(I)*SEC(I))
 132     CONTINUE 
 133  CONTINUE
C 
C-----------------------------------------------------------------------
C   1.4. OPTICAL THICKNESS FOR RAYLEIGH SCATTERING
C-----------------------------------------------------------------------
      IF(IFIXCLD.EQ.0) THEN 
         DO 141 I=IL1,IL2 
            RAYL(I)=CRAY(INU,6) 
 141     CONTINUE 
         DO 143 JR=1,5
            IR=6-JR 
            DO 142 I=IL1,IL2
               RAYL(I)=RAYL(I)*RMU(I)+CRAY(INU,IR)
 142        CONTINUE
 143     CONTINUE 
      ELSE
         DO 145 I=IL1,IL2 
            RAYL(I)=0.06
  145    CONTINUE 
      ENDIF 
C 
C-----------------------------------------------------------------------
C   1.5. REFLECTIVITIES AND TRANSMISSIVITIES ARE CALCULATED WITHOUT 
C        GASEOUS ABSORPTION, I.E., FOR A PURELY SCATTERING ATMOSPHERE 
C        INCLUDING RAYLEIGH, CLOUDS, AEROSOLS 
C 
C        SURFACE CONDITIONS    ALBS 
C-----------------------------------------------------------------------
      DO 150 I=IL1,IL2
         REFZ(I,2,1)=ALBSUR(I,INU)
         REFZ(I,1,1)=ALBSUR(I,INU)
 150  CONTINUE
C 
      Y=1.66
C 
      DO 155 K=2,KMXP 
         KM1=K-1
         DO 151 I=IL1,IL2 
            RNEB(I)=CNEB(I,KM1) 
            RE1(I)=0. 
            TR1(I)=0. 
            RE2(I)=0. 
            TE2(I)=0. 
C 
C-----------------------------------------------------------------------
C   1.5.1.  EQUIVALENT ZENITH ANGLE OBTAINED AS A LINEAR MIX OF RMU, THE
C           ZENITH ANGLE FOR DIRECT RADIATION (CLEAR-SKY+FORWARD PEAK)
C           AND OF THE DIFFUSE ZENITH ANGLE (SEC=1.66) DEPENDING ON THE 
C           OVERHEAD CLOUDINESS 
C-----------------------------------------------------------------------
      XMUE=(1.-C1I(I,K))*SEC(I) 
     S                 +C1I(I,K)*(TOT(I,K)*SEC(I)+(1.-TOT(I,K))*1.66) 
      RMUE(I,K)=1./XMUE 
C 
C     REFLECTIVITY OF LAYER KM1 DUE TO RAYLEIGH SCATTERING
C 
      R=RAYL(I)*DSIG(I,KM1) 
C 
C-- ORIGINAL FORMULA (FOUQUART) 
C     RAY(I,KM1)=0.5*R/(RMUE(I,K)+R)
C 
C-- TANRE'S FORMULA 
C     RAY(I,KM1)= R / (2.*RMUE(I,K) + R)
C 
C-- BURIEZ' FORMULA, JUIN 82
C     R=DSIG(I,KM1) 
C     RAY(I,KM1)= 0.0536*R/(RMUE(I,K)+0.063+0.055*R)
C 
C-----------------------------------------------------------------------
C   1.5.2.  TANRE'S FORMULA MODIFYING THE CONTRIBUTION OF RAYLEIGH
C           SCATTERING TO ACCOUNT FOR EFFECTS OF AEROSOLS 
C-----------------------------------------------------------------------
      Y=0.6024
      PT=R+SSCA Z(I,KM1)*TAUAZ(I,KM1) 
      XMPT=(1.-SSCA Z(I,KM1))*TAUAZ(I,KM1)
      BPT=0.5*(1.-TAUAZ(I,KM1)*CGAZ(I,KM1)/(R+TAUAZ(I,KM1)))*PT 
      TRA1(I,KM1)=1./(1.+XMPT*XMUE+BPT*XMUE)
      RAY1(I,KM1)=TRA1(I,KM1)*BPT*XMUE
      TRA2(I,KM1)=1./(1.+XMPT*Y   +BPT*Y   )
      RAY2(I,KM1)=TRA1(I,KM1)*BPT*Y 
  
C-----------------------------------------------------------------------
C   1.5.3. INTRODUCING THE EFFECT OF A CLOUD LAYER IN THE SCATTERING
C          PROCESS: REFLECTIVITY AND TRANSMISSIVITY FOR A GIVEN LAYER 
C          ARE COMPUTED WITH A DELTA-EDDINGTON APPROXIMATION. 
C-----------------------------------------------------------------------
      K0=0
      K1=1
      LO1=RNEB(I).GT.ZEPSC
C***  KK(I)=CVMGT(K0,K1,LO1) 
      IF(LO1)         THEN
         KK(I)=K0
      ELSE
         KK(I)=K1
      ENDIF
      TAU(I,INU,KM1)=MERGE( MAX(TAU(I,INU,KM1),ZEPSCT) , 0., LO1) 
C     KK(I)=K1
C     IF(RNEB(I).GT.ZEPSC) KK(I)=K0 
C 
      W(I)=OMEGA(I,INU,KM1) 
      TO1(I)=TAU(I,INU,KM1)/W(I)
      REF(I)=REFZ(I,1,KM1)
      GG(I)=CG(I,INU,KM1) 
      RMUZ(I)=RMUE(I,K) 
 151  CONTINUE
C 
      IF(IFIXCLD.LE.0) THEN 
          CALL DELTA E (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON)
      ELSE
C         CALL AFLUXV (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON) 
      ENDIF 
  
C-----------------------------------------------------------------------
C-- REFZ(2,K), TR(2,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A NON-REFLECTING UNDERLYING LAYER  R=0. 
C-- REFZ(1,K), TR(1,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A REFLECTING UNDERLYING LAYER R=0.
C-----------------------------------------------------------------------
      DO 154 I=IL1,IL2
C 
C     REFZ(I,2,K)=(1.-RNEB(I))*(RAY1(I,KM1) ) 
C    1                                                 +  RNEB(I)*RE1(I)
C 
C     TR(I,2,KM1)=RNEB(I)*TR1(I) +
C    1  (TRA1(I,KM1) )                                 * (1.-RNEB(I)) 
C 
      REFZ(I,1,K)=(1.-RNEB(I))*(RAY1(I,KM1)+REFZ(I,1,KM1)*TRA1(I,KM1)*
     1   TRA2(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1) ))  +  RNEB(I)*RE2(I)
C 
      TR(I,1,KM1)=RNEB(I)*TE2(I) +
     1  (TRA1(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1) )) * (1.-RNEB(I))
C 
  
  
 154  CONTINUE
 155  CONTINUE
C 
C 
C-----------------------------------------------------------------------
C   1.6.    REFLECTIVITY BETWEEN SURFACE AND LEVEL K (I.E., RK) AND 
C           TRANSMISSIVITY BETWEEN TOP AND LEVEL K (I.E., RJ) ARE 
C           COMPUTED
C-----------------------------------------------------------------------
      JJ=2
      DO 160 I=IL1,IL2
         RJ(I,JJ,KMXP)=1. 
         RK(I,JJ,KMXP)=REFZ(I, 1,KMXP)
 160  CONTINUE
C 
      DO 162 K=1,KMX
         L=KMXP-K 
         LP1=L+1
         DO 161 I=IL1,IL2 
            RE11=RJ(I,JJ,LP1)*TR(I, 1,L)
            RJ(I,JJ,L)=RE11 
            RK(I,JJ,L)=RE11*REFZ(I, 1,L)
 161     CONTINUE 
 162  CONTINUE
C 
C-----------------------------------------------------------------------
C   1.7.   ABSORPTION BY O3 IS INTRODUCED AND FLUXES ARE DETERMINED FROM
C          RJ (DOWNWARD) AND RK (UPWARD). 
C-----------------------------------------------------------------------
      DO 173 K=1,KMXP 
         L=KMXP-K+1 
         DO 171 I=IL1,IL2 
            W(I)=UD(I,3,L)
 171     CONTINUE 
C 
         CALL TR FCTN (INU, 3, W, R1, IL1, IL2, LON, GLL0)
C 
         DO 172 I=IL1,IL2 
            FD(I,L)= R1(I) * RJ(I,JJ,L) * SUN(INU) * FACT(I)
 172     CONTINUE 
 173  CONTINUE
C 
      DO 174 I=IL1,IL2
         FU(I,1)=ALBSUR(I,INU)*FD(I,1)
 174  CONTINUE
      DO 177 K=1,KMXP 
         DO 175 I=IL1,IL2 
            W(I)=UM(I,K)
 175     CONTINUE 
C 
         CALL TR FCTN (INU, 3, W, R1, IL1, IL2, LON, GLL0)
C 
         DO 176 I=IL1,IL2 
            FU(I,K)= R1(I) * RK(I,JJ,K) * SUN(INU) * FACT(I)
 176     CONTINUE 
 177  CONTINUE
C 
C 
C***********************************************************************
C             2.0. SECOND SPECTRAL INTERVAL  (0.68-4.00 MICRONS)
C                  ---------------------------------------------
C 
      INU=2 
C 
C-----------------------------------------------------------------------
C   2.1. C1I(*) TOTAL CLOUDINESS ABOVE LEVEL K ASSUMING A RANDOM
C       OVERLAPPING OF THE CLOUD LAYERS, TAKING INTO ACCOUNT THE ACTUAL 
C       OPTICAL THICKNESS OF THE CLOUD (VIA THE FORWARD SCATTERING PEAK)
C-----------------------------------------------------------------------
C 
C     DO 210 I=IL1,IL2
C        R23(I)=0.
C        C1I(I,KMXP)=0. 
C210  CONTINUE
C     DO 212 K=1,KMX
C        L=KMXP-K 
C         K P 1 = K + 1 
C        DO 211 I=IL1,IL2 
C           CNEB(I,L) = C MTX(I,K,KP1)
C-- 0.28 = 1.- OMEGA*G*G
C           CORR(I)=1.-EXP(-0.28*TAU(I,INU,L)*SEC(I)) 
C           FACOR=1.-OMEGA(I,INU,L)*CG(I,INU,L)*CG(I,INU,L) 
C           CORR(I)=1.-EXP(-FACOR*TAU(I,INU,L)*SEC(I))
C           R23(I)=1.-(1.-CNEB(I,L)*CORR(I))*(1.-R23(I))
C           C1I(I,L)=R23(I) 
C211     CONTINUE 
C212  CONTINUE
C 
C-----------------------------------------------------------------------
C--  2.2. OPTICAL PARAMETERS FOR THE MIXING OF AEROSOLS 
C-----------------------------------------------------------------------
      DO 224 K=1,KMX
         DO 220 I=IL1,IL2 
            CGAZ(I,K)=0.
            SSCA Z(I,K)=0.
            TAUAZ(I,K)=0. 
 220     CONTINUE 
         DO 222 IAE=1,5 
            DO 221 I=IL1,IL2
               TAUAZ(I,K)=TAUAZ(I,K)+TAUAE(I,K,IAE)*TAUA(INU,IAE) 
               SSCA Z(I,K)=SSCA Z(I,K)+TAUAE(I,K,IAE) 
     S                   * TAUA(INU,IAE)*SSC ALB(INU,IAE) 
               CGAZ(I,K) = CGAZ(I,K) +TAUAE(I,K,IAE)
     S                   * TAUA(INU,IAE)*SSC ALB(INU,IAE)*CGA(INU,IAE)
 221        CONTINUE
 222     CONTINUE 
         DO 223 I=IL1,IL2 
            CGAZ(I,K)=CGAZ(I,K)/SSCA Z(I,K) 
            SSCA Z(I,K)=SSCA Z(I,K)/TAUAZ(I,K)
 223     CONTINUE 
 224  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.3. TRANSMISSION IN THE FORWARD SCATTERING PEAK
C-----------------------------------------------------------------------
      DO 231 I=IL1,IL2
         SS1(I)=0.
         TOT(I,KMXP)=1. 
 231  CONTINUE
      DO 233 K=1,KMX
         L=KMXP-K 
         DO 232 I=IL1,IL2 
            SS1(I)=SS1(I)+TAU(I,INU,L)*(1.-OMEGA(I,INU,L)*CG(I,INU,L))
            TOT(I,L)=EXP(-SS1(I)*SEC(I))
 232     CONTINUE 
 233  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.4. OPTICAL THICKNESS FOR RAYLEIGH SCATTERING
C-----------------------------------------------------------------------
      IF(IFIXCLD.EQ.0) THEN 
         DO 241 I=IL1,IL2 
            RAYL(I)=CRAY(INU,6) 
 241     CONTINUE 
         DO 243 JR=1,5
            IR=6-JR 
            DO 242 I=IL1,IL2
               RAYL(I)=RAYL(I)*RMU(I)+CRAY(INU,IR)
 242        CONTINUE
 243     CONTINUE 
      ELSE
         DO 245 I=IL1,IL2 
            RAYL(I)=0.06
  245    CONTINUE 
      ENDIF 
C 
C-----------------------------------------------------------------------
C   2.5. REFLECTIVITIES AND TRANSMISSIVITIES ARE FIRST CALCULATED 
C        WITHOUT GASEOUS ABSORPTION, I.E., FOR A PURELY SCATTERING
C        ATMOSPHERE INCLUDING RAYLEIGH, CLOUDS, AEROSOLS
C 
C        SURFACE CONDITIONS    ALBS 
C-----------------------------------------------------------------------
      DO 250 I=IL1,IL2
         REFZ(I,2,1)=ALBSUR(I,INU)
         REFZ(I,1,1)=ALBSUR(I,INU)
 250  CONTINUE
C 
      Y=1.66
C 
      DO 255 K=2,KMXP 
         KM1=K-1
         DO 251 I=IL1,IL2 
            RNEB(I)=CNEB(I,KM1) 
            RE1(I)=0. 
            TR1(I)=0. 
            RE2(I)=0. 
            TE2(I)=0. 
C 
C-----------------------------------------------------------------------
C   2.5.1.  EQUIVALENT ZENITH ANGLE OBTAINED AS A LINEAR MIX OF RMU, THE
C           ZENITH ANGLE FOR DIRECT RADIATION (CLEAR-SKY+FORWARD PEAK)
C           AND OF THE DIFFUSE ZENITH ANGLE (SEC=1.66) DEPENDING ON THE 
C           OVERHEAD CLOUDINESS 
C-----------------------------------------------------------------------
      XMUE=(1.-C1I(I,K))*SEC(I) 
     S                 +C1I(I,K)*(TOT(I,K)*SEC(I)+(1.-TOT(I,K))*1.66) 
      RMUE(I,K)=1./XMUE 
C 
C     REFLECTIVITY OF LAYER KM1 DUE TO RAYLEIGH SCATTERING
C 
      R=RAYL(I)*DSIG(I,KM1) 
C 
C-----------------------------------------------------------------------
C   2.5.2.  TANRE'S FORMULA MODIFYING THE CONTRIBUTION OF RAYLEIGH
C           SCATTERING TO ACCOUNT FOR EFFECTS OF AEROSOLS 
C-----------------------------------------------------------------------
      Y=0.6024
      PT=R+SSCA Z(I,KM1)*TAUAZ(I,KM1) 
      XMPT=(1.-SSCA Z(I,KM1))*TAUAZ(I,KM1)
      BPT=0.5*(1.-TAUAZ(I,KM1)*CGAZ(I,KM1)/(R+TAUAZ(I,KM1)))*PT 
      TRA1(I,KM1)=1./(1.+XMPT*XMUE+BPT*XMUE)
      RAY1(I,KM1)=TRA1(I,KM1)*BPT*XMUE
      TRA2(I,KM1)=1./(1.+XMPT*Y   +BPT*Y   )
      RAY2(I,KM1)=TRA1(I,KM1)*BPT*Y 
C 
C-----------------------------------------------------------------------
C   2.5.3.  INTRODUCING THE EFFECT OF A CLOUD LAYER IN THE SCATTERING 
C           PROCESS: REFLECTIVITY AND TRANSMISSIVITY FOR A GIVEN LAYER
C           ARE COMPUTED WITH A DELTA-EDDINGTON APPROXIMATION.
C-----------------------------------------------------------------------
      K0=0
      K1=1
      LO1=RNEB(I).GT.ZEPSC
C***  KK(I)=CVMGT(K0,K1,LO1) 
      IF(LO1)         THEN
         KK(I)=K0
      ELSE
         KK(I)=K1
      ENDIF
      TAU(I,INU,KM1)=MERGE( MAX(TAU(I,INU,KM1),ZEPSCT) , 0., LO1) 
C     KK(I)=K1
C     IF(RNEB(I).GT.ZEPSC) KK(I)=K0 
C 
      W(I)=OMEGA(I,INU,KM1) 
      TO1(I)=TAU(I,INU,KM1)/W(I)
      REF(I)=REFZ(I,1,KM1)
      GG(I)=CG(I,INU,KM1) 
      RMUZ(I)=RMUE(I,K) 
 251  CONTINUE
C 
      IF(IFIXCLD.LE.0) THEN 
          CALL DELTA E (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON)
      ELSE
C         CALL AFLUXV (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON) 
      ENDIF 
C 
C-----------------------------------------------------------------------
C-- REFZ(2,K), TR(2,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A NON-REFLECTING UNDERLYING LAYER  R=0. 
C-- REFZ(1,K), TR(1,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A REFLECTING UNDERLYING LAYER R=0.
C-----------------------------------------------------------------------
      DO 254 I=IL1,IL2
C 
      REFZ(I,2,K) = (1.-RNEB(I)) * (RAY1(I,KM1) + REFZ(I,2,KM1) * 
     1              TRA1(I,KM1) *  TRA2(I,KM1)) +  RNEB(I)*RE1(I) 
C 
      TR(I,2,KM1) = RNEB(I) * TR1(I) + TRA1(I,KM1) * (1. - RNEB(I)) 
C 
      REFZ(I,1,K)=(1.-RNEB(I))*(RAY1(I,KM1)+REFZ(I,1,KM1)*TRA1(I,KM1)*
     1      TRA2(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1))) + RNEB(I)*RE2(I)
C 
      TR(I,1,KM1) = RNEB(I) * TE2(I) +
     1  (TRA1(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1))) * (1.-RNEB(I)) 
C 
  
 254  CONTINUE
 255  CONTINUE
C 
C 
C-----------------------------------------------------------------------
C   2.6.    REFLECTIVITY BETWEEN SURFACE AND LEVEL K (I.E., RK) AND 
C           TRANSMISSIVITY BETWEEN TOP AND LEVEL K (I.E., RJ) ARE 
C           COMPUTED
C-----------------------------------------------------------------------
      DO 263 JJ=1,2 
         DO 260 I=IL1,IL2 
            RJ(I,JJ,KMXP)=1.
            RK(I,JJ,KMXP)=REFZ(I,JJ,KMXP) 
 260     CONTINUE 
C 
         DO 262 K=1,KMX 
            L=KMXP-K
            LP1=L+1 
            DO 261 I=IL1,IL2
               RE11=RJ(I,JJ,LP1)*TR(I,JJ,L) 
               RJ(I,JJ,L)=RE11
               RK(I,JJ,L)=RE11*REFZ(I,JJ,L) 
 261        CONTINUE
 262     CONTINUE 
 263  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.7.  REFLECTIVITIES AND TRANSMISSIVITIES ARE NOW CALCULATED WITH A 
C        SIMULATED GASEOUS ABSORPTION USING AN EMPIRICAL GREY ABSORPTION
C        COEFFICIENTS  AKI(IABS)
C 
C        IABS=1  WATER VAPOR  ; IABS=2  UNIFORMLY MIXED GASES 
C-----------------------------------------------------------------------
      N=2 
C 
      DO 278 IABS=1,2 
      DO 270 I=IL1,IL2
         REFZ(I,2,1)=ALBSUR(I,INU)
         REFZ(I,1,1)=ALBSUR(I,INU)
 270  CONTINUE
C 
      DO 273 K=2,KMXP 
         KM1=K-1
         DO 271 I=IL1,IL2 
            RNEB(I)=CNEB(I,KM1) 
            AA=UD(I,IABS,KM1) 
            RKI=AKI(I,IABS) 
            S(I)=EXP(-RKI*AA*1.66)
            G(I)=EXP(-RKI*AA/RMUE(I,K)) 
            TR1(I)=0. 
            RE1(I)=0. 
            TE2(I)=0. 
            RE2(I)=0. 
C 
C-----------------------------------------------------------------------
C   2.7.1. INTRODUCING THE EFFECT OF A CLOUD LAYER IN THE SCATTERING
C          PROCESS: REFLECTIVITY AND TRANSMISSIVITY FOR A GIVEN LAYER 
C          ARE COMPUTED WITH A DELTA-EDDINGTON APPROXIMATION. 
C-----------------------------------------------------------------------
           LO1=RNEB(I).GT.ZEPSC 
C***       KK(I)=CVMGT(K0,K1,LO1) 
           IF(LO1)         THEN
              KK(I)=K0
           ELSE
              KK(I)=K1
           ENDIF   
           TAU(I,INU,KM1)=MERGE( MAX(TAU(I,INU,KM1),ZEPSCT) , 0., LO1)
C 
C     KK(I)=K1
C     IF(RNEB(I).GT.ZEPSC) KK(I)=K0 
C 
           TO1(I)=RKI*AA+TAU(I,INU,KM1)/OMEGA(I,INU,KM1)
           W(I)=TAU(I,INU,KM1)/TO1(I) 
           REF(I)=REFZ(I,1,KM1) 
           GG(I)=CG(I,INU,KM1)
           RMUZ(I)=RMUE(I,K)
 271    CONTINUE
C 
      IF(IFIXCLD.LE.0) THEN 
          CALL DELTA E (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON)
      ELSE
C         CALL AFLUXV (REF,GG,W,TO1,RE1,TR1,RMUZ,RE2,TE2,IL1,IL2,LON) 
      ENDIF 
C 
C-----------------------------------------------------------------------
C-- REFZ(2,K), TR(2,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A NON-REFLECTING UNDERLYING LAYER  R=0. 
C-- REFZ(1,K), TR(1,KM1): REFLECTIVITY AND TRANSMISSIVITY OF LAYER KM1
C                         FOR A REFLECTING UNDERLYING LAYER R=0.
C-----------------------------------------------------------------------
      DO 272 I=IL1,IL2
C 
      REFZ(I,2,K)=(1.-RNEB(I))*(RAY1(I,KM1) + REFZ(I,2,KM1)*TRA1(I,KM1)*
     1                       TRA2(I,KM1) * G(I) * S(I)) + RNEB(I)*RE1(I)
C 
      TR(I,2,KM1)=RNEB(I)*TR1(I) +
     1 (TRA1(I,KM1) )                                *G(I)*(1.-RNEB(I)) 
C 
      REFZ(I,1,K)=(1.-RNEB(I))*(RAY1(I,KM1)+REFZ(I,1,KM1)*TRA1(I,KM1)*
     1  TRA2(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1)))*G(I)*S(I) 
     2                                                 +  RNEB(I)*RE2(I)
C 
      TR(I,1,KM1)=RNEB(I)*TE2(I) +
     1  (TRA1(I,KM1)/(1.-RAY2(I,KM1)*REFZ(I,1,KM1) )) *G(I)*(1.-RNEB(I))
C 
  
 272  CONTINUE
 273  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.7.2   REFLECTIVITY BETWEEN SURFACE AND LEVEL K (I.E., RK) AND 
C           TRANSMISSIVITY BETWEEN TOP AND LEVEL K (I.E., RJ) ARE 
C           COMPUTED
C-----------------------------------------------------------------------
      DO 277 KREF=1,2 
         N=N+1
         DO 274 I=IL1,IL2 
            RJ(I,N,KMXP)=1. 
            RK(I,N,KMXP)=REFZ(I,KREF,KMXP)
 274     CONTINUE 
C 
         DO 276 K=1,KMX 
            L=KMXP-K
            LP1=L+1 
            DO 275 I=IL1,IL2
               RE11=RJ(I,N,LP1)*TR(I,KREF,L)
               RJ(I,N,L)=RE11 
               RK(I,N,L)=RE11*REFZ(I,KREF,L)
 275        CONTINUE
 276     CONTINUE 
 277  CONTINUE
 278  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.8.   UPWARD (RK) AND DOWNWARD (RJ) FLUXES ARE NOW CALCULATED
C           WITHOUT GASEOUS ABSORPTION     : N= 1 , 2 
C           WITH GASEOUS ABSORPTION BY H2O : N= 3 , 4 
C           WITH GASEOUS ABSORPTION BY CO2 : N= 5 , 6 
C           N= 1, 3, 5 : NON-REFLECTING UNDERLYING LAYER
C           N= 2, 4, 6 : REFLECTING UNDERLYING LAYER
C-----------------------------------------------------------------------
C-NB:  EE IS AN ADJUSTABLE THRESHOLD DEPENDING ON THE COMPUTER PRECISION
C       ADDED TO 'RJ' TO PREVENT LOGARITHM OF ZERO
C-----------------------------------------------------------------------
      EE=1.E-10 
      DO 283 K=1,KMXP 
         DO 282 JJ=1,5,2
            JJP= JJ+1 
            DO 281 I=IL1,IL2
               RJ(I,JJ,K)=        RJ(I,JJ,K) - RJ(I,JJP ,K) 
               RK(I,JJ,K)=        RK(I,JJ,K) - RK(I,JJP ,K) 
               RJ(I,JJ,K)= MAX( RJ(I,JJ,K) , EE ) 
               RK(I,JJ,K)= MAX( RK(I,JJ,K) , EE ) 
 281        CONTINUE
 282     CONTINUE 
 283  CONTINUE
C 
      DO 286 K=1,KMXP 
         DO 285 JJ=2,6,2
            DO 284 I=IL1,IL2
C              RJ(I,JJ,K)= EE + RJ(I,JJ,K)
C              RK(I,JJ,K)= EE + RK(I,JJ,K)
               RJ(I,JJ,K)= MAX( RJ(I,JJ,K) , EE ) 
               RK(I,JJ,K)= MAX( RK(I,JJ,K) , EE ) 
 284        CONTINUE
 285     CONTINUE 
  
  
 286  CONTINUE
C 
C-----------------------------------------------------------------------
C   2.9. COMPUTES EFFECTIVE ABSORBER AMOUNTS (INVERSE LAPLACE TRANSFORM)
C           UEFF= - LN (F(K)/F(0))/K
C        COMPUTES CORRESPONDING TRANSMISSION FUNCTIONS OF H2O AND CO2 
C        COMPUTES UPWARD AND DOWNWARD FLUXES
C-----------------------------------------------------------------------
C 
      DO 296 K=1,KMXP 
         KKI=1
C 
         DO 294 JJ=1,2
C 
            DO 293 N=1,2
               N2J=N+2*JJ 
               KKP4=KKI+4 
C 
C  EFFECTIVE ABSORBER AMOUNT
C---------------------------
               DO 290 I=IL1,IL2 
                  W(I)=LOG(RJ(I,N,K)/RJ(I,N2J,K))/AKI(I, JJ) 
                  W(I)=MAX(W(I),0.) 
 290           CONTINUE 
C 
C  TRANSMISSION FUNCTION
C-----------------------
               CALL TR FCTN (INU, JJ, W, R1, IL1, IL2, LON, GLL0) 
C 
               DO 291 I=IL1,IL2 
                  RL(I,KKI)=R1(I) 
                  RUEF(I,KKI)=W(I)
                  W(I)=LOG(RK(I,N,K)/RK(I,N2J,K))/AKI(I, JJ) 
                  W(I)=MAX(W(I),0.) 
 291           CONTINUE 
C 
               CALL TR FCTN (INU, JJ, W, R1, IL1, IL2, LON, GLL0) 
C 
               DO 292 I=IL1,IL2 
                  RL(I,KKP4)=R1(I)
                  RUEF(I,KKP4)=W(I) 
 292           CONTINUE 
C 
               KKI=KKI+1
 293        CONTINUE
 294     CONTINUE 
C 
C UPWARD AND DOWNWARD FLUXES WITH H2O AND CO2 ABSORPTION
C-------------------------------------------------------
         DO 295 I=IL1,IL2 
            FSDOWN(I,K)=RJ(I,1,K)*RL(I,1)*RL(I,3) 
     1                +RJ(I,2,K)*RL(I,2)*RL(I,4)
            FSUP(I,K)  =RK(I,1,K)*RL(I,5)*RL(I,7) 
     1                +RK(I,2,K)*RL(I,6)*RL(I,8)
 295     CONTINUE 
 296  CONTINUE
  
  
C 
C-----------------------------------------------------------------------
C   3.0. FLUXES ARE MULTIPLIED BY  OZONE TRANSMISSION FUNCTIONS 
C-----------------------------------------------------------------------
      IABS=3
      DO 304 K=1,KMXP 
         DO 301 I=IL1,IL2 
            W(I)=UD(I,3,K)
 301     CONTINUE 
C 
         CALL TR FCTN (INU, IABS, W, R1, IL1, IL2, LON, GLL0) 
C 
         DO 302 I=IL1,IL2 
            FSDOWN(I,K)=R1(I) * FSDOWN(I,K) * SUN(INU) * FACT(I)
            W(I)=UM(I,K)
 302     CONTINUE 
C 
         CALL TR FCTN (INU, IABS, W, R1, IL1, IL2, LON, GLL0) 
C 
         DO 303 I=IL1,IL2 
            FSUP(I,K) = R1(I) * FSUP(I,K) * SUN(INU) * FACT(I)
 303     CONTINUE 
 304  CONTINUE
  
C 
C-----------------------------------------------------------------------
C   3.1. COMPUTES PARTIAL ALBEDO AND NET FLUX AT SURFACE
C-----------------------------------------------------------------------
      DO 310 I=IL1,IL2
         FD(I,KMXP)=MAX( FD(I,KMXP), EE)
         FSDOWN(I,KMXP)=MAX( FSDOWN(I,KMXP), EE)
         ALB1(I)=MIN( FU(I,KMXP)/FD(I,KMXP), 1.)
         ALB2(I)=MIN( FSUP(I,KMXP)/FSDOWN(I,KMXP), 1.)
         FSU1(I)=FD(I,1)-FU(I,1)
         FSU2(I)=FSDOWN(I,1)-FSUP(I,1)
         FAB1(I)=FD(I,1)-FU(I,1)
         FAB2(I)=FSDOWN(I,1)-FSUP(I,1)
         ALBPLA(I)=(FU(I,KMXP)+FSUP(I,KMXP))/(FD(I,KMXP)+FSDOWN(I,KMXP))
         FSSUR(I)=(FD(I,1)+FSDOWN(I,1))-(FU(I,1)+FSUP(I,1)) 
 310  CONTINUE
C-----------------------------------------------------------------------
C   3.2.    NET FLUXES AND HEATING RATES (DEG.K / DAY)
C-----------------------------------------------------------------------
C 
      DO 321 K=1,KMXP 
         DO 320 I=IL1,IL2 
            FSUP(I,K)=FSUP(I,K)+FU(I,K) 
            FSDOWN(I,K)=FSDOWN(I,K)+FD(I,K) 
            FSNET(I,K)=FSUP(I,K)-FSDOWN(I,K)
 320     CONTINUE 
 321  CONTINUE
  
C 
      DO 323 K=1,KMX
         L=KMXP-K 
         LP1=L+1
         DO 322 I=IL1,IL2 
            HEAT(I,L)=CSEC*(FSNET(I,L)-FSNET(I,LP1))/(PSOL(I)*DSIG(I,L))
 322     CONTINUE 
 323  CONTINUE
  99  RETURN
      END 
