      SUBROUTINE SWR93(FSWC,HSWC,UFSWC,DFSWC,FSWL,HSWL,UFSWL,
     &                 DFSWL,
     1                 PRESS,COSZRO,TAUDAR,RH2O,RRCO2,SSOLAR,QO3,
     2                 NCLDS,KTOPSW,KBTMSW,CAMT,CRR,CTT,
     A                 ALVB,ALNB,ALVD,ALND,GDFVB,GDFNB,GDFVD,GDFND)
C===>    *********************************************************
C --- SWR91SIB --- MODIFIED FROM SWR89-BAND12....YUTAI HOU
C
C     -SW- RADIATION CODE............................
C        INPUTS:PRESS,COSZRO,TAUDAR,RH2O,RRCO2,SSOLAR,QO3,NCLDS,
C               KTOPSW,KBTMSW,CIRAB,CIRRF,CUVRF,CAMT,
C               ALVB,ALVD,ALNB,ALND;
C        OUTPUT:FSWC,HSWC,UFSWC,DFSWC,FSWL,HSWL,UFSWL,DFSWL,
C               GDFVB,GDFVD,GDFNB,GDFND.
C --- SWR91SIB --- MODIFIED BY K. CAMPANA..06 MAR 92
C         INCLUDE HPCON,PARMC CHANGED TO HCON,RDPARM
C         6 Q..... VARIABLE NAMES RESTORED TO ORIGINAL 7,8 CHAR
C         CHANGE O3DIFF,DIFFCC TO O3DIFCTR,DIFFCTR
C --- SWR91SIB --- MODIFIED BY Y. HOU         FEB 93
C         INPUTS 12 BANDS CLD REFLECTTANCE AND TRANSMITTANCE
C         CRR,CTT TO REPLACE CIRAB,CIRRF,CUVRF
C===>    *********************************************************
C
C
        INCLUDE "parmeta"
        INCLUDE "HCON.comm" 
        INCLUDE "rdparm" 
        INCLUDE "mpp.h"
#include "sp.h"
C     PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE:
C          L      =  NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL
C          NB IS A SHORTWAVE PARAMETER; OTHER QUANTITIES ARE DERIVED
C                    FROM THE ABOVE PARAMETERS.
C --- VARIABLES AS IN ARGUMENT LIST
                       D I M E N S I O N
     1  FSWC  (IDIM1:IDIM2,LP1),  HSWC  (IDIM1:IDIM2,LP1),
     2  CRR   (IDIM1:IDIM2,NB,LP1)
     3, FSWL  (IDIM1:IDIM2,LP1),  HSWL  (IDIM1:IDIM2,LP1),
     4  CTT   (IDIM1:IDIM2,NB,LP1)
     5, UFSWC (IDIM1:IDIM2,LP1),  DFSWC (IDIM1:IDIM2,LP1)
     6, UFSWL (IDIM1:IDIM2,LP1),  DFSWL (IDIM1:IDIM2,LP1)
     7, PRESS (IDIM1:IDIM2,LP1),  RH2O  (IDIM1:IDIM2,L),
     8  QO3   (IDIM1:IDIM2,L)
     9, CAMT  (IDIM1:IDIM2,LP1),  KTOPSW(IDIM1:IDIM2,LP1),
     o  KBTMSW(IDIM1:IDIM2,LP1)
     1, COSZRO(IDIM1:IDIM2),      TAUDAR(IDIM1:IDIM2),
     2  NCLDS (IDIM1:IDIM2)
     3, ALVB  (IDIM1:IDIM2),  ALNB  (IDIM1:IDIM2),
     4  ALVD  (IDIM1:IDIM2),  ALND  (IDIM1:IDIM2)
     5, GDFVB (IDIM1:IDIM2),  GDFNB (IDIM1:IDIM2),  GDFVD (IDIM1:IDIM2),
     6  GDFND (IDIM1:IDIM2)
C --- LOCAL VARIABLES
                       D I M E N S I O N
     1  PP    (IDIM1:IDIM2,LP1),  DP    (IDIM1:IDIM2,LP1),
     2  PR2   (IDIM1:IDIM2,LP1)
     3, DU    (IDIM1:IDIM2,LP1),  DUCO2 (IDIM1:IDIM2,LP1),
     4  DUO3  (IDIM1:IDIM2,LP1)
     5, FF    (IDIM1:IDIM2,LP1),  FFCO2 (IDIM1:IDIM2,LP1),
     6  FFO3  (IDIM1:IDIM2,LP1)
     7, RRAY  (IDIM1:IDIM2),      DFNTOP(IDIM1:IDIM2,NB),
     8  SECZ  (IDIM1:IDIM2)
     9, REFL  (IDIM1:IDIM2),      TMP1 (IDIM1:IDIM2),
     o  REFL2 (IDIM1:IDIM2)
     1, CCMAX (IDIM1:IDIM2),      XAMT  (IDIM1:IDIM2,LP1)
                       D I M E N S I O N
     1  UD    (IDIM1:IDIM2,LP1),   UR    (IDIM1:IDIM2,LP1)
     2, UCO2  (IDIM1:IDIM2,LLP2), UDCO2 (IDIM1:IDIM2,LP1),
     3  URCO2 (IDIM1:IDIM2,LP1)
     4, UO3   (IDIM1:IDIM2,LLP2), UDO3  (IDIM1:IDIM2,LP1),
     5  URO3  (IDIM1:IDIM2,LP1)
     6, TCO2  (IDIM1:IDIM2,LLP2), TDCO2 (IDIM1:IDIM2,LP1),
     7  TUCO2 (IDIM1:IDIM2,LP1)
     8, TO3   (IDIM1:IDIM2,LLP2), TDO3  (IDIM1:IDIM2,LP1),
     9  TUO3  (IDIM1:IDIM2,LP1)
                       D I M E N S I O N
     1  DFN   (IDIM1:IDIM2,LP1),  UFN   (IDIM1:IDIM2,LP1),
     2  CR    (IDIM1:IDIM2,LP1)
     3, TTD   (IDIM1:IDIM2,LP1),  TTU   (IDIM1:IDIM2,LP1),
     4  CT    (IDIM1:IDIM2,LP1)
     5, PPTOP (IDIM1:IDIM2,LP1),  DPCLD (IDIM1:IDIM2,LP1)
