C ***************************************************************** C SUBROUTINE FST88 IS THE MAIN COMPUTATION MODULE OF THE C LONG-WAVE RADIATION CODE. IN IT ALL "EMISSIVITY" CALCULATIONS, C INCLUDING CALLS TO TABLE LOOKUP SUBROUTINES. ALSO,AFTER CALLING C SUBROUTINE "SPA88", FINAL COMBINED HEATING RATES AND GROUND C FLUX ARE OBTAINED. C ***************************************************************** C INPUTS: C BETINW,BETAWD,AB15WD BDWIDE C BETAD,BO3RND,AO3RND BANDTA C CLDFAC CLDCOM C QH2O,P,DELP2,DELP,T,VAR1,VAR2, KDACOM C VAR3,VAR4,CNTVAL KDACOM C TOTVO2,TOTO3,TOTPHI,EMPL,EMX1 KDACOM C TPHIO3,EMX2 KDACOM C TEMP,PRESS RADISW C NCLDS,KTOP,KBTM,CAMT RADISW C IND,INDX2,KMAXV,SOURCE,DSRCE TABCOM C SKC1R,SKC3R,KMAXVM,NREP1,NREP2 TABCOM C NST1,NST2,NRP1,NRP2 TABCOM C CO2NBL,CO21 TFCOM C CO2SP1,CO2SP2 TFCOM C OUTPUTS: C HEATRA,GRNFLX,TOPFLX LWOUT C C CALLED BY : RADMN OR MAIN PGM C CALLS : CLO88,E1E288,E3V88,SPA88,NLTE C C PASSED VARIABLES: C IN E3V88: C EMD = E3 FUNCTION FOR H2O LINES (0-560,1200-2200 CM-1) C COMPUTED IN E3V88 C TPL = TEMPERATURE INPUT FOR E3 CALCULATION IN E3V88 C EMPL = H2O AMOUNT,INPUT FOR E3 CALCULATION IN E3V88 C (COMPUTED IN LWR88; STORED IN KDACOM.H) C IN E1E288: C E1CTS1 = E1 FUNCTION FOR THE (I+1)TH LEVEL USING THE C TEMPERATURE OF THE ITH DATA LEVEL,COMPUTED OVER C THE FREQUENCY RANGE 0-560,1200-2200 CM-1. (E1CTS1- C E1CTW1) IS USED IN OBTAINING THE FLUX AT THE TOP C IN THE 0-160,1200-2200 CM-1 RANGE (FLX1E1). C E1CTS2 = E1 FUNCTION FOR THE ITH LEVEL, USING THE TEMP. OF C THE ITH DATA LEVEL,COMPUTED OVER THE FREQUENCY RANGE C 0-560,1200-2200 CM-1. (E1CTS2-E1CTW2) IS ALSO USED C IN OBTAINING THE FLUX AT THE TOP IN THE 0-160,. C 1200-2200 CM-1 RANGE. C E1FLX = E1 FCTN. FOR THE ITH LEVEL,USING THE TEMPERATURE AT C THE TOP OF THE ATMOSPHERE. COMPUTED OVER THE FREQ. C RANGE 0-560,1200-2200 CM-1. USED FOR Q(APPROX) TERM. C (IN COMMON BLOCK TFCOM) C E1CTW1 = LIKE E1CTS1,BUT COMPUTED OVER THE 160-560 CM-1 RANGE C AND USED FOR Q(APPROX,CTS) CALCULATION C E1CTW2 = LIKE E1CTS2,BUT COMPUTED OVER THE 160-560 CM-1 RANGE C AND USED FOR Q(APPROX,CTS) CALCULATION C FXO = TEMPERATURE INDEX USED FOR E1 FUNCTION AND ALSO C USED FOR SOURCE FUNCTION CALC. IN FST88. C DT = TEMP. DIFF.BETWEEN MODEL TEMPS. AND TEMPS. AT C TABULAR VALUES OF E1 AND SOURCE FCTNS. USED IN C FST88 AND IN E1 FUNCTION CALC. C FXOE2 = TEMPERATURE INDEX USED FOR E2 FUNCTION C DTE2 = TEMP. DIFF. BETWEEN MODEL TEMP. AND TEMPS. AT C TABULAR VALUES OF E2 FUNCTION. SUBROUTINE FST88(HEATRA,GRNFLX,TOPFLX, 1 QH2O,PRESS,P,DELP,DELP2,TEMP,T, 2 CLDFAC,NCLDS,KTOP,KBTM,CAMT, 3 CO21,CO2NBL,CO2SP1,CO2SP2, 4 VAR1,VAR2,VAR3,VAR4,CNTVAL, 5 TOTO3,TPHIO3,TOTPHI,TOTVO2, 6 EMX1,EMX2,EMPL) C COMMON/PHYCON/AMOLWT,CSUBP,DIFFCTR,G,GRAVDR,O3DIFCTR,P0, * P0XZP2,P0XZP8,P0X2,RADCON,RGAS,RGASSP,SECPDA COMMON/PHYCON/RATCO2MW,RATH2OMW COMMON/PHYCON/RADCON1 COMMON/PHYCON/GINV,P0INV,GP0INV save /PHYCON/ COMMON/HCON/HUNDRED,HNINETY,SIXTY,FIFTY,TEN,EIGHT,FIVE, * FOUR,THREE,TWO,ONE,HAF,QUARTR,ZERO COMMON/HCON/H83E26,H71E26,H1E15,H1E13,H1E11,H1E8,H4E5, * H165E5,H5725E4,H488E4,H1E4,H24E3,H20788E3, * H2075E3,H1224E3,H5E2,H3082E2,H3E2,H2945E2, * H23E2,H15E2,H35E1,H3P6,H181E1,H18E1,H2P9,H2P8, * H2P5,H1P8,H1P4387,H1P4,H1P25892,HP8,HP518, * HP369,HP1 COMMON/HCON/H44871M2,H559M3,H1M3,H987M4,H285M4,H1M4, * H6938M5,H394M5,H37412M5,H1439M5,H128M5,H1M5, * H7M6,H4999M6,H25452M6,H1M6,H391M7,H1174M7, * H8725M8,H327M8,H257M8,H1M8,H23M10,H14M10, * H11M10,H1M10,H83M11,H82M11,H8M11,H77M11, * H72M11,H53M11,H48M11,H44M11,H42M11,H37M11, * H35M11,H32M11,H3M11,H28M11,H24M11,H23M11, * H2M11,H18M11,H15M11,H14M11,H114M11,H11M11, * H1M11,H96M12,H93M12,H77M12,H74M12,H65M12, * H62M12,H6M12,H45M12,H44M12,H4M12,H38M12, * H37M12,H3M12,H29M12,H28M12,H24M12,H21M12, * H16M12,H14M12,H12M12,H8M13,H46M13,H36M13, * H135M13,H12M13,H1M13,H3M14,H15M14,H14M14, * H1M17,H1M18,H1M19,H1M20,H1M21,H1M22,H1M23, * H1M24,H26M30,H14M30,H25M31,H21M31,H12M31, * H9M32,H55M32,H45M32,H4M33,H62M34,H1M60 COMMON/HCON/HMP575,HM13EZ,HM19EZ,HM1E1,HM181E1,HM1E2 COMMON/HCON/H1E6,H2E6,H1M2,HMP66667,HM6666M2,HP166666, * H41666M2,HMP5,HM2M2,H29316E2,H1226E1,H3116E1, * H9P94,HP6,H625M2,HP228,HP60241,HM1797E1, * H8121E1,H2E2,HM1EZ,H26E2,H44194M2,H1P41819 COMMON/HCON/HP219,HP144,HP816,H69766E5,H235M3,HP26, * H129M2,H75826M4,H1P082,HP805,H1386E2, * H658M2,H1036E2,H2118M2,H42M2,H323M4, * H67390E2,HP3795,HP5048,H102M5,H451M6 COMMON/HCON/H16E1,HM161E1,H161E1,H3M3,H101M16, * HM1597E1,H25E2,HP118666,H15M5,H3P5,H18E3, * H6P08108,HMP805,HP602409,HP526315, * H28571M2,H1M16 COMMON/HCON/H3M4 COMMON/HCON/HM8E1 COMMON/HCON/H28E1 save /HCON/ C----------------------------------------------------------------------- INCLUDE "parmeta" INCLUDE "mpp.h" #include "sp.h" C----------------------------------------------------------------------- C PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE: C IMAX = NO. POINTS ALONG THE LAT. CIRCLE USED IN CALCS. C L = NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL C***NOTE: THE USER NORMALLY WILL MODIFY ONLY THE IMAX AND L PARAMETERS C NBLW = NO. FREQ. BANDS FOR APPROX COMPUTATIONS. SEE C BANDTA FOR DEFINITION C NBLX = NO. FREQ BANDS FOR APPROX CTS COMPUTATIONS C NBLY = NO. FREQ. BANDS FOR EXACT CTS COMPUTATIONS. SEE C BDCOMB FOR DEFINITION C INLTE = NO. LEVELS USED FOR NLTE CALCS. C NNLTE = INDEX NO. OF FREQ. BAND IN NLTE CALCS. C NB,KO2 ARE SHORTWAVE PARAMETERS; OTHER QUANTITIES ARE DERIVED C FROM THE ABOVE PARAMETERS. PARAMETER (L=LM) PARAMETER (IMAX=IM,NCOL=IMAX) PARAMETER (NBLW=163,NBLX=47,NBLY=15) PARAMETER (NBLM=NBLY-1) PARAMETER (LP1=L+1,LP2=L+2,LP3=L+3) PARAMETER (LM1=L-1,LM2=L-2,LM3=L-3) PARAMETER (LL=2*L,LLP1=LL+1,LLP2=LL+2,LLP3=LL+3) PARAMETER (LLM1=LL-1,LLM2=LL-2,LLM3=LL-3) PARAMETER (LP1M=LP1*LP1,LP1M1=LP1M-1) PARAMETER (LP1V=LP1*(1+2*L/2)) PARAMETER (LP121=LP1*NBLY) PARAMETER (LL3P=3*L+2) PARAMETER (NB=12) PARAMETER (INLTE=3,INLTEP=INLTE+1,NNLTE=56) PARAMETER (LP1I=IMAX*LP1,LLP1I=IMAX*LLP1,LL3PI=IMAX*LL3P) PARAMETER (NB1=NB-1) PARAMETER (KO2=12) PARAMETER (KO21=KO2+1,KO2M=KO2-1) C PARAMETER SETTINGS FOR THE LONGWAVE AND SHORTWAVE RADIATION CODE: C IMAX = NO. POINTS SENT TO RADFS C L = NO. VERTICAL LEVELS (ALSO LAYERS) IN MODEL C***NOTE: THE USER NORMALLY WILL MODIFY ONLY THE IMAX AND L PARAMETERS C NBLW = NO. FREQ. BANDS FOR APPROX COMPUTATIONS. SEE C BANDTA FOR DEFINITION C NBLX = NO. FREQ BANDS FOR APPROX CTS COMPUTATIONS C NBLY = NO. FREQ. BANDS FOR EXACT CTS COMPUTATIONS. SEE C BDCOMB FOR DEFINITION C INLTE = NO. LEVELS USED FOR NLTE CALCS. C NNLTE = INDEX NO. OF FREQ. BAND IN NLTE CALCS. C NB,KO2 ARE SHORTWAVE PARAMETERS; OTHER QUANTITIES ARE DERIVED C FROM THE ABOVE PARAMETERS. C COMMON BLOCK BANDTA CONTAINS RANDOM BAND PARAMETERS FOR THE LW C CALCULATIONS USING 10 CM-1 WIDE BANDS.THE 15 UM CO2 COMPLEX C IS 2 BANDS,560-670 AND 670-800 CM-1. OZONE COEFFICIENTS ARE C IN 3 BANDS,670-800 (14.1 UM),990-1070 AND 1070-1200 (9.6 UM). C THE (NBLW) BANDS NOW INCLUDE: C 56 BANDS, 10 CM-1 WIDE 0 - 560 CM-1 C 2 BANDS, 15 UM COMPLEX 560 - 670 CM-1 C 670 - 800 CM-1 C 3 "CONTINUUM" BANDS 800 - 900 CM-1 C 900 - 990 CM-1 C 1070 - 1200 CM-1 C 1 BAND FOR 9.6 UM BAND 990 - 1070 CM-1 C 100 BANDS, 10 CM-1 WIDE 1200 - 2200 CM-1 C 1 BAND FOR 4.3 UM SRC 2270 - 2380 CM-1 C THUS NBLW PRESENTLY EQUALS 163 C ALL BANDS ARE ARRANGED IN ORDER OF INCREASING WAVENUMBER C C ARNDM = RANDOM "A" PARAMETER FOR (NBLW) BANDS C BRNDM = RANDOM "B" PARAMETER FOR (NBLW) BANDS C BETAD = CONTINUUM COEFFICIENTS FOR (NBLW) BANDS C AP,BP = CAPPHI COEFFICIENTS FOR (NBLW) BANDS C ATP,BTP = CAPPSI COEFFICIENTS FOR (NBLW) BANDS C BANDLO = LOWEST FREQUENCY IN EACH OF (NBLW) FREQ. BANDS C BANDHI = HIGHEST FREQUENCY IN EACH OF (NBLW) FREQ. BANDS C AO3RND = RANDOM "A" PARAMETER FOR OZONE IN (3) OZONE C BANDS C BO3RND = RANDOM "B" PARAMETER FOR OZONE IN (3) OZONE C BANDS C AB15 = THE PRODUCT ARNDM*BRNDM FOR THE TWO BANDS C REPRESENTING THE 15 UM BAND COMPLEX OF CO2 C DATA FOR ARNDM,BRNDM,AP,BP,ATP,BTP,AO3RND,BO3RND ARE OBTAINED BY C USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS ARE FROM C ROBERTS (1976). COMMON / BANDTA / ARNDM(NBLW),BRNDM(NBLW),BETAD(NBLW),AP(NBLW), 1 BP(NBLW),ATP(NBLW),BTP(NBLW),BANDLO(NBLW), 2 BANDHI(NBLW),AO3RND(3),BO3RND(3),AB15(2) C C COMMON BLOCK BDWIDE CONTAINS RANDOM BAND PARAMETERS FOR SPECIFIC C WIDE BANDS. AT PRESENT,THE INFORMATION CONSISTS OF 1) RANDOM C MODEL PARAMETERS FOR THE 15 UM BAND,560-800 CM-1; 2) THE C CONTINUUM COEFFICIENT FOR THE 800-990,1070-1200 CM-1 BAND C SPECIFICALLY: C AWIDE = RANDOM "A" PARAMETER FOR BAND C BWIDE = RANDOM "B" PARAMETER FOR BAND C BETAWD = CONTINUUM COEFFICIENTS FOR BAND C APWD,BPWD = CAPPHI COEFFICIENTS FOR BAND C ATPWD,BTPWD = CAPPSI COEFFICIENTS FOR BAND C BDLOWD = LOWEST FREQUENCY IN EACH FREQ BAND C BDHIWD = HIGHEST FREQUENCY IN EACH FREQ BAND C AB15WD = THE PRODUCT ARNDM*BRNDM FOR THE ONE BAND C REPRESENTING THE 15 UM BAND COMPLEX OF CO2 C BETINW = CONT.COEFFICIENT FOR A SPECIFIED WIDE C FREQ.BAND (800-990 AND 1070-1200 CM-1). C SKO2D = 1./BETINW, USED IN SPA88 FOR CONT. COEFFS C SKC1R = BETAWD/BETINW, USED FOR CONT. COEFF. FOR C 15 UM BAND IN FST88 C SKO3R = RATIO OF CONT. COEFF. FOR 9.9 UM BAND TO C BETINW, USED FOR 9.6 UM CONT COEFF IN FST88 C DATA FOR AWIDE,BWIDE,APWD,BPWD,ATPWD,BTPWD,AO3WD,BO3WD ARE C OBTAINED BY USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS C ARE FROM ROBERTS (1976). COMMON /BDWIDE/ AWIDE,BWIDE,BETAWD, 1 APWD,BPWD,ATPWD,BTPWD, 2 BDLOWD,BDHIWD,BETINW, 3 AB15WD,SKO2D,SKC1R,SKO3R save /BDWIDE/ C C COMMON BLOCK BDCOMB CONTAINS RANDOM BAND PARAMETERS FOR THE LW C CALCULATIONS USING COMBINED WIDE FREQUENCY BANDS BETWEEN 160 AND C 1200 CM-1,AS WELL AS THE 2270-2380 BAND FOR SOURCE CALC. C BANDS 1-8: COMBINED WIDE FREQUENCY BANDS FOR 160-560 CM-1 C BANDS 9-14: FREQUENCY BANDS,AS IN BANDTA (NARROW BANDS) C FOR 560-1200 CM-1 C BAND 15: FREQUENCY BAND 2270-2380 CM-1,USED FOR SOURCE C CALCULATION ONLY C THUS NBLY PRESENTLY EQUALS 15 C C BANDS ARE ARRANGED IN ORDER OF INCREASING WAVENUMBER C ACOMB = RANDOM "A" PARAMETER FOR (NBLY) BANDS C BCOMB = RANDOM "B" PARAMETER FOR (NBLY) BANDS C BETACM = CONTINUUM COEFFICIENTS FOR (NBLY) BANDS C APCM,BPCM = CAPPHI COEFFICIENTS FOR (NBLY) BANDS C ATPCM,BTPCM = CAPPSI COEFFICIENTS FOR (NBLY) BANDS C BDLOCM = LOWEST FREQUENCY IN EACH OF (NBLY) FREQ. BANDS C BDHICM = HIGHEST FREQUENCY IN EACH OF (NBLY) FREQ. BANDS C AO3CM = RANDOM "A" PARAMETER FOR OZONE IN (3) OZONE C BANDS C BO3CM = RANDOM "B" PARAMETER FOR OZONE IN (3) OZONE C BANDS C AB15CM = THE PRODUCT ARNDM*BRNDM FOR THE TWO BANDS C REPRESENTING THE 15 UM BAND COMPLEX OF CO2 C BETINC = CONT.COEFFICIENT FOR A SPECIFIED WIDE C FREQ.BAND (800-990 AND 1070-1200 CM-1). C IBAND = INDEX NO OF THE 40 WIDE BANDS USED IN C COMBINED WIDE BAND CALCULATIONS. IN OTHER C WORDS,INDEX TELLING WHICH OF THE 40 WIDE C BANDS BETWEEN 160-560 CM-1 ARE INCLUDED IN C EACH OF THE FIRST 8 COMBINED WIDE BANDS C DATA FOR ACOMB,BCOMB,APCM,BPCM,ATPCM,BTPCM,AO3CM,BO3CM ARE C OBTAINED BY USING THE AFGL 1982 CATALOG. CONTINUUM COEFFICIENTS C ARE FROM ROBERTS (1976). IBAND INDEX VALUES ARE OBTAINED BY C EXPERIMENTATION. COMMON / BDCOMB / IBAND(40),ACOMB(NBLY),BCOMB(NBLY), 1 BETACM(NBLY),APCM(NBLY),BPCM(NBLY),ATPCM(NBLY), 2 BTPCM(NBLY),BDLOCM(NBLY),BDHICM(NBLY),BETINC, 3 AO3CM(3),BO3CM(3),AB15CM(2) save / BDCOMB / C C COMMON BLOCK TABCOM CONTAINS QUANTITIES PRECOMPUTED IN SUBROUTINE C TABLE FOR USE IN THE LONGWAVE RADIATION PROGRAM: C EM1 = E1 FUNCTION, EVALUATED OVER THE 0-560 AND C 1200-2200 CM-1 INTERVALS C EM1WDE = E1 FUNCTION, EVALUATED OVER THE 160-560 CM-1 C INTERVAL C TABLE1 = E2 FUNCTION, EVALUATED OVER THE 0-560 AND C 1200-2200 CM-1 INTERVALS C TABLE2 = TEMPERATURE DERIVATIVE OF TABLE1 C TABLE3 = MASS DERIVATIVE OF TABLE1 C EM3 = E3 FUNCTION, EVALUATED OVER THE 0-560 AND C 1200-2200 CM-1 INTERVALS C SOURCE = PLANCK FUNCTION, EVALUATED AT SPECIFIED TEMPS. FOR C BANDS USED IN CTS CALCULATIONS C DSRCE = TEMPERATURE DERIVATIVE OF SOURCE C IND = INDEX, WITH VALUE IND(I)=I. USED IN FST88 C INDX2 = INDEX VALUES USED IN OBTAINING "LOWER TRIANGLE" C ELEMENTS OF AVEPHI,ETC.,IN FST88 C KMAXV = INDEX VALUES USED IN OBTAINING "UPPER TRIANGLE" C ELEMENTS OF AVEPHI,ETC.,IN FST88 C KMAXVM = KMAXV(L),USED FOR DO LOOP INDICES C COMMON /TABCOM/ IND(IMAX),INDX2(LP1V),KMAXV(LP1), 1 KMAXVM COMMON/TABCOM/EM1(28,180),EM1WDE(28,180),TABLE1(28,180), 1 TABLE2(28,180),TABLE3(28,180),EM3(28,180),SOURCE(28,NBLY), 2 DSRCE(28,NBLY) save /TABCOM/ C DIMENSION QH2O(IDIM1:IDIM2,LP1),PRESS(IDIM1:IDIM2,LP1) DIMENSION P(IDIM1:IDIM2,LP1),DELP(IDIM1:IDIM2,L), & DELP2(IDIM1:IDIM2,L),TEMP(IDIM1:IDIM2,LP1) DIMENSION T(IDIM1:IDIM2,LP1),CLDFAC(IDIM1:IDIM2,LP1,LP1), & CAMT(IDIM1:IDIM2,LP1) DIMENSION NCLDS(IDIM1:IDIM2),KTOP(IDIM1:IDIM2,LP1), & KBTM(IDIM1:IDIM2,LP1) DIMENSION CO21(IDIM1:IDIM2,LP1,LP1),CO2NBL(IDIM1:IDIM2,L) DIMENSION CO2SP1(IDIM1:IDIM2,LP1),CO2SP2(IDIM1:IDIM2,LP1) DIMENSION VAR1(IDIM1:IDIM2,L),VAR2(IDIM1:IDIM2,L), & VAR3(IDIM1:IDIM2,L),VAR4(IDIM1:IDIM2,L) DIMENSION CNTVAL(IDIM1:IDIM2,LP1) DIMENSION HEATRA(IDIM1:IDIM2,L),GRNFLX(IDIM1:IDIM2), & TOPFLX(IDIM1:IDIM2) DIMENSION GXCTS(IDIM1:IDIM2),FLX1E1(IDIM1:IDIM2) DIMENSION AVEPHI(IDIM1:IDIM2,LP1),EMISS(IDIM1:IDIM2,LP1), & EMISSB(IDIM1:IDIM2,LP1) C DIMENSION TOTO3(IDIM1:IDIM2,LP1),TPHIO3(IDIM1:IDIM2,LP1), & TOTPHI(IDIM1:IDIM2,LP1) DIMENSION TOTVO2(IDIM1:IDIM2,LP1),EMX1(IDIM1:IDIM2), & EMX2(IDIM1:IDIM2),EMPL(IDIM1:IDIM2,LLP1) C DIMENSION EXCTS(IDIM1:IDIM2,L),CTSO3(IDIM1:IDIM2,L), & CTS(IDIM1:IDIM2,L),E1FLX(IDIM1:IDIM2,LP1) DIMENSION CO2SP(IDIM1:IDIM2,LP1),TO3SPC(IDIM1:IDIM2,L), & TO3SP(IDIM1:IDIM2,LP1) DIMENSION OSS(IDIM1:IDIM2,LP1),CSS(IDIM1:IDIM2,LP1), & SS1(IDIM1:IDIM2,LP1),SS2(IDIM1:IDIM2,LP1), 1 TC(IDIM1:IDIM2,LP1),DTC(IDIM1:IDIM2,LP1) DIMENSION SORC(IDIM1:IDIM2,LP1,NBLY),CSOUR(IDIM1:IDIM2,LP1) CCC DIMENSION AVVO2(IDIM1:IDIM2,LP1),HEATEM(IDIM1:IDIM2,LP1), 1 OVER1D(IDIM1:IDIM2,LP1), 1 TO31D(IDIM1:IDIM2,LP1),CONT1D(IDIM1:IDIM2,LP1), 2 AVMO3(IDIM1:IDIM2,LP1),AVPHO3(IDIM1:IDIM2,LP1), 2 C(IDIM1:IDIM2,LLP1),C2(IDIM1:IDIM2,LLP1) DIMENSION ITOP(IDIM1:IDIM2),IBOT(IDIM1:IDIM2), & INDTC(IDIM1:IDIM2) DIMENSION 4 DELPTC(IDIM1:IDIM2),PTOP(IDIM1:IDIM2),PBOT(IDIM1:IDIM2), & FTOP(IDIM1:IDIM2), 5 FBOT(IDIM1:IDIM2) ,EMSPEC(IDIM1:IDIM2,2) C---DIMENSION OF VARIABLES EQUIVALENCED TO THOSE IN VTEMP--- DIMENSION VTMP3(IDIM1:IDIM2,LP1),DSORC(IDIM1:IDIM2,LP1) DIMENSION ALP(IDIM1:IDIM2,LLP1),CSUB(IDIM1:IDIM2,LLP1), & CSUB2(IDIM1:IDIM2,LLP1) DIMENSION FAC1(IDIM1:IDIM2,LP1) DIMENSION DELPR1(IDIM1:IDIM2,LP1),DELPR2(IDIM1:IDIM2,LP1) DIMENSION EMISDG(IDIM1:IDIM2,LP1),CONTDG(IDIM1:IDIM2,LP1), & TO3DG(IDIM1:IDIM2,LP1) DIMENSION FLXNET(IDIM1:IDIM2,LP1) DIMENSION IXO(IDIM1:IDIM2,LP1) DIMENSION VSUM1(IDIM1:IDIM2,LP1) DIMENSION FLXTHK(IDIM1:IDIM2,LP1) DIMENSION Z1(IDIM1:IDIM2,LP1) C---DIMENSION OF VARIABLES PASSED TO OTHER SUBROUTINES--- C (AND NOT FOUND IN COMMON BLOCKS) DIMENSION E1CTS1(IDIM1:IDIM2,LP1),E1CTS2(IDIM1:IDIM2,L) DIMENSION E1CTW1(IDIM1:IDIM2,LP1),E1CTW2(IDIM1:IDIM2,L) DIMENSION EMD(IDIM1:IDIM2,LLP1),TPL(IDIM1:IDIM2,LLP1) C IT IS POSSIBLE TO EQUIVALENCE EMD,TPL TO THE ABOVE VARIABLES, C AS THEY GET CALLED AT DIFFERENT TIMES DIMENSION FXO(IDIM1:IDIM2,LP1),DT(IDIM1:IDIM2,LP1) DIMENSION FXOE2(IDIM1:IDIM2,LP1),DTE2(IDIM1:IDIM2,LP1) DIMENSION FXOSP(IDIM1:IDIM2,2),DTSP(IDIM1:IDIM2,2) C C DIMENSION OF LOCAL VARIABLES DIMENSION RLOG(IDIM1:IDIM2,L),FLX(IDIM1:IDIM2,LP1) DIMENSION TOTEVV(IDIM1:IDIM2,LP1),CNTTAU(IDIM1:IDIM2,LP1) C EQUIVALENCE (ALP,C,CSUB),(CSUB2,C2) EQUIVALENCE (FAC1,DSORC,OVER1D,DELPR2,FLXNET) EQUIVALENCE (DELPR1,HEATEM) EQUIVALENCE (IXO,AVVO2,FLXTHK,TO3DG) EQUIVALENCE (Z1,AVMO3,CONTDG) EQUIVALENCE (EMISDG,VSUM1,AVPHO3) EQUIVALENCE (EMD(IDIM1,1),E1CTS1(IDIM1,1)), & (EMD(IDIM1,LP2),E1CTS2(IDIM1,1)) EQUIVALENCE (TPL(IDIM1,1),E1CTW1(IDIM1,1)), & (TPL(IDIM1,LP2),E1CTW2(IDIM1,1)) c C C FIRST SECTION IS TABLE LOOKUP FOR SOURCE FUNCTION AND C DERIVATIVE (B AND DB/DT).ALSO,THE NLTE CO2 SOURCE FUNCTION C IS OBTAINED C C---IN CALCS. BELOW, DECREMENTING THE INDEX BY 9 C ACCOUNTS FOR THE TABLES BEGINNING AT T=100K. C AT T=100K. DO 101 K=1,LP1 DO 101 I=MYIS,MYIE C---TEMP. INDICES FOR E1,SOURCE VTMP3(I,K)=AINT(TEMP(I,K)*HP1) FXO(I,K)=VTMP3(I,K)-9. DT(I,K)=TEMP(I,K)-TEN*VTMP3(I,K) C---INTEGER INDEX FOR SOURCE (USED IMMEDIATELY) C wne IXO(I,K)=FXO(I,K) IXO(I,K)=max(FXO(I,K), 1.0) 101 CONTINUE DO 103 k=1,L DO 103 I=MYIS,MYIE C---TEMP. INDICES FOR E2 (KP=1 LAYER NOT USED IN FLUX CALCULATIONS) VTMP3(I,K)=AINT(T(I,K+1)*HP1) FXOE2(I,K)=VTMP3(I,K)-9. DTE2(I,K)=T(I,K+1)-TEN*VTMP3(I,K) 103 CONTINUE C---SPECIAL CASE TO HANDLE KP=LP1 LAYER AND SPECIAL E2 CALCS. DO 105 I=MYIS,MYIE FXOE2(I,LP1)=FXO(I,L) DTE2(I,LP1)=DT(I,L) FXOSP(I,1)=FXOE2(I,LM1) FXOSP(I,2)=FXO(I,LM1) DTSP(I,1)=DTE2(I,LM1) DTSP(I,2)=DT(I,LM1) 105 CONTINUE C C---SOURCE FUNCTION FOR COMBINED BAND 1 DO 4114 I=MYIS,MYIE DO 4114 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),1) DSORC(I,K)=DSRCE(IXO(I,K),1) 4114 CONTINUE DO 4112 K=1,LP1 DO 4112 I=MYIS,MYIE SORC(I,K,1)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4112 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 2 DO 4214 I=MYIS,MYIE DO 4214 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),2) DSORC(I,K)=DSRCE(IXO(I,K),2) 4214 CONTINUE DO 4212 K=1,LP1 DO 4212 I=MYIS,MYIE SORC(I,K,2)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4212 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 3 DO 4314 I=MYIS,MYIE DO 4314 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),3) DSORC(I,K)=DSRCE(IXO(I,K),3) 4314 CONTINUE DO 4312 K=1,LP1 DO 4312 I=MYIS,MYIE SORC(I,K,3)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4312 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 4 DO 4414 I=MYIS,MYIE DO 4414 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),4) DSORC(I,K)=DSRCE(IXO(I,K),4) 4414 CONTINUE DO 4412 K=1,LP1 DO 4412 I=MYIS,MYIE SORC(I,K,4)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4412 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 5 DO 4514 I=MYIS,MYIE DO 4514 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),5) DSORC(I,K)=DSRCE(IXO(I,K),5) 4514 CONTINUE DO 4512 K=1,LP1 DO 4512 I=MYIS,MYIE SORC(I,K,5)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4512 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 6 DO 4614 I=MYIS,MYIE DO 4614 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),6) DSORC(I,K)=DSRCE(IXO(I,K),6) 4614 CONTINUE DO 4612 K=1,LP1 DO 4612 I=MYIS,MYIE SORC(I,K,6)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4612 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 7 DO 4714 I=MYIS,MYIE DO 4714 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),7) DSORC(I,K)=DSRCE(IXO(I,K),7) 4714 CONTINUE DO 4712 K=1,LP1 DO 4712 I=MYIS,MYIE SORC(I,K,7)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4712 CONTINUE C---SOURCE FUNCTION FOR COMBINED BAND 8 DO 4814 I=MYIS,MYIE DO 4814 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),8) DSORC(I,K)=DSRCE(IXO(I,K),8) 4814 CONTINUE DO 4812 K=1,LP1 DO 4812 I=MYIS,MYIE SORC(I,K,8)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4812 CONTINUE C---SOURCE FUNCTION FOR BAND 9 (560-670 CM-1) DO 4914 I=MYIS,MYIE DO 4914 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),9) DSORC(I,K)=DSRCE(IXO(I,K),9) 4914 CONTINUE DO 4912 K=1,LP1 DO 4912 I=MYIS,MYIE SORC(I,K,9)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 4912 CONTINUE C---SOURCE FUNCTION FOR BAND 10 (670-800 CM-1) DO 5014 I=MYIS,MYIE DO 5014 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),10) DSORC(I,K)=DSRCE(IXO(I,K),10) 5014 CONTINUE DO 5012 K=1,LP1 DO 5012 I=MYIS,MYIE SORC(I,K,10)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 5012 CONTINUE C---SOURCE FUNCTION FOR BAND 11 (800-900 CM-1) DO 5114 I=MYIS,MYIE DO 5114 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),11) DSORC(I,K)=DSRCE(IXO(I,K),11) 5114 CONTINUE DO 5112 K=1,LP1 DO 5112 I=MYIS,MYIE SORC(I,K,11)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 5112 CONTINUE C---SOURCE FUNCTION FOR BAND 12 (900-990 CM-1) DO 5214 I=MYIS,MYIE DO 5214 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),12) DSORC(I,K)=DSRCE(IXO(I,K),12) 5214 CONTINUE DO 5212 K=1,LP1 DO 5212 I=MYIS,MYIE SORC(I,K,12)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 5212 CONTINUE C---SOURCE FUNCTION FOR BAND 13 (990-1070 CM-1) DO 5314 I=MYIS,MYIE DO 5314 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),13) DSORC(I,K)=DSRCE(IXO(I,K),13) 5314 CONTINUE DO 5312 K=1,LP1 DO 5312 I=MYIS,MYIE SORC(I,K,13)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 5312 CONTINUE C---SOURCE FUNCTION FOR BAND 14 (1070-1200 CM-1) DO 5414 I=MYIS,MYIE DO 5414 K=1,LP1 VTMP3(I,K)=SOURCE(IXO(I,K),14) DSORC(I,K)=DSRCE(IXO(I,K),14) 5414 CONTINUE DO 5412 K=1,LP1 DO 5412 I=MYIS,MYIE SORC(I,K,14)=VTMP3(I,K)+DT(I,K)*DSORC(I,K) 5412 CONTINUE C C THE FOLLOWING SUBROUTINE OBTAINS NLTE SOURCE FUNCTION FOR CO2 C C C CALL NLTE C C C---OBTAIN SPECIAL SOURCE FUNCTIONS FOR THE 15 UM BAND (CSOUR) C AND THE WINDOW REGION (SS1) DO 131 K=1,LP1 DO 131 I=MYIS,MYIE SS1(I,K)=SORC(I,K,11)+SORC(I,K,12)+SORC(I,K,14) 131 CONTINUE DO 143 K=1,LP1 DO 143 I=MYIS,MYIE CSOUR(I,K)=SORC(I,K,9)+SORC(I,K,10) 143 CONTINUE C C---COMPUTE TEMP**4 (TC) AND VERTICAL TEMPERATURE DIFFERENCES C (OSS,CSS,SS2,DTC). ALL THESE WILL BE USED LATER IN FLUX COMPUTA- C TIONS. C DO 901 K=1,LP1 DO 901 I=MYIS,MYIE TC(I,K)=(TEMP(I,K)*TEMP(I,K))**2 c if(mype.eq.13.and.i.eq.40) then c print*,'i,k,temp(i,k)=',i,k,temp(i,k) c endif 901 CONTINUE DO 903 K=1,L DO 903 I=MYIS,MYIE OSS(I,K+1)=SORC(I,K+1,13)-SORC(I,K,13) CSS(I,K+1)=CSOUR(I,K+1)-CSOUR(I,K) DTC(I,K+1)=TC(I,K+1)-TC(I,K) SS2(I,K+1)=SS1(I,K+1)-SS1(I,K) 903 CONTINUE C C C---THE FOLLOWIMG IS A DRASTIC REWRITE OF THE RADIATION CODE TO C (LARGELY) ELIMINATE THREE-DIMENSIONAL ARRAYS. THE CODE WORKS C ON THE FOLLOWING PRINCIPLES: C C LET K = FIXED FLUX LEVEL, KP = VARYING FLUX LEVEL C THEN FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K)) C OVER ALL KP'S, FROM 1 TO LP1. C C WE CAN BREAK DOWN THE CALCULATIONS FOR ALL K'S AS FOLLOWS: C C FOR ALL K'S K=1 TO LP1: C FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K)) (1) C OVER ALL KP'S, FROM K+1 TO LP1 C AND C FOR KP FROM K+1 TO LP1: C FLUX(KP) = DELTAB(K)*TAU(K,KP) (2) C C NOW IF TAU(K,KP)=TAU(KP,K) (SYMMETRICAL ARRAYS) C WE CAN COMPUTE A 1-DIMENSIONAL ARRAY TAU1D(KP) FROM C K+1 TO LP1, EACH TIME K IS INCREMENTED. C EQUATIONS (1) AND (2) THEN BECOME: C C TAU1D(KP) = (VALUES FOR TAU(KP,K) AT THE PARTICULAR K) C FLUX(K) = SUM OVER KP : (DELTAB(KP)*TAU1D(KP)) (3) C FLUX(KP) = DELTAB(K)*TAU1D(KP) (4) C C THE TERMS FOR TAU (K,K) AND OTHER SPECIAL TERMS (FOR C NEARBY LAYERS) MUST, OF COURSE, BE HANDLED SEPARATELY, AND C WITH CARE. C C COMPUTE "UPPER TRIANGLE" TRANSMISSION FUNCTIONS FOR C THE 9.6 UM BAND (TO3SP) AND THE 15 UM BAND (OVER1D). ALSO, C THE C STAGE 1...COMPUTE O3 ,OVER TRANSMISSION FCTNS AND AVEPHI C---DO K=1 CALCULATION (FROM FLUX LAYER KK TO THE TOP) SEPARATELY C AS VECTORIZATION IS IMPROVED,AND OZONE CTS TRANSMISSIVITY C MAY BE EXTRACTED HERE. DO 3021 K=1,L DO 3021 I=MYIS,MYIE AVEPHI(I,K)=TOTPHI(I,K+1) 3021 CONTINUE C---IN ORDER TO PROPERLY EVALUATE EMISS INTEGRATED OVER THE (LP1) C LAYER, A SPECIAL EVALUATION OF EMISS IS DONE. THIS REQUIRES C A SPECIAL COMPUTATION OF AVEPHI, AND IT IS STORED IN THE C (OTHERWISE VACANT) LP1'TH POSITION C DO 803 I=MYIS,MYIE AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I) 803 CONTINUE C COMPUTE FLUXES FOR K=1 CALL E1E290(E1CTS1,E1CTS2,E1FLX,E1CTW1,E1CTW2,EMISS, 1 FXO,DT,FXOE2,DTE2,AVEPHI,TEMP,T) DO 302 K=1,L DO 302 I=MYIS,MYIE FAC1(I,K)=BO3RND(2)*TPHIO3(I,K+1)/TOTO3(I,K+1) TO3SPC(I,K)=HAF*(FAC1(I,K)* 1 (SQRT(ONE+(FOUR*AO3RND(2)*TOTO3(I,K+1))/FAC1(I,K))-ONE)) C FOR K=1, TO3SP IS USED INSTEAD OF TO31D (THEY ARE EQUAL IN THIS C CASE); TO3SP IS PASSED TO SPA90, WHILE TO31D IS A WORK-ARRAY. TO3SP(I,K)=EXP(HM1EZ*(TO3SPC(I,K)+SKO3R*TOTVO2(I,K+1))) OVER1D(I,K)=EXP(HM1EZ*(SQRT(AB15WD*TOTPHI(I,K+1))+ 1 SKC1R*TOTVO2(I,K+1))) C---BECAUSE ALL CONTINUUM TRANSMISSIVITIES ARE OBTAINED FROM THE C 2-D QUANTITY CNTTAU (AND ITS RECIPROCAL TOTEVV) WE STORE BOTH C OF THESE HERE. FOR K=1, CONT1D EQUALS CNTTAU CNTTAU(I,K)=EXP(HM1EZ*TOTVO2(I,K+1)) TOTEVV(I,K)=1./CNTTAU(I,K) 302 CONTINUE DO 3022 K=1,L DO 3022 I=MYIS,MYIE CO2SP(I,K+1)=OVER1D(I,K)*CO21(I,1,K+1) 3022 CONTINUE DO 3023 K=1,L DO 3023 I=MYIS,MYIE CO21(I,K+1,1)=CO21(I,K+1,1)*OVER1D(I,K) 3023 CONTINUE C---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION DO 1808 I=MYIS,MYIE RLOG(I,1)=OVER1D(I,1)*CO2NBL(I,1) 1808 CONTINUE C---THE TERMS WHEN KP=1 FOR ALL K ARE THE PHOTON EXCHANGE WITH C THE TOP OF THE ATMOSPHERE, AND ARE OBTAINED DIFFERENTLY THAN C THE OTHER CALCULATIONS DO 305 K=2,LP1 DO 305 I=MYIS,MYIE FLX(I,K)= (TC(I,1)*E1FLX(I,K) 1 +SS1(I,1)*CNTTAU(I,K-1) 2 +SORC(I,1,13)*TO3SP(I,K-1) 3 +CSOUR(I,1)*CO2SP(I,K)) 4 *CLDFAC(I,1,K) 305 CONTINUE DO 307 I=MYIS,MYIE FLX(I,1)= TC(I,1)*E1FLX(I,1)+SS1(I,1)+SORC(I,1,13) 1 +CSOUR(I,1) 307 CONTINUE C---THE KP TERMS FOR K=1... DO 303 KP=2,LP1 DO 303 I=MYIS,MYIE FLX(I,1)=FLX(I,1)+(OSS(I,KP)*TO3SP(I,KP-1) 1 +SS2(I,KP)*CNTTAU(I,KP-1) 2 +CSS(I,KP)*CO21(I,KP,1) 3 +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,1) 303 CONTINUE C SUBROUTINE SPA88 IS CALLED TO OBTAIN EXACT CTS FOR WATER C CO2 AND O3, AND APPROXIMATE CTS CO2 AND O3 CALCULATIONS. C CALL SPA88(EXCTS,CTSO3,GXCTS,SORC,CSOUR, 1 CLDFAC,TEMP,PRESS,VAR1,VAR2, 2 P,DELP,DELP2,TOTVO2,TO3SP,TO3SPC, 3 CO2SP1,CO2SP2,CO2SP) C C THIS SECTION COMPUTES THE EMISSIVITY CTS HEATING RATES FOR 2 C EMISSIVITY BANDS: THE 0-160,1200-2200 CM-1 BAND AND THE 800- C 990,1070-1200 CM-1 BAND. THE REMAINING CTS COMTRIBUTIONS ARE C CONTAINED IN CTSO3, COMPUTED IN SPA88. C DO 998 I=MYIS,MYIE VTMP3(I,1)=1. 998 CONTINUE DO 999 K=1,L DO 999 I=MYIS,MYIE VTMP3(I,K+1)=CNTTAU(I,K)*CLDFAC(I,K+1,1) 999 CONTINUE DO 1001 K=1,L DO 1001 I=MYIS,MYIE CTS(I,K)=RADCON*DELP(I,K)*(TC(I,K)* 1 (E1CTW2(I,K)*CLDFAC(I,K+1,1)-E1CTW1(I,K)*CLDFAC(I,K,1)) + 2 SS1(I,K)*(VTMP3(I,K+1)-VTMP3(I,K))) 1001 CONTINUE C DO 1011 K=1,L DO 1011 I=MYIS,MYIE VTMP3(I,K)=TC(I,K)*(CLDFAC(I,K,1)*(E1CTS1(I,K)-E1CTW1(I,K)) - 1 CLDFAC(I,K+1,1)*(E1CTS2(I,K)-E1CTW2(I,K))) 1011 CONTINUE DO 1012 I=MYIS,MYIE FLX1E1(I)=TC(I,LP1)*CLDFAC(I,LP1,1)* 1 (E1CTS1(I,LP1)-E1CTW1(I,LP1)) c if(mype.eq.13.and.i.eq.40) then c print*,'i,lp1=',i,lp1 c print*,'tc(i,lp1)=',tc(i,lp1) c print*,'cldfac(i,lp1,1)=',cldfac(i,lp1,1) c print*,'e1cts1(i,lp1)=',e1cts1(i,lp1) c print*,'e1ctw1(i,lp1)=',e1ctw1(i,lp1) c endif 1012 CONTINUE DO 1014 K=1,L DO 1013 I=MYIS,MYIE FLX1E1(I)=FLX1E1(I)+VTMP3(I,K) 1013 CONTINUE 1014 CONTINUE C C---NOW REPEAT FLUX CALCULATIONS FOR THE K=2..LM1 CASES. C CALCULATIONS FOR FLUX LEVEL L AND LP1 ARE DONE SEPARATELY, AS ALL C EMISSIVITY AND CO2 CALCULATIONS ARE SPECIAL CASES OR NEARBY LAYERS. C DO 321 K=2,LM1 KLEN=K C DO 3218 KK=1,LP1-K DO 3218 I=MYIS,MYIE AVEPHI(I,KK+K-1)=TOTPHI(I,KK+K)-TOTPHI(I,K) 3218 CONTINUE DO 1803 I=MYIS,MYIE AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I) 1803 CONTINUE C---COMPUTE EMISSIVITY FLUXES (E2) FOR THIS CASE. NOTE THAT C WE HAVE OMITTED THE NEARBY LATER CASE (EMISS(I,K,K)) AS WELL C AS ALL CASES WITH K=L OR LP1. BUT THESE CASES HAVE ALWAYS C BEEN HANDLED AS SPECIAL CASES, SO WE MAY AS WELL COMPUTE C THEIR FLUXES SEPARASTELY. C CALL E290(EMISSB,EMISS,AVEPHI,KLEN,FXOE2,DTE2) DO 322 KK=1,LP1-K DO 322 I=MYIS,MYIE AVMO3(I,KK+K-1)=TOTO3(I,KK+K)-TOTO3(I,K) AVPHO3(I,KK+K-1)=TPHIO3(I,KK+K)-TPHIO3(I,K) AVVO2(I,KK+K-1)=TOTVO2(I,KK+K)-TOTVO2(I,K) CONT1D(I,KK+K-1)=CNTTAU(I,KK+K-1)*TOTEVV(I,K-1) 322 CONTINUE C DO 3221 KK=1,LP1-K DO 3221 I=MYIS,MYIE FAC1(I,K+KK-1)=BO3RND(2)*AVPHO3(I,K+KK-1)/AVMO3(I,K+KK-1) VTMP3(I,K+KK-1)=HAF*(FAC1(I,K+KK-1)* 1 (SQRT(ONE+(FOUR*AO3RND(2)*AVMO3(I,K+KK-1))/ 2 FAC1(I,K+KK-1))-ONE)) TO31D(I,K+KK-1)=EXP(HM1EZ*(VTMP3(I,K+KK-1) 1 +SKO3R*AVVO2(I,K+KK-1))) OVER1D(I,K+KK-1)=EXP(HM1EZ*(SQRT(AB15WD*AVEPHI(I,K+KK-1))+ 1 SKC1R*AVVO2(I,K+KK-1))) CO21(I,K+KK,K)=OVER1D(I,K+KK-1)*CO21(I,K+KK,K) 3221 CONTINUE DO 3223 KP=K+1,LP1 DO 3223 I=MYIS,MYIE CO21(I,K,KP)=OVER1D(I,KP-1)*CO21(I,K,KP) 3223 CONTINUE C---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION DO 1804 I=MYIS,MYIE RLOG(I,K)=OVER1D(I,K)*CO2NBL(I,K) 1804 CONTINUE C---THE KP TERMS FOR ARBIRRARY K.. DO 3423 KP=K+1,LP1 DO 3423 I=MYIS,MYIE FLX(I,K)=FLX(I,K)+(OSS(I,KP)*TO31D(I,KP-1) 1 +SS2(I,KP)*CONT1D(I,KP-1) 2 +CSS(I,KP)*CO21(I,KP,K) 3 +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,K) 3423 CONTINUE DO 3425 KP=K+1,LP1 DO 3425 I=MYIS,MYIE FLX(I,KP)=FLX(I,KP)+(OSS(I,K)*TO31D(I,KP-1) 1 +SS2(I,K)*CONT1D(I,KP-1) 2 +CSS(I,K)*CO21(I,K,KP) 3 +DTC(I,K)*EMISSB(I,KP-1))*CLDFAC(I,K,KP) 3425 CONTINUE 321 CONTINUE C C NOW DO K=L CASE. SINCE THE KP LOOP IS LENGTH 1, MANY SIMPLIFI- C CATIONS OCCUR. ALSO, THE CO2 QUANTITIES (AS WELL AS THE EMISS C QUANTITIES) ARE COMPUTED IN THE NBL SEDCTION; THEREFORE, WE WANT C ONLY OVER,TO3 AND CONT1D (OVER(I,L),TO31D(I,L) AND CONT1D(I,L) C ACCORDING TO THE NOTATION. THUS NO CALL IS MADE TO THE E290 C SUBROUTINE. C THE THIRD SECTION CALCULATES BOUNDARY LAYER AND NEARBY LAYER C CORRECTIONS TO THE TRANSMISSION FUNCTIONS OBTAINED ABOVE. METHODS C ARE GIVEN IN REF. (4). C THE FOLLOWING RATIOS ARE USED IN VARIOUS NBL CALCULATIONS: C C THE REMAINING CALCULATIONS ARE FOR : C 1) THE (K,K) TERMS, K=2,LM1; C 2) THE (L,L) TERM C 3) THE (L,LP1) TERM C 4) THE (LP1,L) TERM C 5) THE (LP1,LP1) TERM. C EACH IS UNIQUELY HANDLED; DIFFERENT FLUX TERMS ARE COMPUTED C DIFFERENTLY C C C FOURTH SECTION OBTAINS WATER TRANSMISSION FUNCTIONS C USED IN Q(APPROX) CALCULATIONS AND ALSO MAKES NBL CORRECTIONS: C 1) EMISS (I,J) IS THE TRANSMISSION FUNCTION MATRIX OBTAINED C BY CALLING SUBROUTINE E1E288; C 2) "NEARBY LAYER" CORRECTIONS (EMISS(I,I)) ARE OBTAINED C USING SUBROUTINE E3V88; C 3) SPECIAL VALUES AT THE SURFACE (EMISS(L,LP1),EMISS(LP1,L), C EMISS(LP1,LP1)) ARE CALCULATED. C C C OBTAIN ARGUMENTS FOR E1E288 AND E3V88: C DO 821 I=MYIS,MYIE TPL(I,1)=TEMP(I,L) TPL(I,LP1)=HAF*(T(I,LP1)+TEMP(I,L)) TPL(I,LLP1)=HAF*(T(I,L)+TEMP(I,L)) 821 CONTINUE DO 823 K=2,L DO 823 I=MYIS,MYIE TPL(I,K)=T(I,K) TPL(I,K+L)=T(I,K) 823 CONTINUE C C---E2 FUNCTIONS ARE REQUIRED IN THE NBL CALCULATIONS FOR 2 CASES, C DENOTED (IN OLD CODE) AS (L,LP1) AND (LP1,LP1) DO 833 I=MYIS,MYIE AVEPHI(I,1)=VAR2(I,L) AVEPHI(I,2)=VAR2(I,L)+EMPL(I,L) 833 CONTINUE CALL E2SPEC(EMISS,AVEPHI,FXOSP,DTSP) C C CALL E3V88 FOR NBL H2O TRANSMISSIVITIES CALL E3V88(EMD,TPL,EMPL) C C COMPUTE NEARBY LAYER AND SPECIAL-CASE TRANSMISSIVITIES FOR EMISS C USING METHODS FOR H2O GIVEN IN REF. (4) DO 851 K=2,L DO 851 I=MYIS,MYIE EMISDG(I,K)=EMD(I,K+L)+EMD(I,K) 851 CONTINUE C C NOTE THAT EMX1/2 (PRESSURE SCALED PATHS) ARE NOW COMPUTED IN C LWR88 DO 861 I=MYIS,MYIE EMSPEC(I,1)=(EMD(I,1)*EMPL(I,1)-EMD(I,LP1)*EMPL(I,LP1))/ 1 EMX1(I) + QUARTR*(EMISS(I,1)+EMISS(I,2)) EMISDG(I,LP1)=TWO*EMD(I,LP1) EMSPEC(I,2)=TWO*(EMD(I,1)*EMPL(I,1)-EMD(I,LLP1)*EMPL(I,LLP1))/ * EMX2(I) 861 CONTINUE DO 331 I=MYIS,MYIE FAC1(I,L)=BO3RND(2)*VAR4(I,L)/VAR3(I,L) VTMP3(I,L)=HAF*(FAC1(I,L)* 1 (SQRT(ONE+(FOUR*AO3RND(2)*VAR3(I,L))/FAC1(I,L))-ONE)) TO31D(I,L)=EXP(HM1EZ*(VTMP3(I,L)+SKO3R*CNTVAL(I,L))) OVER1D(I,L)=EXP(HM1EZ*(SQRT(AB15WD*VAR2(I,L))+ 1 SKC1R*CNTVAL(I,L))) CONT1D(I,L)=CNTTAU(I,L)*TOTEVV(I,LM1) RLOG(I,L)=OVER1D(I,L)*CO2NBL(I,L) 331 CONTINUE DO 618 K=1,L DO 618 I=MYIS,MYIE RLOG(I,K)=LOG(RLOG(I,K)) 618 CONTINUE DO 601 K=1,LM1 DO 601 I=MYIS,MYIE DELPR1(I,K+1)=DELP(I,K+1)*(PRESS(I,K+1)-P(I,K+1)) ALP(I,LP1+K-1)=-SQRT(DELPR1(I,K+1))*RLOG(I,K+1) 601 CONTINUE DO 603 K=1,L DO 603 I=MYIS,MYIE DELPR2(I,K+1)=DELP(I,K)*(P(I,K+1)-PRESS(I,K)) ALP(I,K)=-SQRT(DELPR2(I,K+1))*RLOG(I,K) 603 CONTINUE DO 625 I=MYIS,MYIE ALP(I,LL)=-RLOG(I,L) ALP(I,LLP1)=-RLOG(I,L)*SQRT(DELP(I,L)*(P(I,LP1)-PRESS(I,LM1))) 625 CONTINUE C THE FIRST COMPUTATION IS FOR THE 15 UM BAND,WITH THE C FOR THE COMBINED H2O AND CO2 TRANSMISSION FUNCTION. C C PERFORM NBL COMPUTATIONS FOR THE 15 UM BAND C***THE STATEMENT FUNCTION SF IN PREV. VERSIONS IS NOW EXPLICITLY C EVALUATED. DO 631 K=1,LLP1 DO 631 I=MYIS,MYIE C(I,K)=ALP(I,K)*(HMP66667+ALP(I,K)*(QUARTR+ALP(I,K)*HM6666M2)) 631 CONTINUE DO 641 I=MYIS,MYIE CO21(I,LP1,LP1)=ONE+C(I,L) CO21(I,LP1,L)=ONE+(DELP2(I,L)*C(I,LL)-(PRESS(I,L)-P(I,L))* 1 C(I,LLM1))/(P(I,LP1)-PRESS(I,L)) CO21(I,L,LP1)=ONE+((P(I,LP1)-PRESS(I,LM1))*C(I,LLP1)- 1 (P(I,LP1)-PRESS(I,L))*C(I,L))/(PRESS(I,L)-PRESS(I,LM1)) 641 CONTINUE DO 643 K=2,L DO 643 I=MYIS,MYIE CO21(I,K,K)=ONE+HAF*(C(I,LM1+K)+C(I,K-1)) 643 CONTINUE C C COMPUTE NEARBY-LAYER TRANSMISSIVITIES FOR THE O3 BAND AND FOR THE C ONE-BAND CONTINUUM BAND (TO3 AND EMISS2). THE SF2 FUNCTION IS C USED. THE METHOD IS THE SAME AS DESCRIBED FOR CO2 IN REF (4). DO 651 K=1,LM1 DO 651 I=MYIS,MYIE CSUB(I,K+1)=CNTVAL(I,K+1)*DELPR1(I,K+1) CSUB(I,LP1+K-1)=CNTVAL(I,K)*DELPR2(I,K+1) 651 CONTINUE C---THE SF2 FUNCTION IN PREV. VERSIONS IS NOW EXPLICITLY EVALUATED DO 655 K=1,LLM2 DO 655 I=MYIS,MYIE CSUB2(I,K+1)=SKO3R*CSUB(I,K+1) C(I,K+1)=CSUB(I,K+1)*(HMP5+CSUB(I,K+1)* 1 (HP166666-CSUB(I,K+1)*H41666M2)) C2(I,K+1)=CSUB2(I,K+1)*(HMP5+CSUB2(I,K+1)* 1 (HP166666-CSUB2(I,K+1)*H41666M2)) 655 CONTINUE DO 661 I=MYIS,MYIE CONTDG(I,LP1)=1.+C(I,LLM1) TO3DG(I,LP1)=1.+C2(I,LLM1) 661 CONTINUE DO 663 K=2,L DO 663 I=MYIS,MYIE CONTDG(I,K)=ONE+HAF*(C(I,K)+C(I,LM1+K)) TO3DG(I,K)=ONE+HAF*(C2(I,K)+C2(I,LM1+K)) 663 CONTINUE C---NOW OBTAIN FLUXES C C FOR THE DIAGONAL TERMS... DO 871 K=2,LP1 DO 871 I=MYIS,MYIE FLX(I,K)=FLX(I,K)+(DTC(I,K)*EMISDG(I,K) 1 +SS2(I,K)*CONTDG(I,K) 2 +OSS(I,K)*TO3DG(I,K) 3 +CSS(I,K)*CO21(I,K,K))*CLDFAC(I,K,K) 871 CONTINUE C FOR THE TWO OFF-DIAGONAL TERMS... DO 873 I=MYIS,MYIE FLX(I,L)=FLX(I,L)+(CSS(I,LP1)*CO21(I,LP1,L) 1 +DTC(I,LP1)*EMSPEC(I,2) 2 +OSS(I,LP1)*TO31D(I,L) 3 +SS2(I,LP1)*CONT1D(I,L))*CLDFAC(I,LP1,L) FLX(I,LP1)=FLX(I,LP1)+(CSS(I,L)*CO21(I,L,LP1) 1 +OSS(I,L)*TO31D(I,L) 2 +SS2(I,L)*CONT1D(I,L) 3 +DTC(I,L)*EMSPEC(I,1))*CLDFAC(I,L,LP1) 873 CONTINUE C C FINAL SECTION OBTAINS EMISSIVITY HEATING RATES, C TOTAL HEATING RATES AND THE FLUX AT THE GROUND C C .....CALCULATE THE EMISSIVITY HEATING RATES DO 1101 K=1,L DO 1101 I=MYIS,MYIE HEATEM(I,K)=RADCON*(FLX(I,K+1)-FLX(I,K))*DELP(I,K) 1101 CONTINUE C .....CALCULATE THE TOTAL HEATING RATES DO 1103 K=1,L DO 1103 I=MYIS,MYIE c if(mype.eq.13.and.i.eq.40) then c print*,'k=',k c print*,'cts(i,k)=',cts(i,k) c print*,'ctso(i,k)=',ctso3(i,k) c print*,'excts(i,k)=',excts(i,k) c endif HEATRA(I,K)=HEATEM(I,K)-CTS(I,K)-CTSO3(I,K)+EXCTS(I,K) 1103 CONTINUE C .....CALCULATE THE FLUX AT EACH FLUX LEVEL USING THE FLUX AT THE C TOP (FLX1E1+GXCTS) AND THE INTEGRAL OF THE HEATING RATES (VSUM1) DO 1111 K=1,L DO 1111 I=MYIS,MYIE VSUM1(I,K)=HEATRA(I,K)*DELP2(I,K)*RADCON1 1111 CONTINUE DO 1115 I=MYIS,MYIE TOPFLX(I)=FLX1E1(I)+GXCTS(I) c if(mype.eq.13.and.i.eq.40) then c print*,'flx1e1(i),gxcts(i)=',flx1e1(i),gxcts(i) c print*,'topflx(i)=',topflx(i) c endif FLXNET(I,1)=TOPFLX(I) 1115 CONTINUE C---ONLY THE SURFACE VALUE OF FLUX (GRNFLX) IS NEEDED UNLESS C THE THICK CLOUD SECTION IS INVOKED. DO 1123 K=2,LP1 DO 1123 I=MYIS,MYIE c if(mype.eq.13.and.i.eq.40) then c print*,'k,k-1,flxnet(i,k-1),vsum1(i,k-1)=', c * k,k-1,flxnet(i,k-1),vsum1(i,k-1) c endif FLXNET(I,K)=FLXNET(I,K-1)+VSUM1(I,K-1) 1123 CONTINUE DO 1125 I=MYIS,MYIE GRNFLX(I)=FLXNET(I,LP1) 1125 CONTINUE C C THIS IS THE THICK CLOUD SECTION.OPTIONALLY,IF THICK CLOUD C FLUXES ARE TO BE "CONVECTIVELY ADJUSTED",IE,DF/DP IS CONSTANT, C FOR CLOUDY PART OF GRID POINT, THE FOLLOWING CODE IS EXECUTED. C***FIRST,COUNT THE NUMBER OF CLOUDS ALONG THE LAT. ROW. SKIP THE C ENTIRE THICK CLOUD COMPUTATION OF THERE ARE NO CLOUDS. ICNT=0 DO 1301 I=MYIS,MYIE ICNT=ICNT+NCLDS(I) 1301 CONTINUE IF (ICNT.EQ.0) GO TO 6999 C---FIND THE MAXIMUM NUMBER OF CLOUDS IN THE LATITUDE ROW KCLDS=NCLDS(1) DO 2106 I=MYIS,MYIE KCLDS=MAX(NCLDS(I),KCLDS) 2106 CONTINUE C C C***OBTAIN THE PRESSURES AND FLUXES OF THE TOP AND BOTTOM OF C THE NC'TH CLOUD (IT IS ASSUMED THAT ALL KTOP AND KBTM'S HAVE C BEEN DEFINED!). DO 1361 KK=1,KCLDS KMIN=LP1 KMAX=0 DO 1362 I=MYIS,MYIE J1=KTOP(I,KK+1) C IF (J1.EQ.1) GO TO 1362 J3=KBTM(I,KK+1) IF (J3.GT.J1) THEN PTOP(I)=P(I,J1) PBOT(I)=P(I,J3+1) FTOP(I)=FLXNET(I,J1) FBOT(I)=FLXNET(I,J3+1) C***OBTAIN THE "FLUX DERIVATIVE" DF/DP (DELPTC) DELPTC(I)=(FTOP(I)-FBOT(I))/(PTOP(I)-PBOT(I)) KMIN=MIN(KMIN,J1) KMAX=MAX(KMAX,J3) ENDIF 1362 CONTINUE KMIN=KMIN+1 C***CALCULATE THE TOT. FLUX CHG. FROM THE TOP OF THE CLOUD, FOR C ALL LEVELS. DO 1365 K=KMIN,KMAX DO 1363 I=MYIS,MYIE C IF (KTOP(I,KK+1).EQ.1) GO TO 1363 IF(KTOP(I,KK+1).LT.K .AND. K.LE.KBTM(I,KK+1)) THEN Z1(I,K)=(P(I,K)-PTOP(I))*DELPTC(I)+FTOP(I) CORIGINAL FLXNET(I,K)=FLXNET(I,K)*(ONE-CAMT(I,KK+1)) + CORIGINAL1 Z1(I,K)*CAMT(I,KK+1) c if(mype.eq.13.and.i.eq.40) then c print*,'k,z1(i,k)=',k,z1(i,k) c endif FLXNET(I,K)=Z1(I,K) ENDIF 1363 CONTINUE 1365 CONTINUE 1361 CONTINUE C***USING THIS FLUX CHG. IN THE CLOUDY PART OF THE GRID BOX, OBTAIN C THE NEW FLUXES, WEIGHTING THE CLEAR AND CLOUDY FLUXES:AGAIN, ONLY C THE FLUXES IN THICK-CLOUD LEVELS WILL EVENTUALLY BE USED. C DO 6051 K=1,LP1 C DO 6051 I=MYIS,MYIE C FLXNET(I,K)=FLXNET(I,K)*(ONE-CAMT(I,NC)) + C 1 Z1(I,K)*CAMT(I,NC) C051 CONTINUE C***MERGE FLXTHK INTO FLXNET FOR APPROPRIATE LEVELS. C DO 1401 K=1,LP1 C DO 1401 I=MYIS,MYIE C IF (K.GT.ITOP(I) .AND. K.LE.IBOT(I) C 1 .AND. (NC-1).LE.NCLDS(I)) THEN if(mype.eq.13.and.i.eq.40) then print*,'k,flxthk(i,k)=',k,flxthk(i,k) endif C FLXNET(I,K)=FLXTHK(I,K) C ENDIF C401 CONTINUE C C******END OF CLOUD LOOP***** 6001 CONTINUE 6999 CONTINUE C***THE FINAL STEP IS TO RECOMPUTE THE HEATING RATES BASED ON THE C REVISED FLUXES: DO 6101 K=1,L DO 6101 I=MYIS,MYIE c if(mype.eq.13.and.i.eq.40) then c print*,'i,k=',i,k c print*,'radcon=',radcon c print*,'flxnet(i,k+1)=',flxnet(i,k+1) c print*,'flxnet(i,k)=',flxnet(i,k) c print*,'delp(i,k)=',delp(i,k) c endif HEATRA(I,K)=RADCON*(FLXNET(I,K+1)-FLXNET(I,K))*DELP(I,K) 6101 CONTINUE C THE THICK CLOUD SECTION ENDS HERE. RETURN END