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