C --- EQUIVALENCED LOCAL VARIABLES
                       D I M E N S I O N
     1  TTUB1 (IDIM1:IDIM2,LP1),  TUCL1 (IDIM1:IDIM2,LP1)
     2, TTDB1 (IDIM1:IDIM2,LP1),  TDCL1 (IDIM1:IDIM2,LP1),
     3  TDCL2 (IDIM1:IDIM2,LP1)
     4, UFNTRN(IDIM1:IDIM2,LP1),  UFNCLU(IDIM1:IDIM2,LP1),
     5  TCLU  (IDIM1:IDIM2,LP1)
     6, DFNTRN(IDIM1:IDIM2,LP1),  DFNCLU(IDIM1:IDIM2,LP1),
     7  TCLD  (IDIM1:IDIM2,LP1)
     8, ALFA  (IDIM1:IDIM2,LP1),  ALFAU (IDIM1:IDIM2,LP1)
                     E Q U I V A L E N C E
     1  (UDO3,UO3(IDIM1,1),DFNCLU), (URO3,UO3(IDIM1,LP2), UFNCLU)
     2, (UDCO2,UCO2(IDIM1,1),TCLD), (URCO2,UCO2(IDIM1,LP2), TCLU)
     3, (TDO3 ,TO3(IDIM1,1),DFNTRN),(TUO3,TO3(IDIM1,LP2), UFNTRN)
     4, (TDCO2,TCO2(IDIM1,1)      ),(TUCO2,TCO2(IDIM1,LP2)        )
     5, (FF   , ALFA ),   (FFCO2 , ALFAU ),   (FFO3  , TTDB1 )
     6, (DU   , TTUB1),   (DUCO2 , TUCL1 ),   (DUO3  , TDCL1 )
     7, (PR2  , TDCL2)
C
C---COMMON FOR LOCAL DATA VARIABLES---
      COMMON /SWRSAV/ ABCFF(NB),PWTS(NB),CFCO2,CFO3,REFLO3,RRAYAV
C                           D A T A
C    1  ABCFF / 2*4.0E-5, 0.002, 0.035, 0.377, 1.95, 9.40, 44.6,
C    1          190.0,    989.0, 2706.0, 39011.0 /
C    2, PWTS  / 0.5000, 0.121416, 0.0698, 0.1558, 0.0631, 0.0362,
C    2          0.0243, 0.0158, 0.0087, 0.001467, 0.002342, 0.001075 /
C    3, CFCO2, CFO3, REFLO3, RRAYAV / 508.96, 466.64, 1.9, 0.144 /
C    1  ABCFF / 2*4.0E-5, .002, .035, .377, 1.95, 9.40, 44.6, 190. /
C    2, PWTS  /.5000,.1470,.698,.1443,.0584,.0335,.0225,.0158,.0087/
C    3, CFCO2, CFO3, REFLO3, RRAYAV / 508.96, 466.64, 1.9, 0.144 /
C
C     CALCULATE SECANT OF ZENITH ANGLE (SECZ),FLUX PRESSURES(PP),
C     LAYER WIDTH (DP) AND PRESSURE SCALING FACTOR (PR2).
      DO 100 I=MYIS,MYIE
        SECZ(I) = H35E1/SQRT(H1224E3*COSZRO(I)*COSZRO(I)+ONE)
        PP(I,1)   = ZERO
        PP(I,LP1) = PRESS(I,LP1)
        TMP1(I)  = ONE/PRESS(I,LP1)
  100  CONTINUE
      DO 110 K=1,LM1
      DO 110 I=MYIS,MYIE
        PP(I,K+1) = HAF*(PRESS(I,K+1)+PRESS(I,K))
  110  CONTINUE
      DO 120 K=1,L
      DO 120 I=MYIS,MYIE
        DP (I,K) = PP(I,K+1)-PP(I,K)
        PR2(I,K) = HAF*(PP(I,K)+PP(I,K+1))
  120  CONTINUE
      DO 130 K=1,L
      DO 130 I=MYIS,MYIE
        PR2(I,K) = PR2(I,K)*TMP1(I)
  130  CONTINUE
C     CALCULATE ENTERING FLUX AT THE TOP FOR EACH BAND(IN CGS UNITS)
      DO 140 N=1,NB
      DO 140 IP=MYIS,MYIE
        DFNTOP(IP,N) = SSOLAR*H69766E5*COSZRO(IP)*TAUDAR(IP)*PWTS(N)
  140  CONTINUE
