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) !===> ********************************************************* ! --- SWR91SIB --- MODIFIED FROM SWR89-BAND12....YUTAI HOU ! -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 ! INCLUDE HPCON,PARMC CHANGED TO HCON,RDPARM ! 6 Q..... VARIABLE NAMES RESTORED TO ORIGINAL 7,8 CHAR ! CHANGE O3DIFF,DIFFCC TO O3DIFCTR,DIFFCTR ! --- SWR91SIB --- MODIFIED BY Y. HOU FEB 93 ! INPUTS 12 BANDS CLD REFLECTTANCE AND TRANSMITTANCE ! CRR,CTT TO REPLACE CIRAB,CIRRF,CUVRF !===> ********************************************************* INCLUDE "parmeta.f90" INCLUDE "COMM_HCON.f90" INCLUDE "rdparm.f90" INCLUDE "mpp.h" #include "sp.h" ! PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE: ! L = NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL ! NB IS A SHORTWAVE PARAMETER; OTHER QUANTITIES ARE DERIVED ! FROM THE ABOVE PARAMETERS. ! --- VARIABLES AS IN ARGUMENT LIST DIMENSION & FSWC (IDIM1:IDIM2,LP1), HSWC (IDIM1:IDIM2,LP1), & CRR (IDIM1:IDIM2,NB,LP1) & , FSWL (IDIM1:IDIM2,LP1), HSWL (IDIM1:IDIM2,LP1), & CTT (IDIM1:IDIM2,NB,LP1) & , UFSWC (IDIM1:IDIM2,LP1), DFSWC (IDIM1:IDIM2,LP1) & , UFSWL (IDIM1:IDIM2,LP1), DFSWL (IDIM1:IDIM2,LP1) & , PRESS (IDIM1:IDIM2,LP1), RH2O (IDIM1:IDIM2,L), & QO3 (IDIM1:IDIM2,L) & , CAMT (IDIM1:IDIM2,LP1), KTOPSW(IDIM1:IDIM2,LP1), & KBTMSW(IDIM1:IDIM2,LP1) & , COSZRO(IDIM1:IDIM2), TAUDAR(IDIM1:IDIM2), & NCLDS (IDIM1:IDIM2) & , ALVB (IDIM1:IDIM2), ALNB (IDIM1:IDIM2), & ALVD (IDIM1:IDIM2), ALND (IDIM1:IDIM2) & , GDFVB (IDIM1:IDIM2), GDFNB (IDIM1:IDIM2), GDFVD (IDIM1:IDIM2), & GDFND (IDIM1:IDIM2) ! --- LOCAL VARIABLES DIMENSION & PP (IDIM1:IDIM2,LP1), DP (IDIM1:IDIM2,LP1), & PR2 (IDIM1:IDIM2,LP1) & , DU (IDIM1:IDIM2,LP1), DUCO2 (IDIM1:IDIM2,LP1), & DUO3 (IDIM1:IDIM2,LP1) & , FF (IDIM1:IDIM2,LP1), FFCO2 (IDIM1:IDIM2,LP1), & FFO3 (IDIM1:IDIM2,LP1) & , RRAY (IDIM1:IDIM2), DFNTOP(IDIM1:IDIM2,NB), & SECZ (IDIM1:IDIM2) & , REFL (IDIM1:IDIM2), TMP1 (IDIM1:IDIM2), & REFL2 (IDIM1:IDIM2) & , CCMAX (IDIM1:IDIM2), XAMT (IDIM1:IDIM2,LP1) DIMENSION & UD (IDIM1:IDIM2,LP1), UR (IDIM1:IDIM2,LP1) & , UCO2 (IDIM1:IDIM2,LLP2), UDCO2 (IDIM1:IDIM2,LP1), & URCO2 (IDIM1:IDIM2,LP1) & , UO3 (IDIM1:IDIM2,LLP2), UDO3 (IDIM1:IDIM2,LP1), & URO3 (IDIM1:IDIM2,LP1) & , TCO2 (IDIM1:IDIM2,LLP2), TDCO2 (IDIM1:IDIM2,LP1), & TUCO2 (IDIM1:IDIM2,LP1) & , TO3 (IDIM1:IDIM2,LLP2), TDO3 (IDIM1:IDIM2,LP1), & TUO3 (IDIM1:IDIM2,LP1) DIMENSION & DFN (IDIM1:IDIM2,LP1), UFN (IDIM1:IDIM2,LP1), & CR (IDIM1:IDIM2,LP1) & , TTD (IDIM1:IDIM2,LP1), TTU (IDIM1:IDIM2,LP1), & CT (IDIM1:IDIM2,LP1) & , PPTOP (IDIM1:IDIM2,LP1), DPCLD (IDIM1:IDIM2,LP1) ! --- EQUIVALENCED LOCAL VARIABLES DIMENSION & TTUB1 (IDIM1:IDIM2,LP1), TUCL1 (IDIM1:IDIM2,LP1) & , TTDB1 (IDIM1:IDIM2,LP1), TDCL1 (IDIM1:IDIM2,LP1), & TDCL2 (IDIM1:IDIM2,LP1) & , UFNTRN(IDIM1:IDIM2,LP1), UFNCLU(IDIM1:IDIM2,LP1), & TCLU (IDIM1:IDIM2,LP1) & , DFNTRN(IDIM1:IDIM2,LP1), DFNCLU(IDIM1:IDIM2,LP1), & TCLD (IDIM1:IDIM2,LP1) & , ALFA (IDIM1:IDIM2,LP1), ALFAU (IDIM1:IDIM2,LP1) EQUIVALENCE & (UDO3,UO3(IDIM1,1),DFNCLU), (URO3,UO3(IDIM1,LP2), UFNCLU) & , (UDCO2,UCO2(IDIM1,1),TCLD), (URCO2,UCO2(IDIM1,LP2), TCLU) & , (TDO3 ,TO3(IDIM1,1),DFNTRN),(TUO3,TO3(IDIM1,LP2), UFNTRN) & , (TDCO2,TCO2(IDIM1,1) ),(TUCO2,TCO2(IDIM1,LP2) ) & , (FF , ALFA ), (FFCO2 , ALFAU ), (FFO3 , TTDB1 ) & , (DU , TTUB1), (DUCO2 , TUCL1 ), (DUO3 , TDCL1 ) & , (PR2 , TDCL2) !---COMMON FOR LOCAL DATA VARIABLES--- COMMON /SWRSAV/ ABCFF(NB),PWTS(NB),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 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 END DO DO 110 K=1,LM1 DO 110 I=MYIS,MYIE PP(I,K+1) = HAF*(PRESS(I,K+1)+PRESS(I,K)) 110 END DO 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 END DO DO 130 K=1,L DO 130 I=MYIS,MYIE PR2(I,K) = PR2(I,K)*TMP1(I) 130 END DO ! 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) !if (mype == 2) write(*,*) "TESTES2", IP,N,DFNTOP(IP,N),SSOLAR,H69766E5,COSZRO(IP),TAUDAR(IP),PWTS(N) 140 END DO ! EXECUTE THE LACIS-HANSEN REFLECTIVITY PARAMETERIZATION ! 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)/ & (ONE-ALVD(I)*RRAYAV) 150 END DO 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)/ & (ONE-ALVD(I)*0.093) 155 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 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) !if (mype == 2) write(*,*) "TESTEWW1",I,K,DUO3 (I,K),GINV,CFO3,QO3(I,K),DP(I,K) 160 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 200 IP=MYIS,MYIE UD (IP,1) = ZERO UDCO2(IP,1) = ZERO UDO3 (IP,1) = ZERO 200 END DO 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) !if (mype == 2) write(*,*) "TESTESWR1",I,K,UDO3 (I,K),UDO3 (I,K-1),DUO3 (I,K-1),SECZ(I) 210 END DO DO 220 IP=MYIS,MYIE UR (IP,LP1) = UD (IP,LP1) URCO2(IP,LP1) = UDCO2(IP,LP1) URO3 (IP,LP1) = UDO3 (IP,LP1) !if (mype == 2) write(*,*) "TESTESWR0",IP,URO3(IP,LP1),UDO3 (IP,LP1) 220 END DO 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 !if (mype == 2) write(*,*) "TESTEWR0",IP,K,URO3(IP,K),URO3 (IP,K+1),DUO3 (IP,K),O3DIFCTR,DIFFCTR 230 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 240 K=1,LL DO 240 I=MYIS,MYIE TCO2(I,K+1)=ONE-TWO*(H235M3*EXP(HP26*LOG(UCO2(I,K+1)+H129M2)) & -H75826M4) 240 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 250 K=1,LL DO 250 I=MYIS,MYIE TO3(I,K+1)=ONE-TWO*UO3(I,K+1)* & (H1P082*EXP(HMP805*LOG(ONE+H1386E2*UO3(I,K+1)))+ & H658M2/(ONE+HTEMP*UO3(I,K+1)*UO3(I,K+1)*UO3(I,K+1))+ & H2118M2/(ONE+UO3(I,K+1)*(H42M2+H323M4*UO3(I,K+1)))) !if (mype == 2) write(*,*) "TESTEWR1",TO3(I,K+1),ONE,TWO,UO3(I,K+1),H1P082,EXP(HMP805*LOG(ONE+H1386E2*UO3(I,K+1))),HMP805,H1386E2,H658M2,HTEMP,H2118M2,H42M2,H323M4 250 END DO ! START FREQUENCY LOOP (ON N) HERE !--- 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) !if (mype == 2) write(*,*) "TESTEWR", I,K,TTD(I,K+1),TDO3(I,K+1),EXP(HM1EZ*MIN(FIFTY,ABCFF(1)*UD(I,K+1))),HM1EZ,FIFTY,ABCFF(1),UD(I,K+1),MIN(FIFTY,ABCFF(1)*UD(I,K+1)) UFN(I,K) = TTU(I,K)*TUO3(I,K) 260 END DO DO 270 I=MYIS,MYIE DFN(I,1) = ONE UFN(I,LP1) = DFN(I,LP1) 270 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 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) !if (mype == 2) write(*,*) "TESTES1", DFN(I,K),DFNTOP(I,1) 280 END DO DO 285 I=MYIS,MYIE 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 !if (mype == 2) write(*,*) "TESTESW", GDFVB(I),DFSWL(I,LP1),SECZ(I),EXP(-0.15746*SECZ(I)) 285 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 350 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 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 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 300 K=1,L DO 300 I=MYIS,MYIE DFN(I,K+1)=EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UD(I,K+1))) & *TDCO2(I,K+1) UFN(I,K)=EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UR(I,K))) & *TUCO2(I,K) 300 END DO ENDIF !---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 310 I=MYIS,MYIE DFN(I,1) = ONE UFN(I,LP1) = DFN(I,LP1) 310 END DO ! SCALE THE PREVIOUSLY COMPUTED FLUXES BY THE FLUX AT THE TOP ! 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 END DO DO 330 I=MYIS,MYIE GDFNB(I) = GDFNB(I) + DFN(I,LP1)*DFNTOP(I,N) 330 END DO 350 END DO DO 360 K=1,LP1 DO 360 I=MYIS,MYIE FSWL(I,K) = UFSWL(I,K)-DFSWL(I,K) 360 END DO 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 END DO !---END OF FREQUENCY LOOP (OVER N) ! CALCULATE CLOUDY SKY SW FLUX KCLDS=NCLDS(MYIS) DO 400 I=MYIS1,MYIE KCLDS=MAX(NCLDS(I),KCLDS) 400 END DO 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 END DO DO 420 K=1,L DO 420 I=MYIS,MYIE HSWC(I,K) = HSWL(I,K) 420 END DO !******************************************************************* IF (KCLDS == 0) RETURN !******************************************************************* DO 430 K=1,LP1 DO 430 I=MYIS,MYIE XAMT(I,K) = CAMT(I,K) 430 END DO DO 470 I=MYIS,MYIE NNCLDS = NCLDS(I) CCMAX(I) = ZERO IF (NNCLDS <= 0) GO TO 470 CCMAX(I) = ONE DO 450 K=1,NNCLDS CCMAX(I) = CCMAX(I) * (ONE - CAMT(I,K+1)) 450 END DO CCMAX(I) = ONE - CCMAX(I) IF (CCMAX(I) > ZERO) THEN DO 460 K=1,NNCLDS XAMT(I,K+1) = CAMT(I,K+1)/CCMAX(I) 460 END DO END IF 470 END DO DO 480 K=1,LP1 DO 480 I=MYIS,MYIE FF (I,K) = DIFFCTR FFCO2(I,K) = DIFFCTR FFO3 (I,K) = O3DIFCTR 480 END DO 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 END DO DO 500 I=MYIS,MYIE RRAY(I) = HP219/(ONE+HP816*COSZRO(I)) REFL(I) = RRAY(I) + (ONE-RRAY(I))*(ONE-RRAYAV)*ALVD(I)/ & (ONE-ALVD(I)*RRAYAV) 500 END DO DO 510 IP=MYIS,MYIE UD (IP,1) = ZERO UDCO2(IP,1) = ZERO UDO3 (IP,1) = ZERO 510 END DO 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 END DO DO 530 IP=MYIS,MYIE UR (IP,LP1) = UD (IP,LP1) URCO2(IP,LP1) = UDCO2(IP,LP1) URO3 (IP,LP1) = UDO3 (IP,LP1) 530 END DO 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 END DO 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)) & -H75826M4) 550 END DO DO 560 K=1,LL DO 560 I=MYIS,MYIE TO3(I,K+1)=ONE-TWO*UO3(I,K+1)* & (H1P082*EXP(HMP805*LOG(ONE+H1386E2*UO3(I,K+1)))+ & H658M2/(ONE+HTEMP*UO3(I,K+1)*UO3(I,K+1)*UO3(I,K+1))+ & H2118M2/(ONE+UO3(I,K+1)*(H42M2+H323M4*UO3(I,K+1)))) 560 END DO !******************************************************************** !---THE FIRST CLOUD IS THE GROUND; ITS PROPERTIES ARE GIVEN ! BY REFL (THE TRANSMISSION (0) IS IRRELEVANT FOR NOW!). !******************************************************************** DO 570 I=MYIS,MYIE CR(I,1) = REFL(I) 570 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 581 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 581 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 591 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 591 DO 590 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))) ENDIF 590 END DO 591 END DO 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 END DO DO 610 I=MYIS,MYIE TTD(I,1) = ONE TTU(I,LP1) = TTD(I,LP1) 610 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 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 END DO DO 631 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 631 END DO !---COMPUTE INVERSES DO 641 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 641 DO 640 KK=2,KCLDS DFNTRN(I,KK) = ONE/TDCL1(I,KK) UFNTRN(I,KK) = ONE/TUCL1(I,KK) 640 END DO 641 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 651 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 651 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 660 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 660 ALFA (I,1)=CR(I,1) ALFAU(I,1)=ZERO 660 END DO !---AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER! DO 671 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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)/ & (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 END DO 671 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 680 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO !---THIS CALCULATION IS THE REVERSE OF THE RECURSION RELATION USED ! ABOVE DO 691 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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)* & TCLU(I,KK)) DFNCLU(I,KK) = UFNCLU(I,KK)/ALFA(I,KK) 690 END DO 691 END DO DO 701 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 701 END DO !---CASE OF KK=1( FROM THE GROUND TO THE BOTTOM OF THE LOWEST CLOUD) DO 720 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 720 END DO !---REMAINING LEVELS (IF ANY!) DO 760 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 760 DO 755 KK=2,KCLDS+1 J1=KTOPSW(I,KK) J2=KBTMSW(I,KK+1) IF (J1 == 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 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 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 END DO ENDIF 755 END DO 760 END DO DO 770 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 770 END DO DO 780 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 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 1000 N=2,NB ! H93 DO 791 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 791 END DO ! H93 IF (N == 2) THEN ! THE WATER VAPOR TRANSMISSION FUNCTION FOR BAND 2 IS EQUAL TO ! THAT OF BAND 1 (SAVED AS TTDB1,TTUB1) DO 800 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 800 DO 801 KK=2,LP1 TTD(I,KK) = TTDB1(I,KK)*TDCO2(I,KK) 801 END DO DO 802 KK=1,L TTU(I,KK) = TTUB1(I,KK)*TUCO2(I,KK) 802 END DO 800 END DO ELSE DO 810 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 810 DO 811 KK=2,LP1 TTD(I,KK) = EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UD(I,KK))) & * TDCO2(I,KK) 811 END DO DO 812 KK=1,L TTU(I,KK) = EXP(HM1EZ*MIN(FIFTY,ABCFF(N)*UR(I,KK))) & * TUCO2(I,KK) 812 END DO 810 END DO ENDIF !---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 820 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 820 TTU(I,LP1) = TTD(I,LP1) TTD(I,1) = ONE 820 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 830 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO DO 841 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 841 END DO DO 851 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 851 END DO DO 861 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 861 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 870 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 870 ALFA (I,1) = CR(I,1) ALFAU(I,1) = ZERO 870 END DO !---AGAIN,EXCESSIVE CALCULATIONS-MAY BE CHANGED LATER! DO 881 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 - & 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 END DO 881 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 890 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO DO 901 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 901 DO 900 KK=KCLDS,1,-1 !--- Ferrier, 6/17/02: Emergency change to eliminate problematic ! features of unrealistically small cloud amounts DENOM=ALFA(I,KK+1)*TCLU(I,KK) IF (DENOM > 1.E-15) THEN UFNCLU(I,KK)=UFNCLU(I,KK+1)*ALFAU(I,KK+1)/DENOM ELSE ! print * ! &, ' xnum=',UFNCLU(I,KK+1)*ALFAU(I,KK+1) ! &, ' xden=',DENOM ! &, ' ALFA(I,KK+1)=',ALFA(I,KK+1) ! &, ' TCLU(I,KK)=',TCLU(I,KK) ! &, ' UFNCLU(I,KK+1)=',UFNCLU(I,KK+1) ! &, ' ALFAU(I,KK+1)=',ALFAU(I,KK+1) ! &, ' UFNCLU(I,KK)=',UFNCLU(I,KK) ! &, ' ALFA(I,KK)=',ALFA(I,KK) ! &, ' I=',I ! &, ' KK+1=',KK+1 UFNCLU(I,KK)=0. ENDIF DFNCLU(I,KK) = UFNCLU(I,KK)/ALFA(I,KK) 900 END DO 901 END DO ! NOW OBTAIN DFN AND UFN FOR LEVELS BETWEEN THE CLOUDS DO 911 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 911 END DO DO 930 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 930 END DO DO 970 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 970 DO 965 KK=2,KCLDS+1 J1 = KTOPSW(I,KK) J2 = KBTMSW(I,KK+1) IF (J1 == 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 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 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 END DO ENDIF 965 END DO 970 END DO DO 980 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 980 END DO DO 990 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 990 GDFND(I) = GDFND(I) + CCMAX(I)*DFN(I,LP1)*DFNTOP(I,N) 990 END DO 1000 END DO DO 1100 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 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 END DO 1100 END DO DO 1200 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 1200 DO 1201 KK=1,LP1 FSWC(I,KK) = UFSWC(I,KK)-DFSWC(I,KK) 1201 END DO 1200 END DO DO 1250 I=MYIS,MYIE KCLDS=NCLDS(I) IF(KCLDS == 0) GO TO 1250 DO 1251 KK=1, L HSWC(I,KK) = RADCON*(FSWC(I,KK+1)-FSWC(I,KK))/DP(I,KK) 1251 END DO 1250 END DO RETURN END SUBROUTINE SWR93