SUBROUTINE swr93(fswc,hswc,ufswc,dfswc,fswl,hswl,ufswl,dfswl,press,coszro, &
     taudar,rh2o,rrco2,ssolar,qo3,nclds,ktopsw,kbtmsw,camt,crr, &
     ctt, alvb,alnb,alvd,alnd,gdfvb,gdfnb,gdfvd,gdfnd)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !SW- RADIATION CODE
  !    INPUTS:PRESS,COSZRO,TAUDAR,RH2O,RRCO2,SSOLAR,QO3,NCLDS,
  !           KTOPSW,KBTMSW,CIRAB,CIRRF,CUVRF,CAMT,
  !           ALVB,ALVD,ALNB,ALND;
  !    OUTPUT:FSWC,HSWC,UFSWC,DFSWC,FSWL,HSWL,UFSWL,DFSWL,
  !           GDFVB,GDFVD,GDFNB,GDFND.
  ! SWR91SIB --- MODIFIED BY K. CAMPANA..06 MAR 92
  ! SWR91SIB --- MODIFIED BY Y. HOU         FEB 93

  IMPLICIT NONE

  INCLUDE "parmeta.f90"
  INCLUDE "hcon.f90"
  INCLUDE "rdparm.f90"
  ! PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE:
  !   IMAX - NO. POINTS ALONG THE LAT. CIRCLE USED IN CALCS.
  !   L    - NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL
  !NOTE: THE USER NORMALLY WILL MODIFY ONLY THE IMAX AND L PARAMETERS
  !      NB IS A SHORTWAVE PARAMETER; OTHER QUANTITIES ARE DERIVED
  !      FROM THE ABOVE PARAMETERS.

  REAL,DIMENSION(imax,lp1),INTENT(OUT) :: fswc,hswc,ufswc,dfswc,fswl,hswl,ufswl
  REAL,DIMENSION(imax,lp1),INTENT(OUT) :: dfswl
  REAL,DIMENSION(imax,lp1),INTENT(IN)  :: press
  REAL,DIMENSION(imax),INTENT(IN)      :: coszro,taudar
  REAL,DIMENSION(imax,l),INTENT(IN)    :: rh2o,qo3
  REAL, INTENT(IN OUT)                     :: rrco2
  REAL, INTENT(IN)                         :: ssolar
  REAL, DIMENSION(nb)          :: abcff,pwts
  REAL, DIMENSION(imax)        :: alvb,alnb,alvd,alnd,gdfvb,gdfnb
  REAL, DIMENSION(imax)        :: gdfvd,gdfnd,rray,secz,refl,temp1,refl2,ccmax
  REAL, DIMENSION(imax,nb)     :: dfntop
  REAL, DIMENSION(imax,lp1)    :: camt,pp,dp,pr2,du,duco2,duo3,ff,ffco2
  REAL, DIMENSION(imax,lp1)    :: ffo3,xamt,ud,ur,udco2,urco2,udo3,uro3,tdco2
  REAL, DIMENSION(imax,lp1)    :: tuco2,tdo3,tuo3,dfn,ufn,cr,ttd,ttu,ct,pptop
  REAL, DIMENSION(imax,lp1)    :: dpcld,ttub1,tucl1,ttdb1,tdcl1,tdcl2,ufntrn
  REAL, DIMENSION(imax,lp1)    :: ufnclu,tclu,dfntrn,dfnclu,tcld,alfa,alfau
  REAL, DIMENSION(imax,llp2)   :: uco2,uo3,tco2,to3
  REAL, DIMENSION(imax,nb,lp1) :: crr,ctt   

  INTEGER, DIMENSION(imax)     :: nclds 
  INTEGER, DIMENSION(imax,lp1) :: ktopsw,kbtmsw

  REAL :: cfco2,cfo3,reflo3,rrayav,tempf,tempg,htemp
  INTEGER :: i,k,n,ip,KK,j1,j2,j3,jtop,nnclds,kclds

  EQUIVALENCE (udo3,uo3(1,1),dfnclu),(uro3,uo3(1,lp2),ufnclu),     &
       (udco2,uco2(1,1),tcld),(urco2,uco2(1,lp2),tclu),      &
       (tdo3,to3(1,1),dfntrn),(tuo3,to3(1,lp2),ufntrn),      &
       (tdco2,tco2(1,1)),(tuco2,tco2(1,lp2)),(ff,alfa ),     &
       (ffco2,alfau ),(ffo3,ttdb1),(du,ttub1),(duco2,tucl1), &
       (duo3,tdcl1),(pr2,tdcl2)

  COMMON /swrsav/ abcff,pwts,cfco2,cfo3,reflo3,rrayav
  !                           D A T A
  !    1  ABCFF / 2*4.0E-5, 0.002, 0.035, 0.377, 1.95, 9.40, 44.6,
  !    1          190.0,    989.0, 2706.0, 39011.0 /
  !    2, PWTS  / 0.5000, 0.121416, 0.0698, 0.1558, 0.0631, 0.0362,
  !    2          0.0243, 0.0158, 0.0087, 0.001467, 0.002342, 0.001075 /
  !    3, CFCO2, CFO3, REFLO3, RRAYAV / 508.96, 466.64, 1.9, 0.144 /
  !    1  ABCFF / 2*4.0E-5, .002, .035, .377, 1.95, 9.40, 44.6, 190. /
  !    2, PWTS  /.5000,.1470,.698,.1443,.0584,.0335,.0225,.0158,.0087/
  !    3, CFCO2, CFO3, REFLO3, RRAYAV / 508.96, 466.64, 1.9, 0.144 /

  !CALCULATE SECANT OF ZENITH ANGLE (SECZ),FLUX PRESSURES(PP),
  !LAYER WIDTH (DP) AND PRESSURE SCALING FACTOR (PR2).
  DO  i=1,imax
     secz(i) = h35e1/SQRT(h1224e3*coszro(i)*coszro(i)+one)
     pp(i,1)   = zero
     pp(i,lp1) = press(i,lp1)
     temp1(i)  = one/press(i,lp1)
  END DO
  DO  i=1,imax*lm1
     pp(i,2) = haf*(press(i,2)+press(i,1))
  END DO
  DO  i=1,imax*l
     dp (i,1) = pp(i,2)-pp(i,1)
     pr2(i,1) = haf*(pp(i,1)+pp(i,2))
  END DO
  DO  k=1,l
     DO  i=1,imax
        pr2(i,k) = pr2(i,k)*temp1(i)
     END DO
  END DO

  !CALCULATE ENTERING FLUX AT THE TOP FOR EACH BAND(IN CGS UNITS)
  DO  n=1,nb
     DO  ip=1,imax
        dfntop(ip,n) = ssolar*h69766e5*coszro(ip)*taudar(ip)*pwts(n)
     END DO
  END DO

  !EXECUTE THE LACIS-HANSEN REFLECTIVITY PARAMETERIZATION
  !FOR THE VISIBLE BAND
  DO  i=1,imax
     rray(i) = hp219/(one+hp816*coszro(i))
     refl(i) = rray(i) + (one-rray(i))*(one-rrayav)*alvb(i)/ (one-alvd(i)*rrayav)
  END DO
  DO  i=1,imax
     rray(i) = 0.104/(one+4.8*coszro(i))
     refl2(i)= rray(i) + (one-rray(i))*(one-0.093)*alvb(i)/ (one-alvd(i)*0.093)
  END DO

  !CALCULATE PRESSURE-WEIGHTED OPTICAL PATHS FOR EACH LAYER
  !IN UNITS OF CM-ATM. PRESSURE WEIGHTING IS USING PR2.
  !DU= VALUE FOR H2O;DUCO2 FOR CO2;DUO3 FOR O3.
  DO  i=1,imax*l
     du   (i,1) = ginv*rh2o(i,1)*dp(i,1)*pr2(i,1)
     duco2(i,1) = (rrco2*ginv*cfco2)*dp(i,1)*pr2(i,1)
     duo3 (i,1) = (ginv*cfo3)*qo3(i,1)*dp(i,1)
  END DO

  !CALCULATE CLEAR SKY SW FLUX
  !OBTAIN THE OPTICAL PATH FROM THE TOP OF THE ATMOSPHERE TO THE
  !FLUX PRESSURE. ANGULAR FACTORS ARE NOW INCLUDED. UD=DOWNWARD
  !PATH FOR H2O,WIGTH UR THE UPWARD PATH FOR H2O. CORRESPONDING
  !QUANTITIES FOR CO2,O3 ARE UDCO2/URCO2 AND UDO3/URO3.
  DO  ip=1,imax
     ud   (ip,1) = zero
     udco2(ip,1) = zero
     udo3 (ip,1) = zero
  END DO
  DO  k=2,lp1
     DO  i=1,imax
        ud   (i,k) = ud   (i,k-1)+du   (i,k-1)*secz(i)
        udco2(i,k) = udco2(i,k-1)+duco2(i,k-1)*secz(i)
        udo3 (i,k) = udo3 (i,k-1)+duo3 (i,k-1)*secz(i)
     END DO
  END DO
  DO  ip=1,imax
     ur   (ip,lp1) = ud   (ip,lp1)
     urco2(ip,lp1) = udco2(ip,lp1)
     uro3 (ip,lp1) = udo3 (ip,lp1)
  END DO
  DO  k=l,1,-1
     DO  ip=1,imax
        ur   (ip,k) = ur   (ip,k+1)+du   (ip,k)*diffctr
        urco2(ip,k) = urco2(ip,k+1)+duco2(ip,k)*diffctr
        uro3 (ip,k) = uro3 (ip,k+1)+duo3 (ip,k)*o3difctr
     END DO
  END DO

  !CALCULATE CO2 ABSORPTIONS . THEY WILL BE USED IN NEAR INFRARED
  !BANDS.SINCE THE ABSORPTION AMOUNT IS GIVEN (IN THE FORMULA USED
  !BELOW, DERIVED FROM SASAMORI) IN TERMS OF THE TOTAL SOLAR FLUX,
  !AND THE ABSORPTION IS ONLY INCLUDED IN THE NEAR IR (50 PERCENT
  !OF THE SOLAR SPECTRUM), THE ABSORPTIONS ARE MULTIPLIED BY 2.
  !SINCE CODE ACTUALLY REQUIRES TRANSMISSIONS, THESE ARE THE
  !VALUES ACTUALLY STORED IN TCO2.
  DO  i=1,imax*ll
     tco2(i,2) = one - two*(h235m3 * EXP(hp26*LOG(uco2(i,2)+h129m2)) -h75826m4)
  END DO

  !NOW CALCULATE OZONE ABSORPTIONS. THESE WILL BE USED IN
  !THE VISIBLE BAND.JUST AS IN THE CO2 CASE, SINCE THIS BAND IS
  !50 PERCENT OF THE SOLAR SPECTRUM,THE ABSORPTIONS ARE MULTIPLIED
  !BY 2. THE TRANSMISSIONS ARE STORED IN TO3.
  htemp = h1036e2*h1036e2*h1036e2
  DO  i=1,imax*ll
     to3(i,2) = one - two*uo3(i,2)*  &
          (h1p082*EXP(hmp805*LOG(one+h1386e2*uo3(i,2)))+  &
          h658m2/(one+htemp*uo3(i,2)*uo3(i,2)*uo3(i,2))+  &
          h2118m2/(one+uo3(i,2)*(h42m2+h323m4*uo3(i,2))))
  END DO

  !START FREQUENCY LOOP (ON N) HERE

  !BAND 1 (VISIBLE) INCLUDES O3 AND H2O ABSORPTION
  DO  i=1,imax*l
     ttd(i,2) = EXP(hm1ez*MIN(fifty,abcff(1)*ud(i,2)))
     ttu(i,1) = EXP(hm1ez*MIN(fifty,abcff(1)*ur(i,1)))
     dfn(i,2) = ttd(i,2)*tdo3(i,2)
     ufn(i,1) = ttu(i,1)*tuo3(i,1)
  END DO
  DO  i=1,imax
     dfn(i,1)   = one
     ufn(i,lp1) = dfn(i,lp1)
  END DO

  !SCALE VISIBLE BAND FLUXES BY SOLAR FLUX AT THE TOP OF THE
  !ATMOSPHERE (DFNTOP(I,1))
  !DFSW/UFSW WILL BE THE FLUXES, SUMMED OVER ALL BANDS
  DO   k=1,lp1
     DO   i=1,imax
        dfswl(i,k) =         dfn(i,k)*dfntop(i,1)
        ufswl(i,k) = refl(i)*ufn(i,k)*dfntop(i,1)
     END DO
  END DO
  DO  i=1,imax
     gdfvb(i) = dfswl(i,lp1)*EXP(-0.15746*secz(i))
     gdfvd(i) = ((one-refl2(i))*dfswl(i,lp1) -  &
          (one-alvb(i)) *gdfvb(i)) / (one-alvd(i))
     gdfnb(i) = zero
     gdfnd(i) = zero
  END DO

  !NOW OBTAIN FLUXES FOR THE NEAR IR BANDS. THE METHODS ARE THE SAME
  !AS FOR THE VISIBLE BAND, EXCEPT THAT THE REFLECTION AND
  !TRANSMISSION COEFFICIENTS (OBTAINED BELOW) ARE DIFFERENT, AS
  !RAYLEIGH SCATTERING NEED NOT BE CONSIDERED.
  DO  n=2,nb
     IF (n == 2) THEN
        !THE WATER VAPOR TRANSMISSION FUNCTION FOR BAND 2 IS EQUAL TO
        !THAT OF BAND 1 (SAVED AS TTD,TTU)
        !BAND 2-9 (NEAR-IR) INCLUDES O3, CO2 AND H2O ABSORPTION
        DO  i=1,imax*l
           dfn(i,2) = ttd(i,2)*tdco2(i,2)
           ufn(i,1) = ttu(i,1)*tuco2(i,1)
        END DO
     ELSE
        !CALCULATE WATER VAPOR TRANSMISSION FUNCTIONS FOR NEAR INFRARED
        !BANDS. INCLUDE CO2 TRANSMISSION (TDCO2/TUCO2), WHICH
        !IS THE SAME FOR ALL INFRARED BANDS.
        DO  i=1,imax*l
           dfn(i,2)=EXP(hm1ez*MIN(fifty,abcff(n)*ud(i,2)))*tdco2(i,2)
           ufn(i,1)=EXP(hm1ez*MIN(fifty,abcff(n)*ur(i,1)))*tuco2(i,1)
        END DO
     END IF
     !AT THIS POINT,INCLUDE DFN(1),UFN(LP1), NOTING THAT DFN(1)=1 FOR
     !ALL BANDS, AND THAT UFN(LP1)=DFN(LP1) FOR ALL BANDS.
     DO  i=1,imax
        dfn(i,1)   = one
        ufn(i,lp1) = dfn(i,lp1)
     END DO
     !SCALE THE PREVIOUSLY COMPUTED FLUXES BY THE FLUX AT THE TOP
     !AND SUM OVER BANDS
     DO  k=1,lp1
        DO  i=1,imax
           dfswl(i,k) = dfswl(i,k) +         dfn(i,k)*dfntop(i,n)
           ufswl(i,k) = ufswl(i,k) + alnb(i)*ufn(i,k)*dfntop(i,n)
        END DO
     END DO
     DO  i=1,imax
        gdfnb(i) = gdfnb(i) + dfn(i,lp1)*dfntop(i,n)
     END DO
  END DO
  DO  i=1,imax*lp1
     fswl(i,1) = ufswl(i,1)-dfswl(i,1)
  END DO
  DO  i=1,imax*l
     hswl(i,1) = radcon*(fswl(i,2)-fswl(i,1))/dp(i,1)
  END DO
  !END OF FREQUENCY LOOP (OVER N)

  !CALCULATE CLOUDY SKY SW FLUX
  kclds=nclds(1)
  DO  i=2,imax
     kclds=MAX(nclds(i),kclds)
  END DO
  DO  i=1,imax*lp1
     dfswc(i,1) = dfswl(i,1)
     ufswc(i,1) = ufswl(i,1)
     fswc (i,1) = fswl (i,1)
  END DO
  DO  i=1,imax*l
     hswc(i,1) = hswl(i,1)
  END DO
  IF (kclds == 0)  RETURN
  DO  i=1,imax*lp1
     xamt(i,1) = camt(i,1)
  END DO
  DO  i=1,imax
     nnclds   = nclds(i)
     ccmax(i) = zero
     IF (nnclds <= 0) CYCLE
     ccmax(i) = one
     DO  k=1,nnclds
        ccmax(i) = ccmax(i) * (one - camt(i,k+1))
     END DO
     ccmax(i) = one - ccmax(i)
     IF (ccmax(i) > zero) THEN
        DO  k=1,nnclds
           xamt(i,k+1) = camt(i,k+1)/ccmax(i)
        END DO
     END IF
  END DO
  DO  i=1,imax*lp1
     ff   (i,1) = diffctr
     ffco2(i,1) = diffctr
     ffo3 (i,1) = o3difctr
  END DO
  DO  ip=1,imax
     jtop = ktopsw(ip,nclds(ip)+1)
     DO  k=1,jtop
        ff   (ip,k) = secz(ip)
        ffco2(ip,k) = secz(ip)
        ffo3 (ip,k) = secz(ip)
     END DO
  END DO
  DO  i=1,imax
     rray(i) = hp219/(one+hp816*coszro(i))
     refl(i) = rray(i) + (one-rray(i))*(one-rrayav)*alvd(i)/ (one-alvd(i)*rrayav)
  END DO
  DO  ip=1,imax
     ud   (ip,1) = zero
     udco2(ip,1) = zero
     udo3 (ip,1) = zero
  END DO
  DO  k=2,lp1
     DO  i=1,imax
        ud   (i,k) = ud   (i,k-1)+du   (i,k-1)*ff   (i,k)
        udco2(i,k) = udco2(i,k-1)+duco2(i,k-1)*ffco2(i,k)
        udo3 (i,k) = udo3 (i,k-1)+duo3 (i,k-1)*ffo3 (i,k)
     END DO
  END DO
  DO  ip=1,imax
     ur   (ip,lp1) = ud   (ip,lp1)
     urco2(ip,lp1) = udco2(ip,lp1)
     uro3 (ip,lp1) = udo3 (ip,lp1)
  END DO
  DO  k=l,1,-1
     DO  ip=1,imax
        ur   (ip,k) = ur   (ip,k+1)+du   (ip,k)*diffctr
        urco2(ip,k) = urco2(ip,k+1)+duco2(ip,k)*diffctr
        uro3 (ip,k) = uro3 (ip,k+1)+duo3 (ip,k)*o3difctr
     END DO
  END DO
  DO  i=1,imax*ll
     tco2(i,2) = one - two*(h235m3 * EXP(hp26*LOG(uco2(i,2)+h129m2)) -h75826m4)
  END DO
  DO  i=1,imax*ll
     to3(i,2) = one - two*uo3(i,2)*  &
          (h1p082*EXP(hmp805*LOG(one+h1386e2*uo3(i,2)))+  &
          h658m2/(one+htemp*uo3(i,2)*uo3(i,2)*uo3(i,2))+  &
          h2118m2/(one+uo3(i,2)*(h42m2+h323m4*uo3(i,2))))
  END DO

  !THE FIRST CLOUD IS THE GROUND; ITS PROPERTIES ARE GIVEN
  !BY REFL (THE TRANSMISSION (0) IS IRRELEVANT FOR NOW!).
  DO  i=1,imax
     cr(i,1) = refl(i)
  END DO

  !OBTAIN CLOUD REFLECTION AND TRANSMISSION COEFFICIENTS FOR
  !REMAINING CLOUDS (IF ANY) IN THE VISIBLE BAND
  !THE MAXIMUM NO OF CLOUDS IN THE ROW (KCLDS) IS USED. THIS CREATES
  !EXTRA WORK (MAY BE REMOVED IN A SUBSEQUENT UPDATE).
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=2,kclds+1
        cr(i,kk) = crr(i,1,kk)*xamt(i,kk)
        ct(i,kk) = one - (one-ctt(i,1,kk))*xamt(i,kk)
     END DO
  END DO

  !OBTAIN THE PRESSURE AT THE TOP,BOTTOM AND THE THICKNESS OF
  !"THICK" CLOUDS (THOSE AT LEAST 2 LAYERS THICK). THIS IS USED
  !LATER IS OBTAINING FLUXES INSIDE THE THICK CLOUDS, FOR ALL
  !FREQUENCY BANDS.
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=1,kclds
        IF ((kbtmsw(i,kk+1)-1) > ktopsw(i,kk+1)) THEN
           pptop(i,kk)=pp(i,ktopsw(i,kk+1))
           dpcld(i,kk)=one/(pptop(i,kk)-pp(i,kbtmsw(i,kk+1)))
        END IF
     END DO
  END DO
  DO  i=1,imax*l
     ttdb1(i,2) = EXP(hm1ez*MIN(fifty,abcff(1)*ud(i,2)))
     ttub1(i,1) = EXP(hm1ez*MIN(fifty,abcff(1)*ur(i,1)))
     ttd  (i,2) = ttdb1(i,2)*tdo3(i,2)
     ttu  (i,1) = ttub1(i,1)*tuo3(i,1)
  END DO
  DO  i=1,imax
     ttd(i,1)   = one
     ttu(i,lp1) = ttd(i,lp1)
  END DO

  !FOR EXECUTION OF THE CLOUD LOOP, IT IS NECESSARY TO SEPARATE OUT
  !TRANSMISSION FCTNS AT THE TOP AND BOTTOM OF THE CLOUDS, FOR
  !EACH BAND N. THE REQUIRED QUANTITIES ARE:
  !   TTD(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
  !   TTU(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
  !   TTD(I,KBTMSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
  !   AND INVERSES OF THE FIRST TWO. THE ABOVE QUANTITIES ARE
  !   STORED IN TDCL1,TUCL1,TDCL2, AND DFNTRN,UFNTRN, RESPECTIVELY,
  !   AS THEY HAVE MULTIPLE USE IN THE PGM.
  !FOR FIRST CLOUD LAYER (GROUND) TDCL1,TUCL1 ARE KNOWN:
  DO  i=1,imax
     tdcl1 (i,1) = ttd(i,lp1)
     tucl1 (i,1) = ttu(i,lp1)
     tdcl2 (i,1) = tdcl1(i,1)
     dfntrn(i,1) = one/tdcl1(i,1)
     ufntrn(i,1) = dfntrn(i,1)
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=2,kclds+1
        tdcl1(i,kk) = ttd(i,ktopsw(i,kk))
        tucl1(i,kk) = ttu(i,ktopsw(i,kk))
        tdcl2(i,kk) = ttd(i,kbtmsw(i,kk))
     END DO
  END DO
  !COMPUTE INVERSES
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=2,kclds
        dfntrn(i,kk) = one/tdcl1(i,kk)
        ufntrn(i,kk) = one/tucl1(i,kk)
     END DO
  END DO

  !COMPUTE THE TRANSMISSIVITY FROM THE TOP OF CLOUD (K+1) TO THE
  !TOP OF CLOUD (K). THE CLOUD TRANSMISSION (CT) IS INCLUDED. THIS
  !QUANTITY IS CALLED TCLU (INDEX K). ALSO, OBTAIN THE TRANSMISSIVITY
  !FROM THE BOTTOM OF CLOUD (K+1) TO THE TOP OF CLOUD (K)(A PATH
  !ENTIRELY OUTSIDE CLOUDS). THIS QUANTITY IS CALLED TCLD (INDEX K).
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=1,kclds
        tclu(i,kk) = tdcl1(i,kk)*dfntrn(i,kk+1)*ct(i,kk+1)
        tcld(i,kk) = tdcl1(i,kk)/tdcl2(i,kk+1)
     END DO
  END DO

  !THE FOLLOWING IS THE RECURSION RELATION FOR ALFA: THE REFLECTION
  !COEFFICIENT FOR A SYSTEM INCLUDING THE CLOUD IN QUESTION AND THE
  !FLUX COMING OUT OF THE CLOUD SYSTEM INCLUDING ALL CLOUDS BELOW
  !THE CLOUD IN QUESTION.
  !ALFAU IS ALFA WITHOUT THE REFLECTION OF THE CLOUD IN QUESTION
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     alfa (i,1)=cr(i,1)
     alfau(i,1)=zero
  END DO
  !AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER!
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=2,kclds+1
        alfau(i,kk)= tclu(i,kk-1)*tclu(i,kk-1)*alfa(i,kk-1)/  &
             (one - tcld(i,kk-1)*tcld(i,kk-1)*alfa(i,kk-1)*cr(i,kk))
        alfa (i,kk)= alfau(i,kk)+cr(i,kk)
     END DO
  END DO

  !CALCULATE UFN AT CLOUD TOPS AND DFN AT CLOUD BOTTOMS
  !NOTE THAT UFNCLU(I,KCLDS+1) GIVES THE UPWARD FLUX AT THE TOP
  !OF THE HIGHEST REAL CLOUD (IF NCLDS(I)=KCLDS). IT GIVES THE FLUX
  !AT THE TOP OF THE ATMOSPHERE IF NCLDS(I) < KCLDS. IN THE FIRST
  !CASE, TDCL1 EQUALS THE TRANSMISSION FCTN TO THE TOP OF THE
  !HIGHEST CLOUD, AS WE WANT. IN THE SECOND CASE, TDCL1=1, SO UFNCLU
  !EQUALS ALFA. THIS IS ALSO CORRECT.
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     ufnclu(i,kclds+1) = alfa(i,kclds+1)*tdcl1(i,kclds+1)
     dfnclu(i,kclds+1) = tdcl1(i,kclds+1)
  END DO
  !THIS CALCULATION IS THE REVERSE OF THE RECURSION RELATION USED
  !ABOVE
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=kclds,1,-1
        ufnclu(i,kk) = ufnclu(i,kk+1)*alfau(i,kk+1)/(alfa(i,kk+1)* tclu(i,kk))
        dfnclu(i,kk) = ufnclu(i,kk)/alfa(i,kk)
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=1,kclds+1
        ufntrn(i,kk) = ufnclu(i,kk)*ufntrn(i,kk)
        dfntrn(i,kk) = dfnclu(i,kk)*dfntrn(i,kk)
     END DO
  END DO
  !CASE OF KK=1( FROM THE GROUND TO THE BOTTOM OF THE LOWEST CLOUD)
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     j2=kbtmsw(i,2)
     DO  k=j2,lp1
        ufn(i,k) = ufntrn(i,1)*ttu(i,k)
        dfn(i,k) = dfntrn(i,1)*ttd(i,k)
     END DO
  END DO
  !REMAINING LEVELS (IF ANY!)
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=2,kclds+1
        j1=ktopsw(i,kk)
        j2=kbtmsw(i,kk+1)
        IF (j1 == 1) CYCLE
        DO  k=j2,j1
           ufn(i,k) = ufntrn(i,kk)*ttu(i,k)
           dfn(i,k) = dfntrn(i,kk)*ttd(i,k)
        END DO
        !FOR THE THICK CLOUDS, THE FLUX DIVERGENCE THROUGH THE CLOUD
        !LAYER IS ASSUMED TO BE CONSTANT. THE FLUX DERIVATIVE IS GIVEN BY
        !TEMPF (FOR THE UPWARD FLUX) AND TEMPG (FOR THE DOWNWARD FLUX).
        j3=kbtmsw(i,kk)
        IF ((j3-j1) > 1) THEN
           tempf = (ufnclu(i,kk)-ufn(i,j3))*dpcld(i,kk-1)
           tempg = (dfnclu(i,kk)-dfn(i,j3))*dpcld(i,kk-1)
           DO  k=j1+1,j3-1
              ufn(i,k) = ufnclu(i,kk)+tempf*(pp(i,k)-pptop(i,kk-1))
              dfn(i,k) = dfnclu(i,kk)+tempg*(pp(i,k)-pptop(i,kk-1))
           END DO
        END IF
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  k=1,lp1
        dfswc(i,k) = dfn(i,k)*dfntop(i,1)
        ufswc(i,k) = ufn(i,k)*dfntop(i,1)
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     temp1(i) = one - ccmax(i)
     gdfvb(i) = temp1(i)*gdfvb(i)
     gdfnb(i) = temp1(i)*gdfnb(i)
     gdfvd(i) = temp1(i)*gdfvd(i) + ccmax(i)*dfswc(i,lp1)
  END DO

  !NOW OBTAIN FLUXES FOR THE NEAR IR BANDS. THE METHODS ARE THE SAME
  !AS FOR THE VISIBLE BAND, EXCEPT THAT THE REFLECTION AND
  !TRANSMISSION COEFFICIENTS ARE DIFFERENT, AS
  !RAYLEIGH SCATTERING NEED NOT BE CONSIDERED.
  DO  n=2,nb
     !YH93
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  k=1,kclds+1
           cr(i,k) = crr(i,n,k)*xamt(i,k)
           ct(i,k) = one - (one-ctt(i,n,k))*xamt(i,k)
        END DO
     END DO
     !YH93
     IF (n == 2) THEN
        !THE WATER VAPOR TRANSMISSION FUNCTION FOR BAND 2 IS EQUAL TO
        !THAT OF BAND 1 (SAVED AS TTDB1,TTUB1)
        DO  i=1,imax
           kclds=nclds(i)
           IF(kclds == 0) CYCLE
           DO  kk=2,lp1
              ttd(i,kk) = ttdb1(i,kk)*tdco2(i,kk)
           END DO
           DO  kk=1,l
              ttu(i,kk) = ttub1(i,kk)*tuco2(i,kk)
           END DO
        END DO
     ELSE
        DO  i=1,imax
           kclds=nclds(i)
           IF(kclds == 0) CYCLE
           DO  kk=2,lp1
              ttd(i,kk) = EXP(hm1ez*MIN(fifty,abcff(n)*ud(i,kk))) * tdco2(i,kk)
           END DO
           DO  kk=1,l
              ttu(i,kk) = EXP(hm1ez*MIN(fifty,abcff(n)*ur(i,kk))) * tuco2(i,kk)
           END DO
        END DO
     END IF
     !AT THIS POINT,INCLUDE TTD(1),TTU(LP1), NOTING THAT TTD(1)=1 FOR
     !ALL BANDS, AND THAT TTU(LP1)=TTD(LP1) FOR ALL BANDS.
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        ttu(i,lp1) = ttd(i,lp1)
        ttd(i,1)   = one
     END DO
     !FOR EXECUTION OF THE CLOUD LOOP, IT IS NECESSARY TO SEPARATE OUT
     !TRANSMISSION FCTNS AT THE TOP AND BOTTOM OF THE CLOUDS, FOR
     !EACH BAND N. THE REQUIRED QUANTITIES ARE:
     !  TTD(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
     !  TTD(I,KBTMSW(I,K),N)  K RUNS FROM 2 TO NCLDS(I)+1:
     !  TTU(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
     !  AND INVERSES OF THE ABOVE. THE ABOVE QUANTITIES ARE STORED
     !  IN TDCL1,TDCL2,TUCL1,AND DFNTRN,UFNTRN,RESPECTIVELY, AS
     !   THEY HAVE MULTIPLE USE IN THE PGM.
     !FOR FIRST CLOUD LAYER (GROUND) TDCL1,TUCL1 ARE KNOWN:
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        tdcl1 (i,1) = ttd(i,lp1)
        tucl1 (i,1) = ttu(i,lp1)
        tdcl2 (i,1) = tdcl1(i,1)
        dfntrn(i,1) = one/tdcl1(i,1)
        ufntrn(i,1) = dfntrn(i,1)
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=2,kclds+1
           tdcl1(i,kk) = ttd(i,ktopsw(i,kk))
           tucl1(i,kk) = ttu(i,ktopsw(i,kk))
           tdcl2(i,kk) = ttd(i,kbtmsw(i,kk))
        END DO
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=2,kclds+1
           dfntrn(i,kk) = one/tdcl1(i,kk)
           ufntrn(i,kk) = one/tucl1(i,kk)
        END DO
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=1,kclds
           tclu(i,kk) = tdcl1(i,kk)*dfntrn(i,kk+1)*ct(i,kk+1)
           tcld(i,kk) = tdcl1(i,kk)/tdcl2(i,kk+1)
        END DO
     END DO
     !THE FOLLOWING IS THE RECURSION RELATION FOR ALFA: THE REFLECTION
     !COEFFICIENT FOR A SYSTEM INCLUDING THE CLOUD IN QUESTION AND THE
     !FLUX COMING OUT OF THE CLOUD SYSTEM INCLUDING ALL CLOUDS BELOW
     !THE CLOUD IN QUESTION.
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        alfa (i,1) = cr(i,1)
        alfau(i,1) = zero
     END DO
     !AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER!
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=2,kclds+1
           alfau(i,kk) = tclu(i,kk-1)*tclu(i,kk-1)*alfa(i,kk-1)/(one -  &
                tcld(i,kk-1)*tcld(i,kk-1)*alfa(i,kk-1)*cr(i,kk))
           alfa (i,kk) = alfau(i,kk)+cr(i,kk)
        END DO
     END DO
     !CALCULATE UFN AT CLOUD TOPS AND DFN AT CLOUD BOTTOMS
     !NOTE THAT UFNCLU(I,KCLDS+1) GIVES THE UPWARD FLUX AT THE TOP
     !OF THE HIGHEST REAL CLOUD (IF NCLDS(I)=KCLDS). IT GIVES THE FLUX
     !AT THE TOP OF THE ATMOSPHERE IF NCLDS(I) < KCLDS. IT THE FIRST
     !CASE, TDCL1 EQUALS THE TRANSMISSION FCTN TO THE TOP OF THE
     !HIGHEST CLOUD, AS WE WANT. IN THE SECOND CASE, TDCL1=1, SO UFNCLU
     !EQUALS ALFA. THIS IS ALSO CORRECT.
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        ufnclu(i,kclds+1) = alfa(i,kclds+1)*tdcl1(i,kclds+1)
        dfnclu(i,kclds+1) = tdcl1(i,kclds+1)
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=kclds,1,-1
           ufnclu(i,kk)=ufnclu(i,kk+1)*alfau(i,kk+1)/(alfa(i,kk+1)* tclu(i,kk))
           dfnclu(i,kk) = ufnclu(i,kk)/alfa(i,kk)
        END DO
     END DO
     !NOW OBTAIN DFN AND UFN FOR LEVELS BETWEEN THE CLOUDS
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  kk=1,kclds+1
           ufntrn(i,kk) = ufnclu(i,kk)*ufntrn(i,kk)
           dfntrn(i,kk) = dfnclu(i,kk)*dfntrn(i,kk)
        END DO
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        j2=kbtmsw(i,2)
        DO  k=j2,lp1
           ufn(i,k) = ufntrn(i,1)*ttu(i,k)
           dfn(i,k) = dfntrn(i,1)*ttd(i,k)
        END DO
     END DO
     DO   i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO   kk=2,kclds+1
           j1 = ktopsw(i,kk)
           j2 = kbtmsw(i,kk+1)
           IF (j1 == 1) CYCLE
           DO  k=j2,j1
              ufn(i,k) = ufntrn(i,kk)*ttu(i,k)
              dfn(i,k) = dfntrn(i,kk)*ttd(i,k)
           END DO
           j3 = kbtmsw(i,kk)
           IF ((j3-j1) > 1) THEN
              tempf = (ufnclu(i,kk)-ufn(i,j3))*dpcld(i,kk-1)
              tempg = (dfnclu(i,kk)-dfn(i,j3))*dpcld(i,kk-1)
              DO  k=j1+1,j3-1
                 ufn(i,k) = ufnclu(i,kk)+tempf*(pp(i,k)-pptop(i,kk-1))
                 dfn(i,k) = dfnclu(i,kk)+tempg*(pp(i,k)-pptop(i,kk-1))
              END DO
           END IF
        END DO
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        DO  k=1,lp1
           dfswc(i,k) = dfswc(i,k) + dfn(i,k)*dfntop(i,n)
           ufswc(i,k) = ufswc(i,k) + ufn(i,k)*dfntop(i,n)
        END DO
     END DO
     DO  i=1,imax
        kclds=nclds(i)
        IF(kclds == 0) CYCLE
        gdfnd(i) = gdfnd(i) + ccmax(i)*dfn(i,lp1)*dfntop(i,n)
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  k=1,lp1
        dfswc(i,k) = temp1(i)*dfswl(i,k) + ccmax(i)*dfswc(i,k)
        ufswc(i,k) = temp1(i)*ufswl(i,k) + ccmax(i)*ufswc(i,k)
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=1,lp1
        fswc(i,kk) = ufswc(i,kk)-dfswc(i,kk)
     END DO
  END DO
  DO  i=1,imax
     kclds=nclds(i)
     IF(kclds == 0) CYCLE
     DO  kk=1, l
        hswc(i,kk) = radcon*(fswc(i,kk+1)-fswc(i,kk))/dp(i,kk)
     END DO
  END DO
END SUBROUTINE swr93