C     EXECUTE THE LACIS-HANSEN REFLECTIVITY PARAMETERIZATION
C     FOR THE VISIBLE BAND
      DO 150 I=MYIS,MYIE
        RRAY(I) = HP219/(ONE+HP816*COSZRO(I))
        REFL(I) = RRAY(I) + (ONE-RRAY(I))*(ONE-RRAYAV)*ALVB(I)/
     1            (ONE-ALVD(I)*RRAYAV)
  150  CONTINUE
      DO 155 I=MYIS,MYIE
        RRAY(I) = 0.104/(ONE+4.8*COSZRO(I))
        REFL2(I)= RRAY(I) + (ONE-RRAY(I))*(ONE-0.093)*ALVB(I)/
     1            (ONE-ALVD(I)*0.093)
  155  CONTINUE
C     CALCULATE PRESSURE-WEIGHTED OPTICAL PATHS FOR EACH LAYER
C     IN UNITS OF CM-ATM. PRESSURE WEIGHTING IS USING PR2.
C     DU= VALUE FOR H2O;DUCO2 FOR CO2;DUO3 FOR O3.
      DO 160 K=1,L
      DO 160 I=MYIS,MYIE
        DU   (I,K) = GINV*RH2O(I,K)*DP(I,K)*PR2(I,K)
        DUCO2(I,K) = (RRCO2*GINV*CFCO2)*DP(I,K)*PR2(I,K)
        DUO3 (I,K) = (GINV*CFO3)*QO3(I,K)*DP(I,K)
  160  CONTINUE
C
C                 CALCULATE CLEAR SKY SW FLUX
C
C     OBTAIN THE OPTICAL PATH FROM THE TOP OF THE ATMOSPHERE TO THE
C     FLUX PRESSURE. ANGULAR FACTORS ARE NOW INCLUDED. UD=DOWNWARD
C     PATH FOR H2O,WIGTH UR THE UPWARD PATH FOR H2O. CORRESPONDING
C     QUANTITIES FOR CO2,O3 ARE UDCO2/URCO2 AND UDO3/URO3.
      DO 200 IP=MYIS,MYIE
        UD   (IP,1) = ZERO
        UDCO2(IP,1) = ZERO
        UDO3 (IP,1) = ZERO
  200  CONTINUE
      DO 210 K=2,LP1
      DO 210 I=MYIS,MYIE
        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)
  210  CONTINUE
      DO 220 IP=MYIS,MYIE
        UR   (IP,LP1) = UD   (IP,LP1)
        URCO2(IP,LP1) = UDCO2(IP,LP1)
        URO3 (IP,LP1) = UDO3 (IP,LP1)
  220  CONTINUE
      DO 230 K=L,1,-1
      DO 230 IP=MYIS,MYIE
        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
  230  CONTINUE
C     CALCULATE CO2 ABSORPTIONS . THEY WILL BE USED IN NEAR INFRARED
C     BANDS.SINCE THE ABSORPTION AMOUNT IS GIVEN (IN THE FORMULA USED
C     BELOW, DERIVED FROM SASAMORI) IN TERMS OF THE TOTAL SOLAR FLUX,
C     AND THE ABSORPTION IS ONLY INCLUDED IN THE NEAR IR (50 PERCENT
C     OF THE SOLAR SPECTRUM), THE ABSORPTIONS ARE MULTIPLIED BY 2.
C       SINCE CODE ACTUALLY REQUIRES TRANSMISSIONS, THESE ARE THE
C     VALUES ACTUALLY STORED IN TCO2.
      DO 240 K=1,LL
      DO 240 I=MYIS,MYIE
       TCO2(I,K+1)=ONE-TWO*(H235M3*EXP(HP26*LOG(UCO2(I,K+1)+H129M2))
     1                       -H75826M4)
  240  CONTINUE
C     NOW CALCULATE OZONE ABSORPTIONS. THESE WILL BE USED IN
C     THE VISIBLE BAND.JUST AS IN THE CO2 CASE, SINCE THIS BAND IS
C     50 PERCENT OF THE SOLAR SPECTRUM,THE ABSORPTIONS ARE MULTIPLIED
C     BY 2. THE TRANSMISSIONS ARE STORED IN TO3.
      HTEMP = H1036E2*H1036E2*H1036E2
      DO 250 K=1,LL
      DO 250 I=MYIS,MYIE
        TO3(I,K+1)=ONE-TWO*UO3(I,K+1)*
     1            (H1P082*EXP(HMP805*LOG(ONE+H1386E2*UO3(I,K+1)))+
     2            H658M2/(ONE+HTEMP*UO3(I,K+1)*UO3(I,K+1)*UO3(I,K+1))+
     3            H2118M2/(ONE+UO3(I,K+1)*(H42M2+H323M4*UO3(I,K+1))))
  250  CONTINUE
C   START FREQUENCY LOOP (ON N) HERE
C
C--- BAND 1 (VISIBLE) INCLUDES O3 AND H2O ABSORPTION
      DO 260 K=1,L
      DO 260 I=MYIS,MYIE
        TTD(I,K+1) = EXP(HM1EZ*MIN(FIFTY,ABCFF(1)*UD(I,K+1)))
        TTU(I,K) = EXP(HM1EZ*MIN(FIFTY,ABCFF(1)*UR(I,K)))
        DFN(I,K+1) = TTD(I,K+1)*TDO3(I,K+1)
        UFN(I,K) = TTU(I,K)*TUO3(I,K)
  260  CONTINUE
      DO 270 I=MYIS,MYIE
        DFN(I,1)   = ONE
        UFN(I,LP1) = DFN(I,LP1)
  270  CONTINUE
C     SCALE VISIBLE BAND FLUXES BY SOLAR FLUX AT THE TOP OF THE
C     ATMOSPHERE (DFNTOP(I,1))
C     DFSW/UFSW WILL BE THE FLUXES, SUMMED OVER ALL BANDS
      DO 280  K=1,LP1
      DO 280  I=MYIS,MYIE
        DFSWL(I,K) =         DFN(I,K)*DFNTOP(I,1)
        UFSWL(I,K) = REFL(I)*UFN(I,K)*DFNTOP(I,1)
  280  CONTINUE
      DO 285 I=MYIS,MYIE
        GDFVB(I) = DFSWL(I,LP1)*EXP(-0.15746*SECZ(I))
	if (ALVD(I) .eq. 1) then
	write(6,*) 'heading for trouble'
	STOP 666
	endif
        GDFVD(I) = ((ONE-REFL2(I))*DFSWL(I,LP1) -
     1              (ONE-ALVB(I)) *GDFVB(I)) / (ONE-ALVD(I))
!	write(6,*) 'GDFVD val: ', GDFVD(I)
Cmp
	if (GDFVD(I) .gt. -1.e+20 .and. GDFVD(I) .lt. 1.e+20) then
	else
	write(6,*) 'crazy GDFVD val!'
	write(6,*) 'pieces: one,refl2,dfswl,alvb,gdfvb,alvd ',
     +	 one,refl2(i),dfswl(i,lp1),alvb(i),gdfvb(i),alvd(i)
	endif
        GDFNB(I) = ZERO
        GDFND(I) = ZERO
  285  CONTINUE
C---NOW OBTAIN FLUXES FOR THE NEAR IR BANDS. THE METHODS ARE THE SAME
C   AS FOR THE VISIBLE BAND, EXCEPT THAT THE REFLECTION AND
C   TRANSMISSION COEFFICIENTS (OBTAINED BELOW) ARE DIFFERENT, AS
C   RAYLEIGH SCATTERING NEED NOT BE CONSIDERED.
      DO 350 N=2,NB
        IF (N.EQ.2) THEN
C   THE WATER VAPOR TRANSMISSION FUNCTION FOR BAND 2 IS EQUAL TO
C   THAT OF BAND 1 (SAVED AS TTD,TTU)
C--- BAND 2-9 (NEAR-IR) INCLUDES O3, CO2 AND H2O ABSORPTION
          DO 290 K=1,L
          DO 290 I=MYIS,MYIE
            DFN(I,K+1) = TTD(I,K+1)*TDCO2(I,K+1)
            UFN(I,K) = TTU(I,K)*TUCO2(I,K)
  290     CONTINUE
        ELSE
C   CALCULATE WATER VAPOR TRANSMISSION FUNCTIONS FOR NEAR INFRARED
C   BANDS. INCLUDE CO2 TRANSMISSION (TDCO2/TUCO2), WHICH
C   IS THE SAME FOR ALL INFRARED BANDS.
          DO 300 K=1,L
          DO 300 I=MYIS,MYIE
            DFN(I,K+1)=EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UD(I,K+1)))
     1                 *TDCO2(I,K+1)
            UFN(I,K)=EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UR(I,K)))
     1               *TUCO2(I,K)
  300     CONTINUE
        ENDIF
C---AT THIS POINT,INCLUDE DFN(1),UFN(LP1), NOTING THAT DFN(1)=1 FOR
C   ALL BANDS, AND THAT UFN(LP1)=DFN(LP1) FOR ALL BANDS.
        DO 310 I=MYIS,MYIE
          DFN(I,1)   = ONE
          UFN(I,LP1) = DFN(I,LP1)
  310   CONTINUE
C     SCALE THE PREVIOUSLY COMPUTED FLUXES BY THE FLUX AT THE TOP
C     AND SUM OVER BANDS
        DO 320 K=1,LP1
        DO 320 I=MYIS,MYIE
          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)
  320   CONTINUE
        DO 330 I=MYIS,MYIE
          GDFNB(I) = GDFNB(I) + DFN(I,LP1)*DFNTOP(I,N)
  330   CONTINUE
  350  CONTINUE
      DO 360 K=1,LP1
      DO 360 I=MYIS,MYIE
        FSWL(I,K) = UFSWL(I,K)-DFSWL(I,K)
  360  CONTINUE
      DO 370 K=1,L
      DO 370 I=MYIS,MYIE
        HSWL(I,K)=RADCON*(FSWL(I,K+1)-FSWL(I,K))/DP(I,K)
  370  CONTINUE
C
C---END OF FREQUENCY LOOP (OVER N)
C
C                 CALCULATE CLOUDY SKY SW FLUX
C
      KCLDS=NCLDS(MYIS)
      DO 400 I=MYIS1,MYIE
        KCLDS=MAX(NCLDS(I),KCLDS)
  400  CONTINUE
        DO 410 K=1,LP1
        DO 410 I=MYIS,MYIE
          DFSWC(I,K) = DFSWL(I,K)
          UFSWC(I,K) = UFSWL(I,K)
          FSWC (I,K) = FSWL (I,K)
  410   CONTINUE
        DO 420 K=1,L
        DO 420 I=MYIS,MYIE
          HSWC(I,K) = HSWL(I,K)
  420   CONTINUE
C*******************************************************************
      IF (KCLDS .EQ. 0)  RETURN
C*******************************************************************
      DO 430 K=1,LP1
      DO 430 I=MYIS,MYIE
        XAMT(I,K) = CAMT(I,K)
  430  CONTINUE
      DO 470 I=MYIS,MYIE
        NNCLDS   = NCLDS(I)
        CCMAX(I) = ZERO
        IF (NNCLDS .LE. 0) GO TO 470
        CCMAX(I) = ONE
        DO 450 K=1,NNCLDS
          CCMAX(I) = CCMAX(I) * (ONE - CAMT(I,K+1))
  450   CONTINUE
        CCMAX(I) = ONE - CCMAX(I)
        IF (CCMAX(I) .GT. ZERO) THEN
          DO 460 K=1,NNCLDS
            XAMT(I,K+1) = CAMT(I,K+1)/CCMAX(I)
  460     CONTINUE
        END IF
  470  CONTINUE
      DO 480 K=1,LP1
      DO 480 I=MYIS,MYIE
        FF   (I,K) = DIFFCTR
        FFCO2(I,K) = DIFFCTR
        FFO3 (I,K) = O3DIFCTR
  480  CONTINUE
      DO 490 IP=MYIS,MYIE
        JTOP = KTOPSW(IP,NCLDS(IP)+1)
      DO 490 K=1,JTOP
        FF   (IP,K) = SECZ(IP)
        FFCO2(IP,K) = SECZ(IP)
        FFO3 (IP,K) = SECZ(IP)
  490  CONTINUE
      DO 500 I=MYIS,MYIE
        RRAY(I) = HP219/(ONE+HP816*COSZRO(I))
        REFL(I) = RRAY(I) + (ONE-RRAY(I))*(ONE-RRAYAV)*ALVD(I)/
     1            (ONE-ALVD(I)*RRAYAV)
  500  CONTINUE
      DO 510 IP=MYIS,MYIE
        UD   (IP,1) = ZERO
        UDCO2(IP,1) = ZERO
        UDO3 (IP,1) = ZERO
  510  CONTINUE
      DO 520 K=2,LP1
      DO 520 I=MYIS,MYIE
        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)
  520  CONTINUE
      DO 530 IP=MYIS,MYIE
        UR   (IP,LP1) = UD   (IP,LP1)
        URCO2(IP,LP1) = UDCO2(IP,LP1)
        URO3 (IP,LP1) = UDO3 (IP,LP1)
  530  CONTINUE
      DO 540 K=L,1,-1
      DO 540 IP=MYIS,MYIE
        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
  540  CONTINUE
      DO 550 K=1,LL
      DO 550 I=MYIS,MYIE
        TCO2(I,K+1)=ONE-TWO*(H235M3*EXP(HP26*LOG(UCO2(I,K+1)+H129M2))
     1                        -H75826M4)
  550  CONTINUE
      DO 560 K=1,LL
      DO 560 I=MYIS,MYIE
        TO3(I,K+1)=ONE-TWO*UO3(I,K+1)*
     1           (H1P082*EXP(HMP805*LOG(ONE+H1386E2*UO3(I,K+1)))+
     2          H658M2/(ONE+HTEMP*UO3(I,K+1)*UO3(I,K+1)*UO3(I,K+1))+
     3          H2118M2/(ONE+UO3(I,K+1)*(H42M2+H323M4*UO3(I,K+1))))
  560  CONTINUE
C********************************************************************
C---THE FIRST CLOUD IS THE GROUND; ITS PROPERTIES ARE GIVEN
C   BY REFL (THE TRANSMISSION (0) IS IRRELEVANT FOR NOW!).
C********************************************************************
      DO 570 I=MYIS,MYIE
        CR(I,1) = REFL(I)
  570  CONTINUE
C***OBTAIN CLOUD REFLECTION AND TRANSMISSION COEFFICIENTS FOR
C   REMAINING CLOUDS (IF ANY) IN THE VISIBLE BAND
C---THE MAXIMUM NO OF CLOUDS IN THE ROW (KCLDS) IS USED. THIS CREATES
C   EXTRA WORK (MAY BE REMOVED IN A SUBSEQUENT UPDATE).
      DO 581 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 581
      DO 580 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)
  580  CONTINUE
  581  CONTINUE
C---OBTAIN THE PRESSURE AT THE TOP,BOTTOM AND THE THICKNESS OF
C   "THICK" CLOUDS (THOSE AT LEAST 2 LAYERS THICK). THIS IS USED
C   LATER IS OBTAINING FLUXES INSIDE THE THICK CLOUDS, FOR ALL
C   FREQUENCY BANDS.
      DO 591 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 591
      DO 590 KK=1,KCLDS
        IF ((KBTMSW(I,KK+1)-1).GT.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)))
        ENDIF
  590  CONTINUE
  591  CONTINUE
      DO 600 K=1,L
      DO 600 I=MYIS,MYIE
        TTDB1(I,K+1) = EXP(HM1EZ*MIN(FIFTY,ABCFF(1)*UD(I,K+1)))
        TTUB1(I,K) = EXP(HM1EZ*MIN(FIFTY,ABCFF(1)*UR(I,K)))
        TTD  (I,K+1) = TTDB1(I,K+1)*TDO3(I,K+1)
        TTU  (I,K) = TTUB1(I,K)*TUO3(I,K)
  600  CONTINUE
      DO 610 I=MYIS,MYIE
        TTD(I,1)   = ONE
        TTU(I,LP1) = TTD(I,LP1)
  610  CONTINUE
C***FOR EXECUTION OF THE CLOUD LOOP, IT IS NECESSARY TO SEPARATE OUT
C   TRANSMISSION FCTNS AT THE TOP AND BOTTOM OF THE CLOUDS, FOR
C   EACH BAND N. THE REQUIRED QUANTITIES ARE:
C      TTD(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
C      TTU(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
C      TTD(I,KBTMSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
C      AND INVERSES OF THE FIRST TWO. THE ABOVE QUANTITIES ARE
C      STORED IN TDCL1,TUCL1,TDCL2, AND DFNTRN,UFNTRN, RESPECTIVELY,
C      AS THEY HAVE MULTIPLE USE IN THE PGM.
C---FOR FIRST CLOUD LAYER (GROUND) TDCL1,TUCL1 ARE KNOWN:
      DO 620 I=MYIS,MYIE
        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)
  620  CONTINUE
      DO 631 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 631
      DO 630 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))
  630  CONTINUE
  631  CONTINUE
C---COMPUTE INVERSES
      DO 641 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 641
      DO 640 KK=2,KCLDS
        DFNTRN(I,KK) = ONE/TDCL1(I,KK)
        UFNTRN(I,KK) = ONE/TUCL1(I,KK)
  640  CONTINUE
  641  CONTINUE
C---COMPUTE THE TRANSMISSIVITY FROM THE TOP OF CLOUD (K+1) TO THE
C   TOP OF CLOUD (K). THE CLOUD TRANSMISSION (CT) IS INCLUDED. THIS
C   QUANTITY IS CALLED TCLU (INDEX K). ALSO, OBTAIN THE TRANSMISSIVITY
C   FROM THE BOTTOM OF CLOUD (K+1) TO THE TOP OF CLOUD (K)(A PATH
C   ENTIRELY OUTSIDE CLOUDS). THIS QUANTITY IS CALLED TCLD (INDEX K).
      DO 651 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 651
      DO 650 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)
  650  CONTINUE
  651  CONTINUE
C***THE FOLLOWING IS THE RECURSION RELATION FOR ALFA: THE REFLECTION
C   COEFFICIENT FOR A SYSTEM INCLUDING THE CLOUD IN QUESTION AND THE
C   FLUX COMING OUT OF THE CLOUD SYSTEM INCLUDING ALL CLOUDS BELOW
C   THE CLOUD IN QUESTION.
C---ALFAU IS ALFA WITHOUT THE REFLECTION OF THE CLOUD IN QUESTION
      DO 660 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 660
        ALFA (I,1)=CR(I,1)
        ALFAU(I,1)=ZERO
  660  CONTINUE
C---AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER!
      DO 671 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 671
      DO 670 KK=2,KCLDS+1
        ALFAU(I,KK)= TCLU(I,KK-1)*TCLU(I,KK-1)*ALFA(I,KK-1)/
     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)
  670  CONTINUE
  671  CONTINUE
C     CALCULATE UFN AT CLOUD TOPS AND DFN AT CLOUD BOTTOMS
C---NOTE THAT UFNCLU(I,KCLDS+1) GIVES THE UPWARD FLUX AT THE TOP
C   OF THE HIGHEST REAL CLOUD (IF NCLDS(I)=KCLDS). IT GIVES THE FLUX
C   AT THE TOP OF THE ATMOSPHERE IF NCLDS(I) < KCLDS. IN THE FIRST
C   CASE, TDCL1 EQUALS THE TRANSMISSION FCTN TO THE TOP OF THE
C   HIGHEST CLOUD, AS WE WANT. IN THE SECOND CASE, TDCL1=1, SO UFNCLU
C   EQUALS ALFA. THIS IS ALSO CORRECT.
      DO 680 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 680
        UFNCLU(I,KCLDS+1) = ALFA(I,KCLDS+1)*TDCL1(I,KCLDS+1)
        DFNCLU(I,KCLDS+1) = TDCL1(I,KCLDS+1)
  680  CONTINUE
C---THIS CALCULATION IS THE REVERSE OF THE RECURSION RELATION USED
C  ABOVE
      DO 691 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 691
      DO 690 KK=KCLDS,1,-1
        UFNCLU(I,KK) = UFNCLU(I,KK+1)*ALFAU(I,KK+1)/(ALFA(I,KK+1)*
     1                 TCLU(I,KK))
        DFNCLU(I,KK) = UFNCLU(I,KK)/ALFA(I,KK)
  690  CONTINUE
  691  CONTINUE
      DO 701 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 701
      DO 700 KK=1,KCLDS+1
        UFNTRN(I,KK) = UFNCLU(I,KK)*UFNTRN(I,KK)
        DFNTRN(I,KK) = DFNCLU(I,KK)*DFNTRN(I,KK)
  700  CONTINUE
  701  CONTINUE
C---CASE OF KK=1( FROM THE GROUND TO THE BOTTOM OF THE LOWEST CLOUD)
      DO 720 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 720
        J2=KBTMSW(I,2)
        DO 710 K=J2,LP1
          UFN(I,K) = UFNTRN(I,1)*TTU(I,K)
          DFN(I,K) = DFNTRN(I,1)*TTD(I,K)
  710   CONTINUE
  720  CONTINUE
C---REMAINING LEVELS (IF ANY!)
      DO 760 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 760
      DO 755 KK=2,KCLDS+1
        J1=KTOPSW(I,KK)
        J2=KBTMSW(I,KK+1)
        IF (J1.EQ.1) GO TO 755
        DO 730 K=J2,J1
          UFN(I,K) = UFNTRN(I,KK)*TTU(I,K)
          DFN(I,K) = DFNTRN(I,KK)*TTD(I,K)
  730   CONTINUE
C---FOR THE THICK CLOUDS, THE FLUX DIVERGENCE THROUGH THE CLOUD
C   LAYER IS ASSUMED TO BE CONSTANT. THE FLUX DERIVATIVE IS GIVEN BY
C   TEMPF (FOR THE UPWARD FLUX) AND TEMPG (FOR THE DOWNWARD FLUX).
        J3=KBTMSW(I,KK)
        IF ((J3-J1).GT.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 740 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))
  740     CONTINUE
        ENDIF
  755  CONTINUE
  760  CONTINUE
      DO 770 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 770
      DO 771 K=1,LP1
        DFSWC(I,K) = DFN(I,K)*DFNTOP(I,1)
        UFSWC(I,K) = UFN(I,K)*DFNTOP(I,1)
  771  CONTINUE
  770  CONTINUE
      DO 780 I=MYIS,MYIE
      KCLDS=NCLDS(I)
      IF(KCLDS.EQ.0) GO TO 780
        TMP1(I) = ONE - CCMAX(I)
        GDFVB(I) = TMP1(I)*GDFVB(I)
        GDFNB(I) = TMP1(I)*GDFNB(I)
        GDFVD(I) = TMP1(I)*GDFVD(I) + CCMAX(I)*DFSWC(I,LP1)
  780  CONTINUE
C---NOW OBTAIN FLUXES FOR THE NEAR IR BANDS. THE METHODS ARE THE SAME
C   AS FOR THE VISIBLE BAND, EXCEPT THAT THE REFLECTION AND
C   TRANSMISSION COEFFICIENTS ARE DIFFERENT, AS
C   RAYLEIGH SCATTERING NEED NOT BE CONSIDERED.
C
      DO 1000 N=2,NB
CYH93
        DO 791 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 791
        DO 790 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)
  790   CONTINUE
  791  CONTINUE
CYH93
        IF (N.EQ.2) THEN
C   THE WATER VAPOR TRANSMISSION FUNCTION FOR BAND 2 IS EQUAL TO
C   THAT OF BAND 1 (SAVED AS TTDB1,TTUB1)
          DO 800 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 800
        DO 801 KK=2,LP1
            TTD(I,KK) = TTDB1(I,KK)*TDCO2(I,KK)
  801     CONTINUE
        DO 802 KK=1,L
            TTU(I,KK) = TTUB1(I,KK)*TUCO2(I,KK)
  802     CONTINUE
  800  CONTINUE
        ELSE
          DO 810 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 810
        DO 811 KK=2,LP1
            TTD(I,KK) = EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UD(I,KK)))
     1               * TDCO2(I,KK)
  811     CONTINUE
        DO 812 KK=1,L
            TTU(I,KK) = EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UR(I,KK)))
     1               * TUCO2(I,KK)
  812     CONTINUE
  810  CONTINUE
        ENDIF
C---AT THIS POINT,INCLUDE TTD(1),TTU(LP1), NOTING THAT TTD(1)=1 FOR
C   ALL BANDS, AND THAT TTU(LP1)=TTD(LP1) FOR ALL BANDS.
        DO 820 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 820
          TTU(I,LP1) = TTD(I,LP1)
          TTD(I,1)   = ONE
  820   CONTINUE
C***FOR EXECUTION OF THE CLOUD LOOP, IT IS NECESSARY TO SEPARATE OUT
C   TRANSMISSION FCTNS AT THE TOP AND BOTTOM OF THE CLOUDS, FOR
C   EACH BAND N. THE REQUIRED QUANTITIES ARE:
C      TTD(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
C      TTD(I,KBTMSW(I,K),N)  K RUNS FROM 2 TO NCLDS(I)+1:
C      TTU(I,KTOPSW(I,K),N)  K RUNS FROM 1 TO NCLDS(I)+1:
C      AND INVERSES OF THE ABOVE. THE ABOVE QUANTITIES ARE STORED
C      IN TDCL1,TDCL2,TUCL1,AND DFNTRN,UFNTRN,RESPECTIVELY, AS
C      THEY HAVE MULTIPLE USE IN THE PGM.
C---FOR FIRST CLOUD LAYER (GROUND) TDCL1,TUCL1 ARE KNOWN:
        DO 830 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 830
          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)
  830   CONTINUE
        DO 841 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 841
        DO 840 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))
  840   CONTINUE
  841  CONTINUE
        DO 851 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 851
        DO 850 KK=2,KCLDS+1
          DFNTRN(I,KK) = ONE/TDCL1(I,KK)
          UFNTRN(I,KK) = ONE/TUCL1(I,KK)
  850   CONTINUE
  851  CONTINUE
        DO 861 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 861
        DO 860 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)
  860   CONTINUE
  861  CONTINUE
C***THE FOLLOWING IS THE RECURSION RELATION FOR ALFA: THE REFLECTION
C   COEFFICIENT FOR A SYSTEM INCLUDING THE CLOUD IN QUESTION AND THE
C   FLUX COMING OUT OF THE CLOUD SYSTEM INCLUDING ALL CLOUDS BELOW
C   THE CLOUD IN QUESTION.
        DO 870 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 870
          ALFA (I,1) = CR(I,1)
          ALFAU(I,1) = ZERO
  870   CONTINUE
C---AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER!
        DO 881 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 881
        DO 880 KK=2,KCLDS+1
          ALFAU(I,KK) = TCLU(I,KK-1)*TCLU(I,KK-1)*ALFA(I,KK-1)/(ONE -
     1             TCLD(I,KK-1)*TCLD(I,KK-1)*ALFA(I,KK-1)*CR(I,KK))
          ALFA (I,KK) = ALFAU(I,KK)+CR(I,KK)
  880   CONTINUE
  881  CONTINUE
C     CALCULATE UFN AT CLOUD TOPS AND DFN AT CLOUD BOTTOMS
C---NOTE THAT UFNCLU(I,KCLDS+1) GIVES THE UPWARD FLUX AT THE TOP
C   OF THE HIGHEST REAL CLOUD (IF NCLDS(I)=KCLDS). IT GIVES THE FLUX
C   AT THE TOP OF THE ATMOSPHERE IF NCLDS(I) < KCLDS. IT THE FIRST
C   CASE, TDCL1 EQUALS THE TRANSMISSION FCTN TO THE TOP OF THE
C   HIGHEST CLOUD, AS WE WANT. IN THE SECOND CASE, TDCL1=1, SO UFNCLU
C   EQUALS ALFA. THIS IS ALSO CORRECT.
        DO 890 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 890
          UFNCLU(I,KCLDS+1) = ALFA(I,KCLDS+1)*TDCL1(I,KCLDS+1)
          DFNCLU(I,KCLDS+1) = TDCL1(I,KCLDS+1)
  890   CONTINUE
        DO 901 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 901
        DO 900 KK=KCLDS,1,-1
          UFNCLU(I,KK)=UFNCLU(I,KK+1)*ALFAU(I,KK+1)/(ALFA(I,KK+1)*
     1                 TCLU(I,KK))
          DFNCLU(I,KK) = UFNCLU(I,KK)/ALFA(I,KK)
  900   CONTINUE
  901  CONTINUE
C     NOW OBTAIN DFN AND UFN FOR LEVELS BETWEEN THE CLOUDS
        DO 911 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 911
        DO 910 KK=1,KCLDS+1
          UFNTRN(I,KK) = UFNCLU(I,KK)*UFNTRN(I,KK)
          DFNTRN(I,KK) = DFNCLU(I,KK)*DFNTRN(I,KK)
  910   CONTINUE
  911  CONTINUE
        DO 930 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 930
          J2=KBTMSW(I,2)
          DO 920 K=J2,LP1
            UFN(I,K) = UFNTRN(I,1)*TTU(I,K)
            DFN(I,K) = DFNTRN(I,1)*TTD(I,K)
  920     CONTINUE
  930  CONTINUE
        DO 970  I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 970
        DO 965  KK=2,KCLDS+1
          J1 = KTOPSW(I,KK)
          J2 = KBTMSW(I,KK+1)
          IF (J1.EQ.1) GO TO 965
          DO 940 K=J2,J1
            UFN(I,K) = UFNTRN(I,KK)*TTU(I,K)
            DFN(I,K) = DFNTRN(I,KK)*TTD(I,K)
  940     CONTINUE
          J3 = KBTMSW(I,KK)
          IF ((J3-J1).GT.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 950 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))
  950       CONTINUE
          ENDIF
  965   CONTINUE
  970  CONTINUE
        DO 980 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 980
        DO 981 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)
  981   CONTINUE
  980  CONTINUE
        DO 990 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 990
          GDFND(I) = GDFND(I) + CCMAX(I)*DFN(I,LP1)*DFNTOP(I,N)
  990   CONTINUE
 1000 CONTINUE
      DO 1100 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 1100
      DO 1101 K=1,LP1
        DFSWC(I,K) = TMP1(I)*DFSWL(I,K) + CCMAX(I)*DFSWC(I,K)
        UFSWC(I,K) = TMP1(I)*UFSWL(I,K) + CCMAX(I)*UFSWC(I,K)
 1101 CONTINUE
 1100 CONTINUE
      DO 1200 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 1200
        DO 1201 KK=1,LP1
        FSWC(I,KK) = UFSWC(I,KK)-DFSWC(I,KK)
 1201 CONTINUE
 1200 CONTINUE
      DO 1250 I=MYIS,MYIE
        KCLDS=NCLDS(I)
        IF(KCLDS.EQ.0) GO TO 1250
        DO 1251 KK=1, L
        HSWC(I,KK) = RADCON*(FSWC(I,KK+1)-FSWC(I,KK))/DP(I,KK)
 1251 CONTINUE
 1250 CONTINUE
      RETURN
      END
