FUNCTION ANTEMP(L,Z) DIMENSION ZB(10,7),C(11,7),DELTA(10,7),TSTAR(7) C ************** TROPICAL SOUNDING ************************** DATA (ZB(N,1),N=1,10)/ 1 2.0, 3.0, 16.5, 21.5, 45.0, 2 51.0, 70.0, 100., 200., 300./ DATA (C(N,1),N=1,11)/ 1 -6.0, -4.0, -6.7, 4.0, 2.2, 2 1.0, -2.8, -.27, 0.0, 0.0, 0.0/ DATA (DELTA(N,1),N=1,10)/ 1 .5, .5, .3, .5, 1.0, 2 1.0, 1.0, 1.0, 1.0, 1.0/ C ************** SUB-TROPICAL SUMMER ************************ DATA (ZB(N,2),N=1,10)/ 1 1.5, 6.5, 13.0, 18.0, 26.0, 2 36.0, 48.0, 50.0, 70.0, 100./ DATA (C(N,2),N=1,11)/ 1 -4.0, -6.0, -6.5, 0.0, 1.2, 2 2.2, 2.5, 0.0, -3.0, -0.25, 0.0/ DATA (DELTA(N,2),N=1,10)/ 1 .5, 1.0, .5, .5, 1.0, 2 1.0, 2.5, .5, 1.0, 1.0/ C ************** SUB-TROPICAL WINTER ************************ DATA (ZB(N,3),N=1,10)/ 1 3.0, 10.0, 19.0, 25.0, 32.0, 2 44.5, 50.0, 71.0, 98.0, 200.0/ DATA (C(N,3),N=1,11)/ 1 -3.5, -6.0, -0.5, 0.0, 0.4, 2 3.2, 1.6, -1.8, -0.7, 0.0, 0.0/ DATA (DELTA(N,3),N=1,10)/ 1 .5, .5, 1.0, 1.0, 1.0, 2 1.0, 1.0, 1.0, 1.0, 1.0/ C ************* SUB-ARCTIC SUMMER ************************* DATA (ZB(N,4),N=1,10)/ 1 4.7, 10.0, 23.0, 31.8, 44.0, 2 50.2, 69.2, 100.0, 102.0, 103.0/ DATA (C(N,4),N=1,11)/ 1 -5.3, -7.0, 0.0, 1.4, 3.0, 2 0.7, -3.3, -0.2, 0.0, 0.0, 0.0/ DATA (DELTA(N,4),N=1,10)/ 1 .5, .3, 1.0, 1.0, 2.0, 2 1.0, 1.5, 1.0, 1.0, 1.0/ C ************ SUB-ARCTIC WINTER ***************************** DATA (ZB(N,5),N=1,10)/ 1 1.0, 3.2, 8.5, 15.5, 25.0, 2 30.0, 35.0, 50.0, 70.0, 100.0/ DATA (C(N,5),N=1,11)/ 1 3.0, -3.2, -6.8, 0.0, -0.6, 2 1.0, 1.2, 2.5, -0.7, -1.2, 0.0/ DATA (DELTA(N,5),N=1,10)/ 1 .4, 1.5, .3 , .5, 1.0, 2 1.0, 1.0, 1.0, 1.0, 1.0/ C ************ US STANDARD 1976 ****************************** DATA (ZB(N,6),N=1,10)/ 1 11.0, 20.0, 32.0, 47.0, 51.0, 2 71.0, 84.8520, 90.0, 91.0, 92.0/ DATA (C(N,6),N=1,11)/ 1 -6.5, 0.0, 1.0, 2.80, 0.0, 2 -2.80, -2.00, 0.0, 0.0, 0.0, 0.0/ DATA (DELTA(N,6),N=1,10)/ 1 0.3, 1.0, 1.0, 1.0, 1.0, 2 1.0, 1.0, 1.0, 1.0, 1.0/ C C ************ ENLARGED US STANDARD 1976 ********************** DATA (ZB(N,7),N=1,10)/ 1 11.0, 20.0, 32.0, 47.0, 51.0, 2 71.0, 84.8520, 90.0, 91.0, 92.0/ DATA (C(N,7),N=1,11)/ 1 -6.5, 0.0, 1.0, 2.80, 0.0, 2 -2.80, -2.00, 0.0, 0.0, 0.0, 0.0/ DATA (DELTA(N,7),N=1,10)/ 1 0.3, 1.0, 1.0, 1.0, 1.0, 2 1.0, 1.0, 1.0, 1.0, 1.0/ C DATA TSTAR/ 1 300.0, 294.0, 272.2, 287.0, 257.1, 2*288.15/ C NLAST=10 TEMP=TSTAR(L)+C(1,L)*Z DO 20 N=1,NLAST EXPO=(Z-ZB(N,L))/DELTA(N,L) EXPP=ZB(N,L)/DELTA(N,L) FAC=EXP(EXPP)+EXP(-EXPP) IF(ABS(EXPO)-100.0) 23,23,24 23 X=EXP(EXPO) Y=X+1.0/X ZLOG=ALOG(Y) GO TO 25 24 ZLOG=ABS(EXPO) 25 IF(EXPP-100.0) 27,27,28 27 FACLOG=ALOG(FAC) GO TO 29 28 FACLOG=EXPP C TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)* C 1 ALOG((EXP(EXPO)+EXP(-EXPO))/FAC)) 29 TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)* 1 (ZLOG-FACLOG)) 20 CONTINUE ANTEMP=TEMP RETURN END SUBROUTINE COEINT(RAT,IR) C ********************************************************************** C C C THE TRANSMISSION FUNCTION BETWEEN P1 AND P2 IS ASSUMED TO C THE FUNCTIONAL FORM C TAU(P1,P2)= 1.0-SQRT(C*LOG(1.0+X*PATH)), C WHERE C PATH(P1,P2)=((P1-P2)**2)*(P1+P2+CORE)/ C (ETA*(P1+P2+CORE)+(P1-P2)) C C C THE PARAMETERS C AND X ARE FUNCTIONS OF P2, AND ARE TO BE DETER C WHILE CORE IS A PRESPECIFIED NUMBER.ETA IS A FUNCTION OF THE TH C PRODUCT (CX);IT IS OBTAITED ITERATIVELY. THE DERIVATION OF ALL C VALUES WILL BE EXPLAINED IN A FORTHCOMING PAPER. C SUBROUTINE COEINT DETERMINES C(I) AND X(I) BY USING THE ACT C VALUES OF TAU(P(I-2),P(I)) AND TAU(P(I-1),P(I)) AND THE PREVIOU C ITERATION VALUE OF ETA. C DEFINE: C PATHA=PATH(P(I),P(I-2),CORE,ETA) C PATHB=PATH(P(I),P(I-1),CORE,ETA); C THEN C R=(1-TAU(P(I),P(I-2)))/(1-TAU(P(I),P(I-1))) C = SQRT(LOG(1+X*PATHA)/LOG(1+X*PATHB)), C SO THAT C R**2= LOG(1+X*PATHA)/LOG(1+X*PATHB). C THIS EQUATION CAN BE SOLVED BY NEWTON S METHOD FOR X AND THEN T C RESULT USED TO FIND C. THIS IS REPEATED FOR EACH VALUE OF I GRE C THAN 2 TO GIVE THE ARRAYS X(I) AND C(I). C NEWTON S METHOD FOR SOLVING THE EQUATION C F(X)=0 C MAKES USE OF THE LOOP XNEW= XOLD-F(XOLD)/FPRIME(XOLD). C THIS IS ITERATED 20 TIMES, WHICH IS PROBABLY EXCESSIVE. C THE FIRST GUESS FOR ETA IS 3.2E-4*EXP(-P(I)/1000),WHICH HAS C BEEN FOUND TO BE FAIRLY REALISTIC BY EXPERIMENT; WE ITERATE 5 T C (AGAIN,PROBABLY EXCESSIVELY) TO OBTAIN THE VALUES FOR C,X,ETA T C USED FOR INTERPOLATION. C THERE ARE SEVERAL POSSIBLE PITFALLS: C 1) IN THE COURSE OF ITERATION, X MAY REACH A VALUE WHICH C 1+X*PATHA NEGATIVE; IN THIS CASE THE ITERATION IS STOP C AND AN ERROR MESSAGE IS PRINTED OUT. C 2) EVEN IF (1) DOES NOT OCCUR, IT IS STILL POSSIBLE THAT C BE NEGATIVE AND LARGE ENOUGH TO MAKE 1+X*PATH(P(I),0,C C NEGATIVE. THIS IS CHECKED FOR IN A FINAL LOOP, AND IF C A WARNING IS PRINTED OUT. C C ********************************************************************* C.... C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/PRESS/PA(109) REAL RAT,SINV REAL PA,CORE,TRANSA,PATH,UEXP,SEXP,ETA,SEXPV REAL PA2 COMMON/TRAN/ TRANSA(109,109) COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP DIMENSION PATH0(109),ETAP(109),XAP(109),CAP(109) DIMENSION SINV(4) DATA SINV/2.74992,2.12731,4.38111,0.0832926/ CNOV89 DIMENSION SINV(3) CNOV89 DATA SINV/2.74992,2.12731,4.38111/ CO222 OLD CODE USED 2.7528 RATHER THAN 2.74992 ---K.A.C. OCTOBER 1988 CO222 WHEN 2.7528 WAS USED,WE EXACTLY REPRODUCED THE MRF CO2 ARRAYS CORE=5.000 UEXP=0.90 P0=0.7 DO 902 I=1,109 PA2=PA(I)*PA(I) SEXPV(I)=.505+2.0E-5*PA(I)+.035*(PA2-.25)/(PA2+.25) 902 CONTINUE DO 900 I=1,109 ETA(I)=3.2E-4*EXP(-PA(I)/500.) ETAP(I)=ETA(I) 900 CONTINUE DO 1200 NP=1,10 DO 1000 I=3,109 SEXP=SEXPV(I) R=(1.0D0-TRANSA(I,I-2))/(1.0D0-TRANSA(I,I-1)) REXP=R**(UEXP/SEXP) PATHA=(PATH(PA(I),PA(I-2),CORE,ETA(I)))**UEXP PATHB=(PATH(PA(I),PA(I-1),CORE,ETA(I)))**UEXP XX=2.0D0*(PATHB*REXP-PATHA)/(PATHB*PATHB*REXP-PATHA*PATHA) DO 1010 LL=1,20 F1=DLOG(1.0D0+XX*PATHA) F2=DLOG(1.0D0+XX*PATHB) F=F1/F2-REXP FPRIME=(F2*PATHA/(1.0D0+XX*PATHA)-F1*PATHB/(1.0D0+XX*PATHB))/ 1 (F2*F2) XX=XX-F/FPRIME CHECK=1.0D0+XX*PATHA IF (CHECK) 1020,1020,1025 1020 CONTINUE C+++ WRITE (6,360) I,LL,CHECK 360 FORMAT (' ERROR,I=',I3,'LL=',I3,'CHECK=',F20.10) STOP 1025 CONTINUE 1010 CONTINUE CA(I)=(1.0D0-TRANSA(I,I-2))**(UEXP/SEXP)/ 1 (DLOG(1.0D0+XX*PATHA)+1.0D-20) XA(I)=XX 1000 CONTINUE XA(2)=XA(3) XA(1)=XA(3) CA(2)=CA(3) CA(1)=CA(3) DO 1100 I=3,109 PATH0(I)=(PATH(PA(I),0.,CORE,ETA(I)))**UEXP PATH0(I)=1.0D0+XA(I)*PATH0(I) C+++ IF (PATH0(I).LT.0.) WRITE (6,361) I,PATH0(I),XA(I) 1100 CONTINUE DO 1035 I=1,109 SEXP=SEXPV(I) ETAP(I)=ETA(I) ETA(I)=(SINV(IR)/RAT)**(1./SEXP)* 1 (CA(I)*XA(I))**(1./UEXP) 1035 CONTINUE C C THE ETA FORMULATION IS DETAILED IN SCHWARZKOPF AND FELS(1985). C THE QUANTITY SINV=(G*DELTANU)/(RCO2*D*S) C IN CGS UNITS,WITH D,THE DIFFUSICITY FACTOR=2, AND C S,THE SUM OF CO2 LINE STRENGTHS OVER THE 15UM CO2 BAND C ALSO,THE DENOMINATOR IS MULTIPLIED BY C 1000 TO PERMIT USE OF MB UNITS FOR PRESSURE. C S IS ACTUALLY WEIGHTED BY B(250) AT 10 CM-1 WIDE INTERVALS,IN C ORDER TO BE CONSISTENT WITH THE METHODS USED TO OBTAIN THE LBL C 1-BAND CONSOLIDATED TRANCMISSION FUNCTIONS. C FOR THE 490-850 INTERVAL (DELTANU=360,IR=1) SINV=2.74992. C (SLIGHTLY DIFFERENT FROM 2.7528 USED IN EARLIER VERSIONS) C FOR THE 490-670 INTERVAL (IR=2) SINV=2.12731 C FOR THE 670-850 INTERVAL (IR=3) SINV=4.38111 C FOR THE 2270-2380 INTERVAL (IR=4) SINV=0.0832926 C SINV HAS BEEN OBTAINED USING THE 1982 AFGL CATALOG FOR CO2 C RAT IS THE ACTUAL CO2 MIXING RATIO IN UNITS OF 330 PPMV, C LETTING USE OF THIS FORMULATION FOR ANY CO2 CONCENTRATION. C C WRITE (6,366) (NP,I,CA(I),XA(I),ETA(I),SEXPV(I),I=1,109) C366 FORMAT (2I4,4E20.12) 1200 CONTINUE 361 FORMAT (' **WARNING:** 1+XA*PATH(PA(I),0) IS NEGATIVE,I= ',I3,/ 1 20X,'PATH0(I)=',F16.6,' XA(I)=',F16.6) RETURN END CCCC PROGRAM CO2INS SUBROUTINE CO2INS(T22,T23,T66,IQ) C ********************************************************* C SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 ****** C ..... K.CAMPANA MARCH 1988,OCTOBER 1988... C ..... K.CAMPANA DECEMBER 1988-CLEANED UP FOR LAUNCHER C ..... K.CAMPANA NOVEMBER 1989-ALTERED FOR NEW RADIATION C ********************************************************* INCLUDE "parmco2" PARAMETER (L=LERIC,LP1=L+1) DIMENSION T22(LP1,LP1,3),T23(LP1,LP1,3),T66(LP1,LP1,6) DIMENSION DCDT8(LP1,LP1),DCDT10(LP1,LP1),CO2PO(LP1,LP1), * CO2800(LP1,LP1),CO2PO1(LP1,LP1),CO2801(LP1,LP1),CO2PO2(LP1,LP1), * CO2802(LP1,LP1),N(LP1),D2CT8(LP1,LP1),D2CT10(LP1,LP1) CCC ITIN=22 CCC ITIN1=23 CO222 LATEST CODE HAD IQ=1 CCC IQ=4 1011 FORMAT (4F20.14) CCC READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2802(I,J),I=1,LP1),J=1,LP1) DO 300 J=1,LP1 DO 300 I=1,LP1 CO2PO(I,J) = T22(I,J,1) CNOV89 IF (IQ.EQ.5) GO TO 300 CNOV89 CO2PO1(I,J) = T22(I,J,2) CO2PO2(I,J) = T22(I,J,3) 300 CONTINUE DO 301 J=1,LP1 DO 301 I=1,LP1 CO2800(I,J) = T23(I,J,1) CNOV89 IF (IQ.EQ.5) GO TO 301 CNOV89 CO2801(I,J) = T23(I,J,2) CO2802(I,J) = T23(I,J,3) 301 CONTINUE C***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS C ARE: C IQ=1 560-800 (CONSOL.=490-850) C IQ=2 560-670 (CONSOL.=490-670) C IQ=3 670-800 (CONSOL.=670-850) C IQ=4 560-760 (ORIGINAL CODE) (CONSOL.=490-850) CNOV89 C IQ=5 2270-2380 (CONSOL.=2270-2380) CNOV89 C THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS C USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT C WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S. CNOV89 C NOTE: ALTHOUGH THE BAND TRANSMISSION FUNCTIONS ARE C COMPUTED FOR ALL RADIATIVE BANDS, AS OF 9/28/88, THEY C ARE WRITTEN OUT IN FULL ONLY FOR THE FULL 15 UM BAND CASES C (IQ=1,4). IN OTHER CASES, THE TRANSMISSIVITIES (1,K) ARE C WRITTEN OUT, AS THESE ARE THE ONLY ONES NEEDED FOR CTS C CALCULATIONS. ALSO, FOR THE 4.3 UM BAND (IQ=5) THE TEMP. C DERIVATIVE TERMS ARE NOT WRITTEN OUT, AS THEY ARE UNUSED. CNOV89 IF (IQ.EQ.1) THEN C1=1.5 C2=0.5 ENDIF IF (IQ.EQ.2) THEN C1=18./11. C2=7./11. ENDIF IF (IQ.EQ.3) THEN C1=18./13. C2=5./13. ENDIF IF (IQ.EQ.4) THEN C1=1.8 C2=0.8 ENDIF CNOV89 IF (IQ.EQ.5) THEN C1=1.0 C2=0.0 ENDIF CNOV89 DO 1021 I=1,LP1 DO 1021 J=1,LP1 CO2PO(J,I)=C1*CO2PO(J,I)-C2 CO2800(J,I)=C1*CO2800(J,I)-C2 CNOV89 IF (IQ.EQ.5) GO TO 1021 CNOV89 CO2PO1(J,I)=C1*CO2PO1(J,I)-C2 CO2801(J,I)=C1*CO2801(J,I)-C2 CO2PO2(J,I)=C1*CO2PO2(J,I)-C2 CO2802(J,I)=C1*CO2802(J,I)-C2 1021 CONTINUE CNOV89 IF (IQ.GE.1.AND.IQ.LE.4) THEN CNOV89 DO 1 J=1,LP1 DO 1 I=1,LP1 DCDT8(I,J)=.02*(CO2801(I,J)-CO2802(I,J))*100. DCDT10(I,J)=.02*(CO2PO1(I,J)-CO2PO2(I,J))*100. D2CT8(I,J)=.0016*(CO2801(I,J)+CO2802(I,J)-2.*CO2800(I,J))*1000. D2CT10(I,J)=.0016*(CO2PO1(I,J)+CO2PO2(I,J)-2.*CO2PO(I,J))*1000. 1 CONTINUE CNOV89 ENDIF CNOV89 CO222 ********************************************************* CCC REWIND 66 C SAVE CDT51,CO251,C2D51,CDT58,CO258,C2D58..ON TEMPO FILE CCC WRITE (66) DCDT10 CCC WRITE (66) CO2PO CCC WRITE (66) D2CT10 CCC WRITE (66) DCDT8 CCC WRITE (66) CO2800 CCC WRITE (66) D2CT8 CCC REWIND 66 CNOV89 IF (IQ.EQ.1.OR.IQ.EQ.4) THEN CNOV89 DO 400 J=1,LP1 DO 400 I=1,LP1 T66(I,J,1) = DCDT10(I,J) T66(I,J,2) = CO2PO(I,J) T66(I,J,3) = D2CT10(I,J) T66(I,J,4) = DCDT8(I,J) T66(I,J,5) = CO2800(I,J) T66(I,J,6) = D2CT8(I,J) 400 CONTINUE CNOV89 ELSE DO 409 I=1,LP1 T66(I,1,2) = CO2PO(1,I) T66(I,1,5) = CO2800(1,I) IF (IQ.EQ.5) GO TO 409 T66(I,1,1) = DCDT10(1,I) T66(I,1,3) = D2CT10(1,I) T66(I,1,4) = DCDT8(1,I) T66(I,1,6) = D2CT8(1,I) 409 CONTINUE ENDIF CNOV89 CO222 ********************************************************* RETURN END CO222 PROGRAM CO2INT(INPUT,TAPE5=INPUT) CNOV89 SUBROUTINE CO2INT(ITAPE,T15A,T15B,T22,RATIO,IR,NMETHD) CNOV89 C ********************************************************* C CHANGES TO DATA READ AND FORMAT SEE CO222 *** C ..... K.CAMPANA MARCH 1988,OCTOBER 1988 C CHANGES TO PASS ITAPE,AND IF IR=4,READ 1 CO2 REC..KAC NOV89 C ********************************************************* C CO2INT INTERPOLATES CARBON DIOXIDE TRANSMISSION FUNCTIONS C FROM THE 109 LEVEL GRID,FOR WHICH THE TRANSMISSION FUNCTIONS C HAVE BEEN PRE-CALCULATED, TO THE GRID STRUCTURE SPECIFIED BY THE C USER. C C METHOD: C C CO2INT IS EMPLOYABLE FOR TWO PURPOSES: 1) TO OBTAIN TRANSMIS- C SIVITIES BETWEEN ANY 2 OF AN ARRAY OF USER-DEFINED PRESSURES; AND C 2) TO OBTAIN LAYER-MEAN TRANSMISSIVITIES BETWEEN ANY 2 OF AN ARRAY C OF USER-DEFINED PRESSURE LAYERS.TO CLARIFY THESE TWO PURPOSES,SEE C THE DIAGRAM AND DISCUSSION BELOW. C CO2INT MAY BE USED TO EXECUTE ONLY ONE PURPOSE AT ONE TIME. C C LET P BE AN ARRAY OF USER-DEFINED PRESSURES C AND PD BE USER-DEFINED PRESSURE LAYERS. C C - - - - - - - - - PD(I-1) --- C ^ C ----------------- P(I) ^ PRESSURE LAYER I (PLM(I)) C ^ C - - - - - - - - - PD(I) --- C ^ C ----------------- P(I+1) ^ PRESSURE LAYER I+1 (PLM(I+1)) C ^ C - - - - - - - - - PD(I+1)--- C ... (THE NOTATION USED IS C ... CONSISTENT WITH THE CODE) C ... C - - - - - - - - - PD(J-1) C C ----------------- P(J) C C - - - - - - - - - PD(J) C C PURPOSE 1: THE TRANSMISSIVITY BETWEEN SPECIFIC PRESSURES C P(I) AND P(J) ,TAU(P(I),P(J)) IS COMPUTED BY THIS PROGRAM. C IN THIS MODE,THERE IS NO REFERENCE TO LAYER PRESSURES PD C (PD,PLM ARE NOT INPUTTED). C C PURPOSE 2: THE LAYER-MEAN TRANSMISSIVITY BETWEEN A LAYER- C MEAN PRESSURE PLM(J) AND PRESSURE LAYER I IS GIVEN BY C TAULM(PLM(I),PLM(J)). IT IS COMPUTED BY THE INTEGRAL C C PD(I) C ---- C 1 ^ C ------------- * ^ TAU ( P',PLM(J) ) DP' C PD(I)-PD(I-1) ^ C ---- C PD(I-1) C C THE LAYER-MEAN PRESSURE PLM(I) IS SPECIFIED BY THE USER. C FOR MANY PURPOSES,PLM WILL BE CHOSEN TO BE THE AVERAGE C PRESSURE IN THE LAYER-IE,PLM(I)=0.5*(PD(I-1)+PD(I)). C FOR LAYER-MEAN TRANSMISSIVITIES,THE USER THUS INPUTS C A PRESSURE ARRAY (PD) DEFINING THE PRESSURE LAYERS AND AN C ARRAY (PLM) DEFINING THE LAYER-MEAN PRESSURES.THE CALCULATION C DOES NOT DEPEND ON THE P ARRAY USED FOR PURPOSE 1 (P IS NOT C INPUTTED). C C THE FOLLOWING PARAGRAPHS DEPICT THE UTILIZATION OF THIS C CODE WHEN USED TO COMPUTE TRANSMISSIVITIES BETWEEN SPECIFIC C PRESSURES. LATER PARAGRAPHS DESCRIBE ADDITIONAL FEATURES NEEDED C FOR LAYER-MEAN TRANSMISSIVITIES. C C FOR A GIVEN CO2 MIXING RATIO AND STANDARD TEMPERATURE C PROFILE,A TABLE OF TRANSMISSION FUNCTIONS FOR A FIXED GRID C OF ATMOSPHERIC PRESSURES HAS BEEN PRE-CALCULATED. C THE STANDARD TEMPERATURE PROFILE IS COMPUTED FROM THE US C STANDARD ATMOSPHERE (1977) TABLE.ADDITIONALLY, THE C SAME TRANSMISSION FUNCTIONS HAVE BEEN PRE-CALCULATED FOR A C TEMPERATURE PROFILE INCREASED AND DECREASED (AT ALL LEVELS) C BY 25 DEGREES. C THIS PROGRAM READS IN THE PRESPECIFIED TRANSMISSION FUNCTIONS C AND A USER-SUPPLIED PRESSURE GRID (P(I)) AND CALCULATES TRANS- C MISSION FUNCTIONS ,TAU(P(I),P(J)), FOR ALL P(I) S AND P(J) S. C A LOGARITHMIC INTERPOLATION SCHEME IS USED. C THIS METHOD IS REPEATED FOR THE THREE TEMPERATURE PROFILES C GIVEN ABOVE .THEREFORE OUTPUTS FROM THE PROGRAM ARE THREE TABLES C OF TRANSMISSION FUNCTIONS FOR THE USER-SUPPLIED PRESSURE GRID. C THE EXISTENCE OF THE THREE TABLES PERMITS SUBSEQUENT INTERPO- C LATION TO A USER-SUPPLIED TEMPERATURE PROFILE USING THE METHOD C DESCRIBED IN THE REFERENCE.SEE LIMITATIONS SECTION IF THE C USER DESIRES TO OBTAIN ONLY 1 TABLE OF TRANSMISSIVITIES. C C MODIFICATIONS FOR LAYER-MEAN TRANSMISSIVITIES: C THE PRESSURES INPUTTED ARE THE LAYER-MEAN PRESSURES,PD, C AND THE LAYER-MEAN PRESSURES ,PLM. A SERIES OF TRANSMISSIVITIES C (TAU(P'',PLM(J)) ARE COMPUTED AND THE INTEGRAL GIVEN IN THE C DISCUSSION OF PURPOSE 2 IS COMPUTED.FOR PLM(I) NOT EQUAL TO C PLM(J) SIMPSON S RULE IS USED WITH 5 POINTS. IF PLM(I)=PLM(J) C (THE -NEARBY LAYER- CASE) A 49-POINT QUADRATURE IS USED FOR C GREATER ACCURACY.THE OUTPUT IS IN TAULM(PLM(I),PLM(J)). C NOTE: C TAULM IS NOT A SYMMETRICAL MATRIX. FOR THE ARRAY ELEMENT C TAULM(PLM(I),PLM(J)),THE INNER(FIRST,MOST RAPIDLY VARYING) C DIMENSION IS THE VARYING LAYER-MEAN PRESSURE,PLM(I);THE OUTER C (SECOND) DIMENSION IS THE FIXED LAYER-MEAN PRESSURE PLM(J). C THUS THE ELEMENT TAULM(2,3) IS THE TRANSMISSION FUNCTION BETWEEN C THE FIXED PRESSURE PLM(3) AND THE PRESSURE LAYER HAVING AN AVERAG C PRESSURE OF PLM(2). C ALSO NOTE THAT NO QUADRATURE IS PERFORMED OVER THE LAYER C BETWEEN THE SMALLEST NONZERO PRESSURE AND ZERO PRESSURE; C TAULM IS TAULM(0,PLM(J)) IN THIS CASE,AND TAULM(0,0)=1. C C C REFERENCE: C S.B.FELS AND M.D.SCHWARZKOPF,-AN EFFICIENT ACCURATE C ALGORITHM FOR CALCULATING CO2 15 UM BAND COOLING RATES-,JOURNAL C OF GEOPHYSICAL RESEARCH,VOL.86,NO. C2, PP.1205-1232,1981. C MODIFICATIONS TO THE ALGORITHM HAVE BEEN MADE BY THE AUTHORS; C CONTACT S.B.F.OR M.D.S. FOR FURTHER DETAILS.A NOTE TO J.G.R. C IS PLANNED TO DOCUMENT THESE CHANGES. C C AUTHOR: M.DANIEL SCHWARZKOPF C C DATE: 14 JULY 1983 C C ADDRESS: C C G.F.D.L. C P.O.BOX 308 C PRINCETON,N.J.08540 C U.S.A. C TELEPHONE: (609) 452-6521 C C INFORMATION ON TAPE: THIS SOURCE IS THE FIRST FILE C ON THIS TAPE.THE SIX FILES THAT FOLLOW ARE CO2 TRANS- C MISSIVITIES FOR THE 500-850 CM-1 INTERVAL FOR CO2 C CONCENTRATIONS OF 330 PPMV (1X) ,660 PPMV (2X), AND C 1320 PPMV (4X). THE FILES ARE ARRANGED AS FOLLOWS: C FILE 2 1X,CONSOLIDATED USING B(250) WEIGHTING FCTN. C FILE 3 1X,CONSOLIDATED WITH NO WEIGHTING FCTN. C FILE 4 2X,CONSOLIDATED USING B(250) WEIGHTING FCTN. C FILE 5 2X,CONSOLIDATED WITH NO WEIGHTING FCTN. C FILE 6 4X,CONSOLIDATED USING B(250) WEIGHTING FCTN. C FILE 7 4X,CONSOLIDATED WITH NO WEIGHTING FCTN. C FILES 2,4,6 ARE RECOMMENDED FOR USE IN OBTAINING C TRANSMISSION FUNCTIONS FOR USE IN HEATING RATE C COMPUTATIONS;THEY CORRESPOND TO THE TRANSMISSIVITIES C DISCUSSED IN THE 1980 PAPER.FILES 3,5,7 ARE PROVIDED C TO FACILITATE COMPARISON WITH OBSERVATION AND WITH OTHER C CALCULATIONS. C C PROGRAM LANGUAGE: FORTRAN 1977,INCLUDING PARAMETER C AND PROGRAM STATEMENTS.THE PROGRAM IS WRITTEN ON A C CYBER 170-730.SEE THE SECTION ON LIMITATIONS FOR C ADAPTATIONS TO OTHER MACHINES. C C INPUT UNITS,FORMATS AND FORMAT STATEMENT NOS: C C UNIT NO VARIABLES FORMAT STATEMENT NO. TYPE C 5 P (PURPOSE 1) (5E16.9) 201 CARDS C 5 PD (PURPOSE 2) (5E16.9) 201 CARDS C 5 PLM(PURPOSE 2) (5E16.9) 201 CARDS C 5 NMETHD (I3) 202 CARDS C 20 TRANSA (4F20.14) 102 TAPE CNOV89 C ITAPE TRANSA (4F20.14) 102 TAPE CNOV89 C C OUTPUT UNITS,FORMATS AND FORMAT STATEMENT NOS: C C UNIT NO VARIABLES FORMAT STATEMENT NO. C 6 TRNFCT (1X,8F15.8) 301 PRINT C 22 TRNFCT (4F20.14) 102 TAPE C C PARAMETER INPUTS: C A) NLEVLS : NLEVLS IS AN (INTEGER) PARAMETER DENOTING C THE NUMBER OF NONZERO PRESSURE LEVELS FOR PURPOSE 1 C OR THE NUMBER OF NONZERO LAYER PRESSURES NEEDED TO C SPECIFY THE PRESSURE LAYERS(PURPOSE 2) IN THE OUTPUT C GRID. FOR EXAMPLE,IN PURPOSE 1,IF P=0,100,1000,NLEVLS=2. C IF,IN PURPOSE 2,PD=0,100,500,1000,THE NUMBER OF NONZERO C PRESSURE LAYERS=2,SO NLEVLS=2 C IN THE CODE AS WRITTEN,NLEVLS=40; THE USER SHOULD C CHANGE THIS VALUE TO A USER-SPECIFIED VALUE. C B) NLP1,NLP2 : INTEGER PARAMETERS DEFINED AS: NLP1=NLEVLS+1; C NLP2=NLEVLS+2. C SEE LIMITATIONS FOR CODE MODIFICATIONS IF PARAMETER C STATEMENTS ARE NOT ALLOWED ON YOUR MACHINE. C C INPUTS: C C A) TRANSA : THE 109X109 GRID OF TRANSMISSION FUNCTIONS C TRANSA IS A DOUBLE PRECISION REAL ARRAY. C C TRANSA IS READ FROM FILE 20. THIS FILE CONTAINS 3 C RECORDS,AS FOLLOWS: C 1) TRANSA, STANDARD TEMPERATURE PROFILE C 3) TRANSA, STANDARD TEMPERATURES + 25 DEG C 5) TRANSA, STANDARD TEMPERATURES - 25 DEG C C B) NMETHD: AN INTEGER WHOSE VALUE IS EITHER 1 (IF CO2INT IS C TO BE USED FOR PURPOSE 1) OR 2 (IF CO2INT IS TO BE USED FOR C PURPOSE 2). C C C) P,PD,PLM : C P IS A REAL ARRAY (LENGTH NLP1) SPECIFYING THE PRESSURE C GRID AT WHICH TRANSMISSION FUNCTIONS ARE TO BE COMPUTED FOR C PURPOSE 1.THE DIMENSION OF P IS IN MILLIBARS.THE C FOLLOWING LIMITATIONS WILL BE EXPLAINED MORE C IN THE SECTION ON LIMITATIONS: P(1) MUST BE ZERO; P(NLP1),THE C LARGEST PRESSURE, MUST NOT EXCEED 1165 MILLIBARS. C PD IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE PRESSURE C LAYERS FOR WHICH LAYER-AVERAGED TRANSMISSION FUNCTIONS ARE C TO BE COMPUTED.THE DIMENSION OF PD IS MILLIBARS.THE LIMITATIONS C FOR PD ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON C LIMITATIONS. C PLM IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE LAYER-MEAN C PRESSURES. THE DIMENSION OF PLM IS MILLIBARS. THE LIMITATIONS C FOR PLM ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON C LIMITATIONS.PD IS READ IN BEFORE PLM. C C NOTE: AGAIN,WE NOTE THAT THE USER WILL INPUT EITHER P (FOR C PURPOSE 1) OR PD AND PLM(FOR PURPOSE 2) BUT NOT BOTH. C C C C C LIMITATIONS: C 1) P(1)=0.,PD(1)=0.,PLM(1)=0. THE TOP PRESSURE LEVEL C MUST BE ZERO,OR THE TOP PRESSURE LAYER MUST BE BOUNDED BY ZERO. C THE TOP LAYER-MEAN PRESSURE (PLM(1)) MUST BE ZERO; NO C QUADRATURE IS DONE ON THE TOP PRESSURE LAYER.EVEN IF ONE IS C NOT INTERESTED IN THE TRANSMISSION FUNCTION BETWEEN 0 AND P(J), C ONE MUST INCLUDE SUCH A LEVEL. C 2) PD(NLP2)=P(NLP1) IS LESS THAN OR EQUAL TO 1165 MB. C EXTRAPOLATION TO HIGHER PRESSURES IS NOT POSSIBLE. C 3) IF PROGRAM IS NOT PERMITTED ON YOUR COMPILER, C SIMPLY DELETE THE LINE. C 4) IF PARAMETER IS NOT PERMITTED,DO THE FOLLOWING: C 1) DELETE ALL PARAMETER STATEMENTS IN CO2INT C 2) AT THE POINT WHERE NMETHOD IS READ IN,ADD: C READ (5,202) NLEVLS C NLP1=NLEVLS+1 C NLP2=NLEVLS+2 C 3) CHANGE DIMENSION AND/OR COMMON STATEMENTS DEFINING C ARRAYS TRNS,DELTA,P,PD,TRNFCT,PS,PDS,PLM IN CO2INT. C THE NUMERICAL VALUE OF (NLEVLS+1) SHOULD BE INSERTED C IN DIMENSION OR COMMON STATEMENTS FOR TRNS,DELTA, C P,TRNFCT,PS,PLM; THE NUMERICAL VALUE OF (NLEVLS+2) C IN DIMENSION OR COMMON STATEMENTS FOR PD,PDS. C 5) PARAMETER (NLEVLS=40) AND THE OTHER PARAMETER C STATEMENTS ARE WRITTEN IN CDC FORTRAN; ON OTHER MACHINES THE C SAME STATEMENT MAY BE WRITTEN DIFFERENTLY,FOR EXAMPLE AS C PARAMETER NLEVLS=40 C 6) -DOUBLE PRECISION- IS USED INSTEAD OF -REAL*8- ,DUE TO C REQUIREMENTS OF CDC FORTAN. C 7) THE STATEMENT -DO 400 KKK=1,3- CONTROLS THE NUMBER OF C TRANSMISSIVITY OUTPUT MATRICES PORDUCED BY THE PROGRAM.TO C PRODUCE 1 OUTPUT MATRIX,DELETE THIS STATEMENT. C C OUTPUT: C A) TRNFCT IS AN (NLP1,NLP1) REAL ARRAY OF THE TRANSMISSION C FUNCTIONS APPROPRIATE TO YOUR ARRAY. IT IS TO BE SAVED ON FILE 22. C THE PROCEDURE FOR SAVING MAY BE MODIFIED; AS GIVEN HERE,THE C OUTPUT IS IN CARD IMAGE FORM WITH A FORMAT OF (4F20.14). C C B) PRINTED OUTPUT IS A LISTING OF TRNFCT ON UNIT 6, IN C THE FORMAT (1X,8F15.8) (FORMAT STATEMENT 301). THE USER MAY C MODIFY OR ELIMINATE THIS AT WILL. C C ************ FUNCTION INTERPOLATER ROUTINE ***************** C C C ****** THE FOLLOWING PARAMETER GIVES THE NUMBER OF ******* C ****** DATA LEVELS IN THE MODEL ******* C **************************************************************** INCLUDE "parmco2" PARAMETER (NLEVLS=LERIC) C **************************************************************** PARAMETER (NLP1=NLEVLS+1,NLP2=NLEVLS+2) COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N COMMON/PRESS/PA(109) COMMON/TRAN/ TRANSA(109,109) COMMON / OUTPUT / TRNS(NLP1,NLP1) COMMON/INPUTP/P(NLP1),PD(NLP2) DIMENSION PS(NLP1),PDS(NLP2),PLM(NLP1) DIMENSION NRTAB(3) DIMENSION T15A(NLP2,2),T15B(NLP1) DIMENSION T22(NLP1,NLP1,3) DATA NRTAB/1,2,4/ C*********************************** C THE FOLLOWING ARE THE INPUT FORMATS C-CP 100 FORMAT (F20.14) 201 FORMAT (5E16.9) 202 FORMAT (I3) CO222 203 FORMAT (F12.6,I2) 203 FORMAT (F12.6) C THE FOLLOWING ARE THE OUTPUT FORMATS 102 FORMAT (4F20.14) 301 FORMAT (1X,8F15.8) C CCC REWIND 15 CCC REWIND 20 CNOV89 REWIND ITAPE CNOV89 CCC REWIND 22 C C CALCULATION OF PA -THE -TABLE- OF 109 GRID PRESSURES C NOTE-THIS CODE MUST NOT BE CHANGED BY THE USER^^^^^^^^^ PA(1)=0. FACT15=10.**(1./15.) FACT30=10.**(1./30.) PA(2)=1.0E-3 DO 231 I=2,76 PA(I+1)=PA(I)*FACT15 231 CONTINUE DO 232 I=77,108 PA(I+1)=PA(I)*FACT30 232 CONTINUE C N=25 NLV=NLEVLS NLP1V=NLP1 NLP2V=NLP2 C READ IN THE CO2 MIXING RATIO(IN UNITS OF 330 PPMV),AND AN INDEX C GIVING THE FREQUENCY RANGE OF THE LBL DATA CO222 READ (5,203) RATIO,IR CCC IR = 1 CCC READ (5,203) RATIO CO222 *********************************** C***VALUES FOR IR***** C IR=1 CONSOL. LBL TRANS. =490-850 C IR=2 CONSOL. LBL TRANS. =490-670 C IR=3 CONSOL. LBL TRANS. =670-850 C IR=4 CONSOL. LBL TRANS. =2270-2380 C*** IR MUST BE 1,2,3 OR 4 FOR THE PGM. TO WORK C ALSO READ IN THE METHOD NO.(1 OR 2) CCC READ (5,202) NMETHD IF (RATIO.EQ.1.0) GO TO 621 IF (RATIO.GE.2.0) GO TO 622 IF (RATIO.EQ.4.0) GO TO 623 IF (RATIO.GT.1.0.AND.RATIO.LT.2.0) GO TO 624 IF (RATIO.GT.2.0.AND.RATIO.LT.4.0) GO TO 625 STOP 8746 CNOV89 621 ITAP1=20 621 ITAP1=ITAPE CNOV89 NTAP=1 GO TO 630 622 ITAP1=21 NTAP=1 GO TO 630 623 ITAP1=22 NTAP=1 GO TO 630 CNOV89 624 ITAP1=20 624 ITAP1=ITAPE CNOV89 NTAP=2 RATSTD=2.0 RATSM=1.0 GO TO 630 625 ITAP1=21 NTAP=2 RATSTD=4.0 RATSM=2.0 630 CONTINUE IF (NMETHD.EQ.2) GO TO 502 C *****CARDS FOR PURPOSE 1(NMETHD=1) CCC READ (15,201) (P(I),I=1,NLP1) DO 300 I=1,NLP1 P(I)=T15B(I) 300 CONTINUE DO 801 I=1,NLP1 PS(I)=P(I) 801 CONTINUE GO TO 503 502 CONTINUE C *****CARDS FOR PURPOSE 2(NMETHD=2) CCC READ (15,201) (PD(I),I=1,NLP2) CCC READ (15,201) (PLM(I),I=1,NLP1) DO 303 I=1,NLP2 PD(I)=T15A(I,1) 303 CONTINUE DO 302 I=1,NLP1 PLM(I)=T15A(I,2) 302 CONTINUE DO 802 I=1,NLP1 PDS(I)=PD(I+1) PS(I)=PLM(I) 802 CONTINUE C 503 CONTINUE C *****DO LOOP CONTROLLING NUMBER OF OUTPUT MATRICES CNOV89 CNOV89 DO 400 KKK=1,3 ICLOOP = 3 IF (IR.EQ.4) ICLOOP = 1 DO 400 KKK=1,ICLOOP CNOV89 IF (NTAP.EQ.2) CALL RCTRNS(ITAP1,RATSTD,RATSM,RATIO) C ********************** IF (NMETHD.EQ.2) GO TO 505 C *****CARDS FOR PURPOSE 1(NMETHD=1) DO 803 I=1,NLP1 P(I)=PS(I) 803 CONTINUE GO TO 506 505 CONTINUE C *****CARDS FOR PURPOSE 2(NMETHD=2) DO 804 I=1,NLP1 PD(I)=PDS(I) P(I)=PS(I) 804 CONTINUE C 506 CONTINUE IA=108 IAP=IA+1 CNOV89 IF (NTAP.EQ.1) READ (20,100) ((TRANSA(I,J),I=1,109),J=1,109) C-CP IF (NTAP.EQ.1) READ (ITAPE,100) ((TRANSA(I,J),I=1,109),J=1,109) IF (NTAP.EQ.1) THEN DO J=1,109 DO I=1,109 READ(ITAPE,100) TRANSA(I,J) END DO END DO ENDIF CNOV89 DO 4 I=1,IAP TRANSA(I,I)=1.0 4 CONTINUE CALL COEINT(RATIO,IR) DO 805 I=1,NLP1 DO 805 J=1,NLP1 TRNS(J,I)=1.00 805 CONTINUE DO 10 I=1,NLP1 DO 20 J=1,I IF (I.EQ.J) GO TO 20 P1=P(J) P2=P(I) CALL SINTR2 TRNS(J,I)=TRNSLO 20 CONTINUE 10 CONTINUE DO 47 I=1,NLP1 DO 47 J=I,NLP1 TRNS(J,I)=TRNS(I,J) 47 CONTINUE C *****THIS IS THE END OF PURPOSE 1 CALCULATIONS IF (NMETHD.EQ.1) GO TO 2872 C DO 51 J=1,NLP1 DO 52 I=2,NLP1 IA=I JA=J N=25 IF (I.NE.J) N=3 CALL QUADSR(NLV,NLP1V,NLP2V,P,PD,TRNS) 52 CONTINUE 51 CONTINUE C *****THIS IS THE END OF PURPOSE 2 CALCULATIONS 2872 CONTINUE C C+++ WRITE (6,301) ((TRNS(I,J),I=1,NLP1),J=1,NLP1) CCC WRITE (22,102) ((TRNS(I,J),I=1,NLP1),J=1,NLP1) DO 304 J=1,NLP1 DO 304 I=1,NLP1 T22(I,J,KKK) = TRNS(I,J) 304 CONTINUE 400 CONTINUE RETURN END CCCC PROGRAM CO2IN1 SUBROUTINE CO2IN1(T20,T21,T66,IQ) C CO2IN1=CO2INS FOR METHOD 1 C ********************************************************* C SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 *** C ..... K.CAMPANA MARCH 1988,OCTOBER 1988 C ..... K.CAMPANA DECEMBER 88 CLEANED UP FOR LAUNCHER C ********************************************************* INCLUDE "parmco2" PARAMETER (L=LERIC,LP1=L+1) DIMENSION T20(LP1,LP1,3),T21(LP1,LP1,3),T66(L,6) DIMENSION DCDT8(LP1,LP1),DCDT10(LP1,LP1),CO2PO(LP1,LP1), * CO2800(LP1,LP1),CO2PO1(LP1,LP1),CO2801(LP1,LP1),CO2PO2(LP1,LP1), * CO2802(LP1,LP1),N(LP1),D2CT8(LP1,LP1),D2CT10(LP1,LP1) ITIN=20 ITIN1=21 CO222 LATEST CODE HAS IQ=1 CCC IQ=4 1011 FORMAT (4F20.14) CCC READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1) CCC READ (ITIN1,1011) ((CO2802(I,J),I=1,LP1),J=1,LP1) DO 300 J=1,LP1 DO 300 I=1,LP1 CO2PO(I,J) = T20(I,J,1) CNOV89 IF (IQ.EQ.5) GO TO 300 CNOV89 CO2PO1(I,J) = T20(I,J,2) CO2PO2(I,J) = T20(I,J,3) 300 CONTINUE DO 301 J=1,LP1 DO 301 I=1,LP1 CO2800(I,J) = T21(I,J,1) CNOV89 IF (IQ.EQ.5) GO TO 301 CNOV89 CO2801(I,J) = T21(I,J,2) CO2802(I,J) = T21(I,J,3) 301 CONTINUE C***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS C ARE: C IQ=1 560-800 (CONSOL.=490-850) C IQ=2 560-670 (CONSOL.=490-670) C IQ=3 670-800 (CONSOL.=670-850) C IQ=4 560-760 (ORIGINAL CODE) (CONSOL.=490-850) CNOV89 C IQ=5 2270-2380 (CONSOL.=2270-2380) CNOV89 C THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS C USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT C WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S. IF (IQ.EQ.1) THEN C1=1.5 C2=0.5 ENDIF IF (IQ.EQ.2) THEN C1=18./11. C2=7./11. ENDIF IF (IQ.EQ.3) THEN C1=18./13. C2=5./13. ENDIF IF (IQ.EQ.4) THEN C1=1.8 C2=0.8 ENDIF CNOV89 IF (IQ.EQ.5) THEN C1=1.0 C2=0.0 ENDIF CNOV89 DO 1021 I=1,LP1 DO 1021 J=1,LP1 CO2PO(J,I)=C1*CO2PO(J,I)-C2 CO2800(J,I)=C1*CO2800(J,I)-C2 CNOV89 IF (IQ.EQ.5) GO TO 1021 CNOV89 CO2PO1(J,I)=C1*CO2PO1(J,I)-C2 CO2801(J,I)=C1*CO2801(J,I)-C2 CO2PO2(J,I)=C1*CO2PO2(J,I)-C2 CO2802(J,I)=C1*CO2802(J,I)-C2 1021 CONTINUE CNOV89 IF (IQ.GE.1.AND.IQ.LE.4) THEN CNOV89 DO 1 J=1,LP1 DO 1 I=1,LP1 DCDT8(I,J)=.02*(CO2801(I,J)-CO2802(I,J))*100. DCDT10(I,J)=.02*(CO2PO1(I,J)-CO2PO2(I,J))*100. D2CT8(I,J)=.0016*(CO2801(I,J)+CO2802(I,J)-2.*CO2800(I,J))*1000. D2CT10(I,J)=.0016*(CO2PO1(I,J)+CO2PO2(I,J)-2.*CO2PO(I,J))*1000. 1 CONTINUE CNOV89 ENDIF CNOV89 CO222 ********************************************************* CCC REWIND 66 C SAVE CDTM51,CO2M51,C2DM51,CDTM58,CO2M58,C2DM58..ON TEMPO FILE CCC WRITE (66) (DCDT10(I,I+1),I=1,L) CCC WRITE (66) (CO2PO(I,I+1),I=1,L) CCC WRITE (66) (D2CT10(I,I+1),I=1,L) CCC WRITE (66) (DCDT8(I,I+1),I=1,L) CCC WRITE (66) (CO2800(I,I+1),I=1,L) CCC WRITE (66) (D2CT8(I,I+1),I=1,L) CCC REWIND 66 CO222 ********************************************************* DO 400 I=1,L T66(I,2) = CO2PO(I,I+1) T66(I,5) = CO2800(I,I+1) CNOV89 IF (IQ.EQ.5) GO TO 400 CNOV89 T66(I,1) = DCDT10(I,I+1) T66(I,3) = D2CT10(I,I+1) T66(I,4) = DCDT8(I,I+1) T66(I,6) = D2CT8(I,I+1) 400 CONTINUE RETURN END CCCC PROGRAM PTZ - COURTESY OF DAN SCHWARZKOPF,GFDL DEC 1987.... SUBROUTINE CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU,LREAD) INCLUDE "parmco2" PARAMETER (NL=LERIC,NLP=NL+1,NLP2=NL+2) C C ** THIS PROGRAM CALCULATES TEMPERATURES ,H2O MIXING RATIOS C ** AND O3 MIXING RATIOS BY USING AN ANALYTICAL C ** FUNCTION WHICH APPROXIMATES C ** THE US STANDARD (1976). THIS IS C ** CALCULATED IN FUNCTION 'ANTEMP', WHICH IS CALLED BY THE C ** MAIN PROGRAM. THE FORM OF THE ANALYTICAL FUNCTION WAS C ** SUGGESTED TO ME IN 1971 BY RICHARD S. LINDZEN. C ****************************************************************** C CODE TO SAVE STEMP,GTEMP ON DATA SET,BRACKETED BY CO222 ** C ....K. CAMPANA MARCH 88,OCTOBER 88 DIMENSION SGTEMP(NLP,2),T41(NLP2,2),T42(NLP), 1 T43(NLP2,2),T44(NLP) DIMENSION SGLVNU(NLP),SIGLNU(NL) C ****************************************************************** C C*****THIS VERSION IS ONLY USABLE FOR 1976 US STD ATM AND OBTAINS C QUANTITIES FOR CO2 INTERPOLATION AND INSERTION INTO OPERA- C TIONAL RADIATION CODES C CHARACTER*20 PROFIL DIMENSION PRESS(NLP),TEMP(NLP),ALT(NLP),WMIX(NLP),O3MIX(NLP) DIMENSION WMXINT(NLP,4),WMXOUT(NLP2),OMXINT(NLP,4),OMXOUT(NLP2) DIMENSION PD(NLP2),GTEMP(NLP) DIMENSION PRS(NLP),TEMPS(NLP),PRSINT(NLP),TMPINT(NLP,4),A(NLP,4) DIMENSION PROUT(NLP2),TMPOUT(NLP2),TMPFLX(NLP2),TMPMID(NLP2) C C DATA PROFIL/ 6 'US STANDARD 1976'/ DATA PSMAX/1013.250/ C C ** NTYPE IS AN INTEGER VARIABLE WHICH HAS THE FOLLOWING C ** VALUES: 0 =SIGMA LEVELS ARE USED; 1= SKYHI L40 LEVELS C ** ARE USED; 2 = SKYHI L80 LEVELS ARE USED. DEFAULT: 0 C NTYPE=0 CO222 READ (*,*) NTYPE 5 NLEV=NL DELZAP=0.5 R=8.31432 G0=9.80665 ZMASS=28.9644 AA=6356.766 ALT(1)=0.0 TEMP(1)=ANTEMP(6,0.0) C*******DETERMINE THE PRESSURES (PRESS) PSTAR=PSMAX C CO222 IF (NTYPE.EQ.1) CALL SKYP(PSTAR,PD,GTEMP) CO222 IF (NTYPE.EQ.2) CALL SKY80P(PSTAR,PD,GTEMP) CO222 IF (NTYPE.EQ.0) CALL SIGP(PSTAR,PD,GTEMP) CCC---- CALL SIGP(PSTAR,PD,GTEMP) CALL SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU,LREAD) PD(NLP2)=PSTAR DO 40 N=1,NLP PRSINT(N)=PD(NLP2+1-N) 40 CONTINUE C *** CALCULATE TEMPS FOR SEVERAL PRESSURES TO DO QUADRATURE DO 504 NQ=1,4 DO 505 N=2,NLP 505 PRESS(N)=PRSINT(N)+0.25*(NQ-1)*(PRSINT(N-1)-PRSINT(N)) PRESS(1)=PRSINT(1) C********************* DO 100 N=1,NLEV C C ** ESTABLISH COMPUTATATIONAL LEVELS BETWEEN USER LEVELS AT C ** INTERVALS OF APPROXIMATELY 'DELZAP' KM. C DLOGP=7.0*ALOG(PRESS(N)/PRESS(N+1)) NINT=DLOGP/DELZAP NINT=NINT+1 ZNINT=NINT G=G0 DZ=R*DLOGP/(7.0*ZMASS*G*ZNINT) HT=ALT(N) C C ** CALCULATE HEIGHT AT NEXT USER LEVEL BY MEANS OF C ** RUNGE-KUTTA INTEGRATION. C DO 200 M=1,NINT RK1=ANTEMP(6,HT)*DZ RK2=ANTEMP(6,HT+0.5*RK1)*DZ RK3=ANTEMP(6,HT+0.5*RK2)*DZ RK4=ANTEMP(6,HT+RK3)*DZ HT=HT+0.16666667*(RK1+RK2+RK2+RK3+RK3+RK4) 200 CONTINUE ALT(N+1)=HT TEMP(N+1)=ANTEMP(6,HT) 100 CONTINUE DO 506 N=1,NLP TMPINT(N,NQ)=TEMP(N) A(N,NQ)=ALT(N) 506 CONTINUE DO 523 N=1,NLP CALL MIXRAT(6,PRESS(N),TMPINT(N,NQ),WMXINT(N,NQ)) WMXINT(N,NQ)=WMXINT(N,NQ)*1.0E-6 CALL OZONE(6,PRESS(N),TMPINT(N,NQ),OMXINT(N,NQ)) 523 CONTINUE 504 CONTINUE C ***APPLY SIMPSON S RULE DO 507 N=2,NLP TEMP(N)=1./12.*(TMPINT(N-1,1)+TMPINT(N,1)+4.*TMPINT(N,2)+ 1 2.*TMPINT(N,3)+4.*TMPINT(N,4)) WMIX(N)=1./12.*(WMXINT(N-1,1)+WMXINT(N,1)+4.*WMXINT(N,2)+ 1 2.*WMXINT(N,3)+4.*WMXINT(N,4)) O3MIX(N)=1./12.*(OMXINT(N-1,1)+OMXINT(N,1)+4.*OMXINT(N,2)+ 1 2.*OMXINT(N,3)+4.*OMXINT(N,4)) 507 CONTINUE C***OUTPUT FOR LINE-BY-LINE CALCS TMPOUT(1)=TMPINT(NLP,1) WMXOUT(1)=WMXINT(NLP,1) OMXOUT(1)=OMXINT(NLP,1) TMPMID(1)=TMPINT(NLP,1) DO 520 I=2,NLP2 TMPOUT(I)=TEMP(NLP2+1-I) WMXOUT(I)=WMIX(NLP2+1-I) OMXOUT(I)=O3MIX(NLP2+1-I) TMPMID(I)=TMPINT(NLP2+1-I,3) 520 CONTINUE DO 521 I=1,NLP2 PROUT(I)=PD(I) 521 CONTINUE DO 5221 I=2,NLP2 TMPFLX(I)=TMPINT(NLP2+1-I,1) 5221 CONTINUE TMPFLX(1)=TMPFLX(2) C C ** CALCULATE WATER MIXING RATIO USING LUTHER PROGRAM C DO 300 N=1,NLP CALL MIXRAT(6,PRSINT(N),TMPINT(N,1),WMIX(N)) WMIX(N)=WMIX(N)*1.0E-06 300 CONTINUE C C CALCULATE OZONE MIXING RATIO USING SCHWARZKOPF PROGRAM C DO 400 N=1,NLP CALL OZONE(6,PRSINT(N),TMPINT(N,1),O3MIX(N)) 400 CONTINUE C C+++ WRITE (6,101) PROFIL 101 FORMAT (1X,A20) C+++ WRITE (6,201) 201 FORMAT(5X,' HEIGHT TEMPERATURE PRESSURE R(H2O) ', 1' R(O3)') C+++ WRITE (6,202) A(1,1),TMPINT(1,1),PRSINT(1),WMIX(1),O3MIX(1) DO 210 N=2,NLP C+++ WRITE (6,203) TMPOUT(NLP2+1-N),WMXOUT(NLP2+1-N),OMXOUT(NLP2+1-N) C+++ WRITE (6,202) A(N,1),TMPINT(N,1),PRSINT(N),WMIX(N),O3MIX(N) 210 CONTINUE 203 FORMAT (1X,14X,F14.6,14X,2E14.6) 202 FORMAT(1X,2F14.6,E14.6,2E14.6) C DETERMINE WHAT FORMATS TO USE*** NREP=NLP/5 NREM=NLP-5*NREP IF (NREM.EQ.0) THEN NREP=NREP-1 NREM=5 ENDIF NREP2=NL/5 NREM2=NL-5*NREP2 IF (NREM2.EQ.0) THEN NREP2=NREP2-1 NREM2=5 ENDIF NREP3=NLP/4 NREM3=NLP-4*NREP3 IF (NREM3.EQ.0) THEN NREP3=NREP3-1 NREM3=4 ENDIF CO222 ***************************************************** CC REWIND 66 CO222 ***************************************************** C***OUTPUT TEMPERATURES DO 800 IOUT=1,2 IF (IOUT.EQ.1) THEN WRITE (16,701) ELSE WRITE (16,702) ENDIF 701 FORMAT (6X,'DATA DTEMP /') 702 FORMAT (6X,'DATA STEMP /') NF=0 CO222 ***************************************************** C SAVE STEMP CC IF(IOUT.EQ.2) WRITE(66) (TMPINT(NLP2-N,1),N=1,NLP) DO 901 N=1,NLP SGTEMP(N,1) = TMPINT(NLP2-N,1) 901 CONTINUE CO222 ***************************************************** IF (NREP.NE.0) THEN DO 801 NR=1,NREP NS=5*(NR-1)+1 NF=NS+4 WRITE (16,656) (TMPINT(NLP2-N,1),N=NS,NF) 801 CONTINUE ENDIF NF=NF+1 IF (NREM.EQ.1) THEN WRITE (16,616) (TMPINT(NLP2-N,1),N=NF,NLP) ENDIF IF (NREM.EQ.2) THEN WRITE (16,626) (TMPINT(NLP2-N,1),N=NF,NLP) ENDIF IF (NREM.EQ.3) THEN WRITE (16,636) (TMPINT(NLP2-N,1),N=NF,NLP) ENDIF IF (NREM.EQ.4) THEN WRITE (16,646) (TMPINT(NLP2-N,1),N=NF,NLP) ENDIF IF (NREM.EQ.5) THEN WRITE (16,606) (TMPINT(NLP2-N,1),N=NF,NLP) ENDIF 616 FORMAT (5X,'*',1X,F12.6,'/') 626 FORMAT (5X,'*',1X,F12.6,',',F12.6,'/') 636 FORMAT (5X,'*',1X,F12.6,',',F12.6,',',F12.6,'/') 646 FORMAT (5X,'*',1X,F12.6,',',F12.6,',',F12.6,',', 1 F12.6,'/') 656 FORMAT (5X,'*',1X,F12.6,',',F12.6,',',F12.6,',', 1 F12.6,',',F12.6,',') 606 FORMAT (5X,'*',1X,F12.6,',',F12.6,',',F12.6,',', 1 F12.6,',',F12.6,'/') 800 CONTINUE C***OUTPUT GTEMP CO222 ***************************************************** C SAVE GTEMP CC WRITE(66) (GTEMP(N),N=1,NLP) DO 902 N=1,NLP SGTEMP(N,2) = GTEMP(N) 902 CONTINUE CO222 ***************************************************** WRITE (16,706) 706 FORMAT (6X,'DATA GTEMP /') NF=0 IF (NREP3.NE.0) THEN DO 805 NR=1,NREP3 NS=4*(NR-1)+1 NF=NS+3 WRITE (16,648) (GTEMP(N),N=NS,NF) 805 CONTINUE ENDIF NF=NF+1 IF (NREM3.EQ.1) THEN WRITE (16,618) (GTEMP(N),N=NF,NLP) ENDIF IF (NREM3.EQ.2) THEN WRITE (16,628) (GTEMP(N),N=NF,NLP) ENDIF IF (NREM3.EQ.3) THEN WRITE (16,638) (GTEMP(N),N=NF,NLP) ENDIF IF (NREM3.EQ.4) THEN WRITE (16,608) (GTEMP(N),N=NF,NLP) ENDIF C***OUTPUT WMIX IN 10**2 GM/GM DO 630 N=2,NLP WMIX(N)=WMIX(N)*1.0E2 630 CONTINUE NF=0 WRITE (16,703) 703 FORMAT (6X,'DATA RR /') IF (NREP2.NE.0) THEN DO 802 NR=1,NREP2 NS=5*(NR-1)+1 NF=NS+4 WRITE (16,657) (WMIX(NLP2-N),N=NS,NF) 802 CONTINUE ENDIF NF=NF+1 IF (NREM2.EQ.1) THEN WRITE (16,617) (WMIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.2) THEN WRITE (16,627) (WMIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.3) THEN WRITE (16,637) (WMIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.4) THEN WRITE (16,647) (WMIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.5) THEN WRITE (16,607) (WMIX(NLP2-N),N=NF,NL) ENDIF 617 FORMAT (5X,'*',1X,E12.6,'/') 627 FORMAT (5X,'*',1X,E12.6,',',E12.6,'/') 637 FORMAT (5X,'*',1X,E12.6,',',E12.6,',',E12.6,'/') 647 FORMAT (5X,'*',1X,E12.6,',',E12.6,',',E12.6,',', 1 E12.6,'/') 657 FORMAT (5X,'*',1X,E12.6,',',E12.6,',',E12.6,',', 1 E12.6,',',E12.6,',') 607 FORMAT (5X,'*',1X,E12.6,',',E12.6,',',E12.6,',', 1 E12.6,',',E12.6,'/') C***OUTPUT O3MIX IN 10**3 GM/GM DO 631 N=2,NLP O3MIX(N)=O3MIX(N)*1.0E3 631 CONTINUE NF=0 WRITE (16,704) 704 FORMAT (6X,'DATA QQO3 /') IF (NREP2.NE.0) THEN DO 803 NR=1,NREP2 NS=5*(NR-1)+1 NF=NS+4 WRITE (16,657) (O3MIX(NLP2-N),N=NS,NF) 803 CONTINUE ENDIF NF=NF+1 IF (NREM2.EQ.1) THEN WRITE (16,617) (O3MIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.2) THEN WRITE (16,627) (O3MIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.3) THEN WRITE (16,637) (O3MIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.4) THEN WRITE (16,647) (O3MIX(NLP2-N),N=NF,NL) ENDIF IF (NREM2.EQ.5) THEN WRITE (16,607) (O3MIX(NLP2-N),N=NF,NL) ENDIF C***OUTPUT PROUT IN CGS DO 629 N=2,NLP2 PROUT(N)=PROUT(N)*1.0E3 629 CONTINUE WRITE (16,705) 705 FORMAT (6X,'DATA PPRESS /') NF=0 IF (NREP3.NE.0) THEN DO 804 NR=1,NREP3 NS=4*(NR-1)+2 NF=NS+3 WRITE (16,648) (PROUT(N),N=NS,NF) 804 CONTINUE ENDIF NF=NF+1 IF (NREM3.EQ.1) THEN WRITE (16,618) (PROUT(N),N=NF,NLP2) ENDIF IF (NREM3.EQ.2) THEN WRITE (16,628) (PROUT(N),N=NF,NLP2) ENDIF IF (NREM3.EQ.3) THEN WRITE (16,638) (PROUT(N),N=NF,NLP2) ENDIF IF (NREM3.EQ.4) THEN WRITE (16,608) (PROUT(N),N=NF,NLP2) ENDIF 618 FORMAT (5X,'*',1X,E15.9,'/') 628 FORMAT (5X,'*',1X,E15.9,',',E15.9,'/') 638 FORMAT (5X,'*',1X,E15.9,',',E15.9,',',E15.9,'/') 648 FORMAT (5X,'*',1X,E15.9,',',E15.9,',',E15.9,',',E15.9,',') 608 FORMAT (5X,'*',1X,E15.9,',',E15.9,',',E15.9,',',E15.9,'/') RETURN END SUBROUTINE MIXRAT(M,P,T,Q) DIMENSION A(6,7),B(6,7),QS(7),PRAT(6,7),PS(7) DIMENSION PMIN(7),QCONST(7) C TO MAKE PGM. FLOW,WE ASSUME THE US STD ATM.=SUBTROPICAL SUMMER IN C WATER AND OZONE DATA QS/19.,14.,3.5,9.1,1.2,2*14./ DATA A/3.55,1.37,3.99,2.79,8.66,4.66, 1 3.31,4.77,4.19,1.30,8.21,4.52, 2 2.72,2.11,3.41,5.44,3.91,0.00, 3 3.76,2.95,0.83,4.71,7.43,4.28, 4 0.00,-.59,4.99,4.18,3.30,0.00, 5 3.31,4.77,4.19,1.30,8.21,4.52, 5 3.31,4.77,4.19,1.30,8.21,4.52/ DATA B/1.91,-7.58,-0.41,-1.52, 1.81, 0.00, 1 -1.83, 0.67, 0.00,-2.26, 1.67, 0.00, 2 0.39,-2.02,-0.81, 0.71, 0.00, 0.00, 3 2.37,-1.01,-3.17, 0.00, 1.41, 0.00, 4 0.00,-4.46, 0.86, 0.38, 0.00, 0.00, 5 -1.83, 0.67, 0.00,-2.26, 1.67, 0.00, 5 -1.83, 0.67, 0.00,-2.26, 1.67, 0.00/ DATA PRAT/0.795 ,0.694, 0.339, 0.172, 0.110, 0.000, 1 0.558, 0.421, 0.278, 0.172, 0.110, 0.000, 2 0.776, 0.342, 0.263, 0.100, 0.000, 0.000, 3 0.787, 0.375, 0.294, 0.145, 0.107, 0.000, 4 0.876, 0.350, 0.185, 0.100, 0.000, 0.000, 5 0.558, 0.421, 0.278, 0.172, 0.110, 0.000, 5 0.558, 0.421, 0.278, 0.172, 0.110, 0.000/ DATA PS/1013.01, 1013.01,1018.01,1010.01,1013.01,2*1013.01/ DATA PMIN/110.99,152.99,117.79,124.99,110.29,2*152.99/ DATA QCONST/3.25,4.0,4.0,4.0,4.0,2*4.0/ CONST=10.0*8314.320/28.9644 PR=P/PS(M) N=1 15 IF (PR.GE.PRAT(N,M)) GO TO 20 IF (P.LE.PMIN(M)) GO TO 25 N=N+1 GO TO 15 20 ETA=A(N,M)+B(N,M)*ALOG(PR) Q=CONST*(T/P)*QS(M)*PR**ETA GO TO 40 25 Q=QCONST(M) GO TO 40 40 RETURN END SUBROUTINE OZONE(L,P,T,R) C L=NO. OF PROFILE;P=USER PRESSURE;R=OUTPUT O3 MIXING RATIO C C TO MAKE PGM FLOW,WE ASSUME US STD ATM=SUBTROPICAL SUMMER IN C WATER AND OZONE DIMENSION O3MAC(33,7),PMAC(33,7) DIMENSION O3MAC1(33),O3MAC2(33),O3MAC3(33),O3MAC4(33),O3MAC5(33), 1 O3MAC6(33),O3MAC7(33) DIMENSION PMAC1(33),PMAC2(33),PMAC3(33),PMAC4(33),PMAC5(33), 1 PMAC6(33),PMAC7(33) EQUIVALENCE (O3MAC1(1),O3MAC(1,1)),(O3MAC2(1),O3MAC(1,2)) EQUIVALENCE (O3MAC3(1),O3MAC(1,3)),(O3MAC4(1),O3MAC(1,4)) EQUIVALENCE (O3MAC5(1),O3MAC(1,5)),(PMAC1(1),PMAC(1,1)) EQUIVALENCE (PMAC2(1),PMAC(1,2)),(PMAC3(1),PMAC(1,3)) EQUIVALENCE (PMAC4(1),PMAC(1,4)),(PMAC5(1),PMAC(1,5)) EQUIVALENCE (O3MAC6(1),O3MAC(1,6)),(PMAC6(1),PMAC(1,6)) EQUIVALENCE (O3MAC7(1),O3MAC(1,7)),(PMAC7(1),PMAC(1,7)) C O3MAC=MC CLATCHY OZONE MIXING RATIO;P=MC CLARTCHY PRESSURESQQ DATA O3MAC1/ 1 5.6E-5,5.6E-5,5.4E-5,5.1E-5,4.7E-5,4.5E-5,4.3E-5,4.1E-5,3.9E-5, 2 3.9E-5,3.9E-5,4.1E-5,4.3E-5,4.5E-5,4.5E-5,4.7E-5,4.7E-5,6.9E-5, 3 9.0E-5,1.4E-4,1.9E-4,2.4E-4,2.8E-4,3.2E-4,3.4E-4,3.4E-4,2.4E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC2/ 1 6.0E-5,6.0E-5,6.0E-5,6.2E-5,6.4E-5,6.6E-5,6.9E-5,7.5E-5,7.9E-5, 2 8.6E-5,9.0E-5,1.1E-4,1.2E-4,1.5E-4,1.8E-4,1.9E-4,2.1E-4,2.4E-4, 3 2.8E-4,3.2E-4,3.4E-4,3.6E-4,3.6E-4,3.4E-4,3.2E-4,3.0E-4,2.0E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC3/ 1 6.0E-5,5.4E-5,4.9E-5,4.9E-5,4.9E-5,5.8E-5,6.4E-5,7.7E-5,9.0E-5, 2 1.2E-4,1.6E-4,2.1E-4,2.6E-4,3.0E-4,3.2E-4,3.4E-4,3.6E-4,3.9E-4, 3 4.1E-4,4.3E-4,4.5E-4,4.3E-4,4.3E-4,3.9E-4,3.6E-4,3.4E-4,1.9E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC4/ 1 4.9E-5,5.4E-5,5.6E-5,5.8E-5,6.0E-5,6.4E-5,7.1E-5,7.5E-5,7.9E-5, 2 1.1E-4,1.3E-4,1.8E-4,2.1E-4,2.6E-4,2.8E-4,3.2E-4,3.4E-4,3.9E-4, 3 4.1E-4,4.1E-4,3.9E-4,3.6E-4,3.2E-4,3.0E-4,2.8E-4,2.6E-4,1.4E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC5/ 1 4.1E-5,4.1E-5,4.1E-5,4.3E-5,4.5E-5,4.7E-5,4.9E-5,7.1E-5,9.0E-5, 2 1.6E-4,2.4E-4,3.2E-4,4.3E-4,4.7E-4,4.9E-4,5.6E-4,6.2E-4,6.2E-4, 3 6.2E-4,6.0E-4,5.6E-4,5.1E-4,4.7E-4,4.3E-4,3.6E-4,3.2E-4,1.5E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC6/ 1 6.0E-5,6.0E-5,6.0E-5,6.2E-5,6.4E-5,6.6E-5,6.9E-5,7.5E-5,7.9E-5, 2 8.6E-5,9.0E-5,1.1E-4,1.2E-4,1.5E-4,1.8E-4,1.9E-4,2.1E-4,2.4E-4, 3 2.8E-4,3.2E-4,3.4E-4,3.6E-4,3.6E-4,3.4E-4,3.2E-4,3.0E-4,2.0E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA O3MAC7/ 1 6.0E-5,6.0E-5,6.0E-5,6.2E-5,6.4E-5,6.6E-5,6.9E-5,7.5E-5,7.9E-5, 2 8.6E-5,9.0E-5,1.1E-4,1.2E-4,1.5E-4,1.8E-4,1.9E-4,2.1E-4,2.4E-4, 3 2.8E-4,3.2E-4,3.4E-4,3.6E-4,3.6E-4,3.4E-4,3.2E-4,3.0E-4,2.0E-4, 4 9.2E-5,4.1E-5,1.3E-5,4.3E-6,8.6E-8,4.3E-11 5 / DATA PMAC1/ 1 1.013E+3,9.040E+2,8.050E+2,7.150E+2,6.330E+2,5.590E+2,4.920E+2, 2 4.320E+2,3.780E+2,3.290E+2,2.860E+2,2.470E+2,2.130E+2,1.820E+2, 3 1.560E+2,1.320E+2,1.110E+2,9.370E+1,7.890E+1,6.660E+1,5.650E+1, 4 4.800E+1,4.090E+1,3.500E+1,3.000E+1,2.570E+1,1.220E+1,6.000E+0, 5 3.050E+0,1.590E+0,8.540E-1,5.790E-2,3.000E-4/ DATA PMAC2/ 1 1.013E+3,9.020E+2,8.020E+2,7.100E+2,6.280E+2,5.540E+2,4.870E+2, 2 4.260E+2,3.720E+2,3.240E+2,2.810E+2,2.430E+2,2.090E+2,1.790E+2, 3 1.530E+2,1.300E+2,1.110E+2,9.500E+1,8.120E+1,6.950E+1,5.950E+1, 4 5.100E+1,4.370E+1,3.760E+1,3.220E+1,2.770E+1,1.320E+1,6.520E+0, 5 3.330E+0,1.760E+0,9.510E-1,6.710E-2,3.000E-4/ DATA PMAC3/ 1 1.018E+3,8.973E+2,7.897E+2,6.938E+2,6.081E+2,5.313E+2,4.627E+2, 2 4.016E+2,3.473E+2,2.992E+2,2.568E+2,2.199E+2,1.882E+2,1.610E+2, 3 1.378E+2,1.178E+2,1.007E+2,8.610E+1,7.350E+1,6.280E+1,5.370E+1, 4 4.580E+1,3.910E+1,3.340E+1,2.860E+1,2.430E+1,1.110E+1,5.180E+0, 5 2.530E+0,1.290E+0,6.820E-1,4.670E-2,3.000E-4/ DATA PMAC4/ 1 1.010E+3,8.960E+2,7.929E+2,7.000E+2,6.160E+2,5.410E+2,4.730E+2, 2 4.130E+2,3.590E+2,3.107E+2,2.677E+2,2.300E+2,1.977E+2,1.700E+2, 3 1.460E+2,1.250E+2,1.080E+2,9.280E+1,7.980E+1,6.860E+1,5.890E+1, 4 5.070E+1,4.360E+1,3.750E+1,3.227E+1,2.780E+1,1.340E+1,6.610E+0, 5 3.400E+0,1.810E+0,9.870E-1,7.070E-2,3.000E-4/ DATA PMAC5/ 1 1.013E+3,8.878E+2,7.775E+2,6.798E+2,5.932E+2,5.158E+2,4.467E+2, 2 3.853E+2,3.308E+2,2.829E+2,2.418E+2,2.067E+2,1.766E+2,1.510E+2, 3 1.291E+2,1.103E+2,9.431E+1,8.058E+1,6.882E+1,5.875E+1,5.014E+1, 4 4.277E+1,3.647E+1,3.109E+1,2.649E+1,2.256E+1,1.020E+1,4.701E+0, 5 2.243E+0,1.113E+0,5.719E-1,4.016E-2,3.000E-4/ DATA PMAC6/ 1 1.013E+3,9.020E+2,8.020E+2,7.100E+2,6.280E+2,5.540E+2,4.870E+2, 2 4.260E+2,3.720E+2,3.240E+2,2.810E+2,2.430E+2,2.090E+2,1.790E+2, 3 1.530E+2,1.300E+2,1.110E+2,9.500E+1,8.120E+1,6.950E+1,5.950E+1, 4 5.100E+1,4.370E+1,3.760E+1,3.220E+1,2.770E+1,1.320E+1,6.520E+0, 5 3.330E+0,1.760E+0,9.510E-1,6.710E-2,3.000E-4/ DATA PMAC7/ 1 1.013E+3,9.020E+2,8.020E+2,7.100E+2,6.280E+2,5.540E+2,4.870E+2, 2 4.260E+2,3.720E+2,3.240E+2,2.810E+2,2.430E+2,2.090E+2,1.790E+2, 3 1.530E+2,1.300E+2,1.110E+2,9.500E+1,8.120E+1,6.950E+1,5.950E+1, 4 5.100E+1,4.370E+1,3.760E+1,3.220E+1,2.770E+1,1.320E+1,6.520E+0, 5 3.330E+0,1.760E+0,9.510E-1,6.710E-2,3.000E-4/ C THE USER S PRESSURES ARE IN MB FOR THIS PGM. DO 10 I=1,33 II=I IF (P.GT.PMAC(I,L)) GO TO 20 10 CONTINUE 20 CONTINUE IF (II.EQ.1) II=2 W1=(PMAC(II-1,L)-P)/(PMAC(II-1,L)-PMAC(II,L)) W2=(P-PMAC(II,L))/(PMAC(II-1,L)-PMAC(II,L)) QO3=O3MAC(II,L)*W1+O3MAC(II-1,L)*W2 C RO3=QO3(IN GM/M**3)/DENSITY OF AIR(IN GM/M**3);HENCE THE 10**9 R=QO3*8.31432E7*T/(1.0E9*P*28.9644) RETURN END FUNCTION PATH(A,B,C,E) C.... C DOUBLE PRECISION XA,CA COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP PEXP=1./SEXP PATH=((A-B)**PEXP*(A+B+C))/(E*(A+B+C)+(A-B)**(PEXP-1.)) RETURN END SUBROUTINE QINTRP(XM,X0,XP,FM,F0,FP,X,F) C.... C DOUBLE PRECISION FM,F0,FP,F,D1,D2,B,A,DEL D1=(FP-F0)/(XP-X0) D2=(FM-F0)/(XM-X0) B=(D1-D2)/(XP-XM) A=D1-B*(XP-X0) DEL=(X-X0) F=F0+DEL*(A+DEL*B) RETURN END SUBROUTINE QUADSR(NLV,NLP1V,NLP2V,P,PD,TRNS) COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N DIMENSION P(NLP1V),PD(NLP2V),TRNS(NLP1V,NLP1V) DIMENSION WT(101) N2=2*N N2P=2*N+1 C *****WEIGHTS ARE CALCULATED WT(1)=1. DO 21 I=1,N WT(2*I)=4. WT(2*I+1)=1. 21 CONTINUE IF (N.EQ.1) GO TO 25 DO 22 I=2,N WT(2*I-1)=2. 22 CONTINUE 25 CONTINUE TRNSNB=0. DP=(PD(IA)-PD(IA-1))/N2 PFIX=P(JA) DO 1 KK=1,N2P PVARY=PD(IA-1)+(KK-1)*DP IF (PVARY.GE.PFIX) P2=PVARY IF (PVARY.GE.PFIX) P1=PFIX IF (PVARY.LT.PFIX) P1=PVARY IF (PVARY.LT.PFIX) P2=PFIX CALL SINTR2 TRNSNB=TRNSNB+TRNSLO*WT(KK) 1 CONTINUE TRNS(IA,JA)=TRNSNB*DP/(3.*(PD(IA)-PD(IA-1))) RETURN END SUBROUTINE RCTRNS(ITAP1,RATSTD,RATSM,RATIO) C RATSTD=VALUE OF HIGHER STD CO2 CONCENTRATION C RATSM=VALUE OF LOWER STD CO2 CONCENTRATION C RATIO=ACTUAL CO2 CONCENTRATION C THE 3 ABOVE QUANTITIES ARE IN UNITS OF 330 PPMV. COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N COMMON/PRESS/PA(109) COMMON/TRAN/TRANSA(109,109) DIMENSION TRNS1(109,109),TRNS2(109,109) 100 FORMAT (F20.14) C READ IN TFS OF LOWER STD CO2 CONCENTRATION C-CP READ (ITAP1,100) ((TRNS1(I,J),I=1,109),J=1,109) DO J=1,109 DO I=1,109 READ(ITAP1,100) TRANSA(I,J) END DO END DO ITAP2=ITAP1+1 C READ IN TFS OF HIGHER STD CO2 CONCENTRATION C-CP READ (ITAP2,100) ((TRANSA(I,J),I=1,109),J=1,109) DO J=1,109 DO I=1,109 READ(ITAP2,100) TRANSA(I,J) END DO END DO CALL COEINT(RATSTD) DO 401 I=1,109 DO 401 J=1,I IF (J.EQ.I) GO TO 401 C USING HIGHER CO2 CONCENTRATION,COMPUTE 1ST GUESS CO2 TFS FOR C ACTUAL CO2 CONCENTRATION. P2=(RATIO+RATSTD)*PA(I)/(2.*RATSTD) + 1 (RATSTD-RATIO)*PA(J)/(2.*RATSTD) P1=(RATSTD-RATIO)*PA(I)/(2.*RATSTD) + 1 (RATIO+RATSTD)*PA(J)/(2.*RATSTD) CALL SINTR2 TRNSPR=TRNSLO C USING HIGHER CO2 CONCENTRATION,COMPUTE 1ST GUESS CO2 TFS FOR C LOWER STD CO2 CONCENTRATION P2=(RATSM+RATSTD)*PA(I)/(2.*RATSTD) + 1 (RATSTD-RATSM)*PA(J)/(2.*RATSTD) P1=(RATSTD-RATSM)*PA(I)/(2.*RATSTD) + 1 (RATSM+RATSTD)*PA(J)/(2.*RATSTD) CALL SINTR2 TRNSPM=TRNSLO C COMPUTE TFS FOR CO2 CONCENTRATION GIVEN BY (RATIO). C STORE TEMPORARILY IN (TRNS2) TRNS2(J,I)=TRNSPR+(RATSTD-RATIO)*(TRNS1(J,I)- 1 TRNSPM)/(RATSTD-RATSM) TRNS2(I,J)=TRNS2(J,I) C WE NOW CAN OVERWRITE (TRNS1) AND STORE IN (TRNS1) THE 1ST GUESS C CO2 TFS FOR LOWER STD CO2 CONCENTRATION TRNS1(J,I)=TRNSLO TRNS1(I,J)=TRNSLO 401 CONTINUE C SET DIAGONAL VALUES OF CO2 TFS TO UNITY DO 402 I=1,109 TRNS1(I,I)=1.0 TRNS2(I,I)=1.0 402 CONTINUE C NOW OUTPUT THE COMPUTED CO2 TFS FOR (RATIO) CO2 CONC. IN (TRANSA) DO 405 I=1,109 DO 405 J=1,109 TRANSA(J,I)=TRNS2(J,I) 405 CONTINUE RETURN END CCCC SUBROUTINE SIGP(PSTAR,PD,GTEMP) SUBROUTINE SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, 1 LREAD) INCLUDE "parmco2" PARAMETER (KD=LERIC) PARAMETER (KP=KD+1,KM=KD-1,KP2=KD+2) DIMENSION Q(KD),QMH(KP),PD(KP2),PLM(KP),GTEMP(KP),PDT(KP2) DIMENSION SIGLY(KD),SIGLV(KP) DIMENSION CI(KP),SGLVNU(KP),DEL(KD),SIGLNU(KD),CL(KD),RPI(KM) DIMENSION IDATE(4) DIMENSION T41(KP2,2),T42(KP), 1 T43(KP2,2),T44(KP) NAMELIST /PERIC/ PPTOP CCC 18 LEVEL SIGMAS FOR NMC MRF(NEW) MODEL CCC DATA Q/.021,.074,.124,.175,.225,.275,.325,.375,.425,.497, CCC 1 .594,.688,.777,.856,.920,.960,.981,.995/ C FOR SIGMA MODELS,Q=SIGMA,QMH=0.5(Q(I)+Q(I+1), C PD=Q*PSS,PLM=QMH*PSS.PSS=SURFACE PRESSURE(SPEC.) C C..... GET NMC SIGMA STRUCTURE CCC IF (LREAD.GT.0) GO TO 914 C--- PPTOP IS MODEL TOP PRESSURE IN CB.... C SIGMA DATA IS BOTTOM OF ATMOSPHERE TO T.O.A..... PPTOP=5.0 READ(5,PERIC,END=12321) 12321 CONTINUE WRITE(6,88221)PPTOP,KD,KP 88221 FORMAT(' ENTER SIGP PPTOP=',E12.5,' KD=',I2,' KP=',I2) REWIND 23 READ(23)SIGLY WRITE(6,88222) 88222 FORMAT(' READ AETA') DO 37821 LLL=1,KD WRITE(6,37820)LLL,SIGLY(LLL) 37820 FORMAT(' L=',I2,' AETA=',E12.5) 37821 CONTINUE READ(23)SIGLV WRITE(6,88223) 88223 FORMAT(' READ ETA') DO 37823 LLL=1,KP WRITE(6,37822)LLL,SIGLV(LLL) 37822 FORMAT(' L=',I2,' ETA=',E12.5) 37823 CONTINUE 701 FORMAT(F6.2) 702 FORMAT(7F10.6) IF (PPTOP.LE.0.) GO TO 708 PSFC=100. C--- IF PTOP NOT EQUAL TO ZERO ADJUST SIGMA SO AS TO GET PROPER STD ATM C VERTICAL LOCATION DO 706 K=1,KD SIGLY(K) = (SIGLY(K)*(PSFC-PPTOP)+PPTOP)/PSFC 706 CONTINUE DO 707 K=1,KP SIGLV(K) = (SIGLV(K)*(PSFC-PPTOP)+PPTOP)/PSFC 707 CONTINUE 708 CONTINUE PRINT 703,PPTOP PRINT 704,(SIGLY(K),K=1,KD) PRINT 704,(SIGLV(K),K=1,KP) 703 FORMAT(1H ,'PTOP =',F6.2) 704 FORMAT(1H ,7F10.6) DO 913 K=1,KP SGLVNU(K) = SIGLV(K) IF (K.LE.KD) SIGLNU(K) = SIGLY(K) 913 CONTINUE CC914 CALL SETSIG(CI,SGLVNU,DEL,SIGLNU,CL,RPI) DO 77 K=1,KD Q(K) = SIGLNU(KD+1-K) 77 CONTINUE PSS= 1013250. QMH(1)=0. QMH(KP)=1. DO 1 K=2,KD QMH(K)=0.5*(Q(K-1)+Q(K)) 1 CONTINUE PD(1)=0. PD(KP2)=PSS DO 2 K=2,KP PD(K)=Q(K-1)*PSS 2 CONTINUE PLM(1)=0. DO 3 K=1,KM PLM(K+1)=0.5*(PD(K+1)+PD(K+2)) 3 CONTINUE PLM(KP)=PSS DO 4 K=1,KD GTEMP(K)=PD(K+1)**0.2*(1.+PD(K+1)/30000.)**0.8/1013250. 4 CONTINUE GTEMP(KP)=0. C+++ WRITE (6,100) (GTEMP(K),K=1,KD) C+++ WRITE (6,100) (PD(K),K=1,KP2) C+++ WRITE (6,100) (PLM(K),K=1,KP) C***TAPES 41,42 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM (PS=1013MB) C THE FOLLOWING PUTS P-DATA INTO MB DO 11 I=1,KP PD(I)=PD(I)*1.0E-3 PLM(I)=PLM(I)*1.0E-3 11 CONTINUE PD(KP2)=PD(KP2)*1.0E-3 CCC WRITE (41,101) (PD(K),K=1,KP2) CCC WRITE (41,101) (PLM(K),K=1,KP) CCC WRITE (42,101) (PLM(K),K=1,KP) DO 300 K=1,KP2 T41(K,1) = PD(K) 300 CONTINUE DO 301 K=1,KP T41(K,2) = PLM(K) T42(K) = PLM(K) 301 CONTINUE C***STORE AS PDT,SO THAT RIGHT PD IS RETURNED TO PTZ DO 12 I=1,KP2 PDT(I)=PD(I) 12 CONTINUE C***SECOND PASS: PSS=810MB,GTEMP NOT COMPUTED PSS=0.8*1013250. QMH(1)=0. QMH(KP)=1. DO 201 K=2,KD QMH(K)=0.5*(Q(K-1)+Q(K)) 201 CONTINUE PD(1)=0. PD(KP2)=PSS DO 202 K=2,KP PD(K)=Q(K-1)*PSS 202 CONTINUE PLM(1)=0. DO 203 K=1,KM PLM(K+1)=0.5*(PD(K+1)+PD(K+2)) 203 CONTINUE PLM(KP)=PSS C+++ WRITE (6,100) (PD(K),K=1,KP2) C+++ WRITE (6,100) (PLM(K),K=1,KP) C***TAPES 43,44 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM(PS=810 MB) C THE FOLLOWING PUTS P-DATA INTO MB DO 211 I=1,KP PD(I)=PD(I)*1.0E-3 PLM(I)=PLM(I)*1.0E-3 211 CONTINUE PD(KP2)=PD(KP2)*1.0E-3 CCC WRITE (43,101) (PD(K),K=1,KP2) CCC WRITE (43,101) (PLM(K),K=1,KP) CCC WRITE (44,101) (PLM(K),K=1,KP) DO 302 K=1,KP2 T43(K,1) = PD(K) 302 CONTINUE DO 303 K=1,KP T43(K,2) = PLM(K) T44(K) = PLM(K) 303 CONTINUE C***RESTORE PD DO 212 I=1,KP2 PD(I)=PDT(I) 212 CONTINUE 100 FORMAT (1X,5E20.13) 101 FORMAT (5E16.9) RETURN END SUBROUTINE SINTR2 C.... C IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL P1,P2,PA,TRNSLO,CORE,TRANSA,PATH,UEXP,SEXP,ETA,SEXPV COMMON/INPUT/P1,P2,TRNSLO,IA,JA,N COMMON/PRESS/ PA(109) COMMON/TRAN/ TRANSA(109,109) COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP DO 70 L=1,109 IP1=L IF (P2-PA(L)) 65,65,70 70 CONTINUE 65 I=IP1-1 IF (IP1.EQ.1) IP1=2 IF (I.EQ.0) I=1 DO 80 L=1,109 JP1=L IF (P1-PA(L)) 75,75,80 80 CONTINUE 75 J=JP1-1 IF (JP1.EQ.1) JP1=2 IF (J.EQ.0) J=1 JJJ=J III=I J=JJJ JP1=J+1 I=III IP1=I+1 C DETERMINE ETAP,THE VALUE OF ETA TO USE BY LINEAR INTERPOLATION C FOR PETA(=0.5*(P1+P2)) PETA=P2 DO 90 L=1,109 IETAP1=L IF (PETA-PA(L)) 85,85,90 90 CONTINUE 85 IETA=IETAP1-1 IF (IETAP1.EQ.1) IETAP1=2 IF (IETA.EQ.0) IETA=1 ETAP=ETA(IETA)+(PETA-PA(IETA))*(ETA(IETAP1)-ETA(IETA))/ 1 (PA(IETAP1)-PA(IETA)) SEXP=SEXPV(IETA)+(PETA-PA(IETA))*(SEXPV(IETAP1)- 1 SEXPV(IETA))/ (PA(IETAP1)-PA(IETA)) PIPMPI=PA(IP1)-PA(I) UP2P1=(PATH(P2,P1,CORE,ETAP))**UEXP IF (I-J) 126,126,127 126 CONTINUE TRIP=(CA(IP1)*DLOG(1.0D0+XA(IP1)*UP2P1))**(SEXP/UEXP) TRI=(CA(I)*DLOG(1.0D0+XA(I)*UP2P1))**(SEXP/UEXP) TRNSLO=1.0D0-((PA(IP1)-P2)*TRI+(P2-PA(I))*TRIP)/PIPMPI GO TO 128 127 TIJ=TRANSA(I,J) TIPJ=TRANSA(I+1,J) TIJP=TRANSA(I,J+1) TIPJP=TRANSA(I+1,J+1) UIJ=(PATH(PA(I),PA(J),CORE,ETAP))**UEXP UIPJ=(PATH(PA(I+1),PA(J),CORE,ETAP))**UEXP UIJP=(PATH(PA(I),PA(J+1),CORE,ETAP))**UEXP UIPJP=(PATH(PA(I+1),PA(J+1),CORE,ETAP))**UEXP PRODI=CA(I)*XA(I) PRODIP=CA(I+1)*XA(I+1) PROD=((PA(I+1)-P2)*PRODI+(P2-PA(I))*PRODIP)/PIPMPI XINT=((PA(I+1)-P2)*XA(I)+(P2-PA(I))*XA(I+1))/PIPMPI CINT=PROD/XINT AIJ=(CINT*DLOG(1.0D0+XINT*UIJ))**(SEXP/UEXP) AIJP=(CINT*DLOG(1.0D0+XINT*UIJP))**(SEXP/UEXP) AIPJ=(CINT*DLOG(1.0D0+XINT*UIPJ))**(SEXP/UEXP) AIPJP=(CINT*DLOG(1.0D0+XINT*UIPJP))**(SEXP/UEXP) EIJ=TIJ+AIJ EIPJ=TIPJ+AIPJ EIJP=TIJP+AIJP EIPJP=TIPJP+AIPJP DTDJ=(EIJP-EIJ)/(PA(J+1)-PA(J)) DTDPJ=(EIPJP-EIPJ)/(PA(J+1)-PA(J)) EPIP1=EIJ+DTDJ*(P1-PA(J)) EPIPP1=EIPJ+DTDPJ*(P1-PA(J)) EPP2P1=((PA(I+1)-P2)*EPIP1+(P2-PA(I))*EPIPP1)/PIPMPI TRNSLO=EPP2P1-(CINT*DLOG(1.0D0+XINT*UP2P1))**(SEXP/UEXP) IF (I.GE.108.OR.J.GE.108) GO TO 350 IF (I-J-2) 350,350,355 355 CONTINUE TIP2J=TRANSA(I+2,J) TIP2JP=TRANSA(I+2,J+1) TI2J2=TRANSA(I+2,J+2) TIJP2=TRANSA(I,J+2) TIPJP2=TRANSA(I+1,J+2) UIP2J=(PATH(PA(I+2),PA(J),CORE,ETAP))**UEXP UIJP2=(PATH(PA(I),PA(J+2),CORE,ETAP))**UEXP UIPJP2=(PATH(PA(I+1),PA(J+2),CORE,ETAP))**UEXP UI2J2=(PATH(PA(I+2),PA(J+2),CORE,ETAP))**UEXP UIP2JP=(PATH(PA(I+2),PA(J+1),CORE,ETAP))**UEXP AIJP2=(CINT*DLOG(1.0D0+XINT*UIJP2))**(SEXP/UEXP) AIPJP2=(CINT*DLOG(1.0D0+XINT*UIPJP2))**(SEXP/UEXP) AIP2J=(CINT*DLOG(1.0D0+XINT*UIP2J))**(SEXP/UEXP) AIP2JP=(CINT*DLOG(1.0D0+XINT*UIP2JP))**(SEXP/UEXP) AI2J2=(CINT*DLOG(1.0D0+XINT*UI2J2))**(SEXP/UEXP) EIP2J=TIP2J+AIP2J EIP2JP=TIP2JP+AIP2JP EIJP2=TIJP2+AIJP2 EIPJP2=TIPJP2+AIPJP2 EI2J2=TI2J2+AI2J2 CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIJ,EIJP,EIJP2,P1,EI) CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIPJ,EIPJP,EIPJP2,P1,EP) CALL QINTRP(PA(J),PA(J+1),PA(J+2),EIP2J,EIP2JP,EI2J2,P1,EP2) CALL QINTRP(PA(I),PA(I+1),PA(I+2),EI,EP,EP2,P2,EPSIL) TRNSLO=EPSIL-(CINT*DLOG(1.0D0+XINT*UP2P1))**(SEXP/UEXP) 350 CONTINUE 128 CONTINUE 205 CONTINUE RETURN END CCCC PROGRAM CO2O3 = CONSOLIDATION OF A NUMBER OF DAN SCHWARZKOPF,GFDL C CODES TO PRODUCE A FILE OF CO2 HGT DATA C FOR ANY VERTICAL COORDINATE (READ BY SUBROUTINE C CONRAD IN THE GFDL RADIATION CODES)-K.A.C. JUN89. CNOV89--UPDATED (NOV 89) FOR LATEST GFDL LW RADIATION.....K.A.C. INCLUDE "parmco2" PARAMETER (L=LERIC,LP1=L+1,LP2=L+2) DIMENSION SGTEMP(LP1,2),CO2D1D(L,6),CO2D2D(LP1,LP1,6) CNOV89 DIMENSION CO2IQ2(LP1,LP1,6),CO2IQ3(LP1,LP1,6),CO2IQ5(LP1,LP1,6) CNOV89 DIMENSION T41(LP2,2),T42(LP1), 1 T43(LP2,2),T44(LP1) DIMENSION T20(LP1,LP1,3),T21(LP1,LP1,3) DIMENSION T22(LP1,LP1,3),T23(LP1,LP1,3) DIMENSION SGLVNU(LP1),SIGLNU(L) DIMENSION STEMP(LP1),GTEMP(LP1) DIMENSION CDTM51(L),CO2M51(L),C2DM51(L) DIMENSION CDTM58(L),CO2M58(L),C2DM58(L) DIMENSION CDT51(LP1,LP1),CO251(LP1,LP1),C2D51(LP1,LP1) DIMENSION CDT58(LP1,LP1),CO258(LP1,LP1),C2D58(LP1,LP1) CNOV89 DIMENSION CDT31(LP1),CO231(LP1),C2D31(LP1) DIMENSION CDT38(LP1),CO238(LP1),C2D38(LP1) DIMENSION CDT71(LP1),CO271(LP1),C2D71(LP1) DIMENSION CDT78(LP1),CO278(LP1),C2D78(LP1) DIMENSION CO211(LP1),CO218(LP1) EQUIVALENCE (CDT31(1),CO2IQ2(1,1,1)),(CO231(1),CO2IQ2(1,1,2)) EQUIVALENCE (C2D31(1),CO2IQ2(1,1,3)),(CDT38(1),CO2IQ2(1,1,4)) EQUIVALENCE (CO238(1),CO2IQ2(1,1,5)),(C2D38(1),CO2IQ2(1,1,6)) EQUIVALENCE (CDT71(1),CO2IQ3(1,1,1)),(CO271(1),CO2IQ3(1,1,2)) EQUIVALENCE (C2D71(1),CO2IQ3(1,1,3)),(CDT78(1),CO2IQ3(1,1,4)) EQUIVALENCE (CO278(1),CO2IQ3(1,1,5)),(C2D78(1),CO2IQ3(1,1,6)) EQUIVALENCE (CO211(1),CO2IQ5(1,1,2)),(CO218(1),CO2IQ5(1,1,5)) CNOV89 EQUIVALENCE (STEMP(1),SGTEMP(1,1)),(GTEMP(1),SGTEMP(1,2)) EQUIVALENCE (CDTM51(1),CO2D1D(1,1)),(CO2M51(1),CO2D1D(1,2)) EQUIVALENCE (C2DM51(1),CO2D1D(1,3)),(CDTM58(1),CO2D1D(1,4)) EQUIVALENCE (CO2M58(1),CO2D1D(1,5)),(C2DM58(1),CO2D1D(1,6)) EQUIVALENCE (CDT51(1,1),CO2D2D(1,1,1)),(CO251(1,1),CO2D2D(1,1,2)) EQUIVALENCE (C2D51(1,1),CO2D2D(1,1,3)),(CDT58(1,1),CO2D2D(1,1,4)) EQUIVALENCE (CO258(1,1),CO2D2D(1,1,5)),(C2D58(1,1),CO2D2D(1,1,6)) C===> GET SGTEMP AND OUTPUT WHICH USED TO BE ON UNITS 41,42,43,44.... LREAD = 0 CALL CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU,LREAD) C===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. C IR=1,IQ=1 IS FOR COMMON /CO2BD3/ IN RADIATION CODE... C FOR THE CONSOLIDATED 490-850 CM-1 BAND... CNOV89 ICO2TP=61 CNOV89 IR = 1 RATIO = 1.0 NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD) IR = 1 RATIO = 1.0 NMETHD = 1 CALL CO2INT(ICO2TP,T41,T42,T20,RATIO,IR,NMETHD) IR = 1 RATIO = 1.0 NMETHD = 2 CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD) IR = 1 RATIO = 1.0 NMETHD = 1 CALL CO2INT(ICO2TP,T43,T44,T21,RATIO,IR,NMETHD) C===> FILL UP THE CO2D1D ARRAY C THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND C THEIR DERIVATIVES FOR TAU(I,I+1),I=1,LEVS, C WHERE THE VALUES ARE NOT OBTAINED BY QUADRATURE BUT ARE THE C ACTUAL TRANSMISSIVITIES,ETC,BETWEEN A PAIR OF PRESSURES. THESE C ARE USED ONLY FOR NEARBY LAYER CALCULATIONS INCLUDING H2O.. C IQ = 1 CALL CO2IN1(T20,T21,CO2D1D,IQ) C C===> FILL UP THE CO2D2D ARRAY C THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND THEIR DERIVATIVES C FROM 109-LEVEL LINE-BY-LINE CALCULATIONS MADE USING THE 1982 C MCCLATCHY TAPE (12511 LINES),CONSOLIDATED,INTERPOLATED C TO THE MRF VERTICAL COORDINATE,AND RE-CONSOLIDATED TO A C 200 CM-1 BANDWIDTH. THE INTERPOLATION METHOD IS DESCRIBED IN C SCHWARZKOPF AND FELS (J.G.R.,1985). C CALL CO2INS(T22,T23,CO2D2D,IQ) C CNOV89 C===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. C IR=2,IQ=2 IS FOR COMMON /CO2BD2/ IN RADIATION CODE... C FOR THE CONSOLIDATED 490-670 CM-1 BAND... ICO2TP=62 IR = 2 RATIO = 1.0 NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD) CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD) IQ = 2 CALL CO2INS(T22,T23,CO2IQ2,IQ) C===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. C IR=3,IQ=3 IS FOR COMMON /CO2BD4/ IN RADIATION CODE... C FOR THE CONSOLIDATED 670-850 CM-1 BAND... ICO2TP=63 IR = 3 RATIO = 1.0 NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD) CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD) IQ = 3 CALL CO2INS(T22,T23,CO2IQ3,IQ) C--- FOLLOWING CODE NOT WORKING AND NOT NEEDED YET C===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. C IR=4,IQ=5 IS FOR COMMON /CO2BD5/ IN RADIATION CODE... C FOR THE 4.3 MICRON BAND... C NOT USED YET ICO2TP=65 C NOT USED YET IR = 4 C NOT USED YET RATIO = 1.0 C DAN SCHWARZ --- USE 300PPMV RATIO = 0.9091 (NOT TESTED YET)..... C NOT USED YET NMETHD = 2 C NOT USED YET CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD) C NOT USED YET CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD) C NOT USED YET IQ = 5 C NOT USED YET CALL CO2INS(T22,T23,CO2IQ5,IQ) CNOV89 C... WRITE DATA TO DISK.. C ...SINCE THESE CODES ARE COMPILED WITH AUTODBL,THE CO2 DATA C IS CONVERTED TO SINGLE PRECISION IN A LATER JOB STEP.. REWIND 66 WRITE(66) STEMP WRITE(66) GTEMP WRITE(66) CDTM51 WRITE(66) CO2M51 WRITE(66) C2DM51 WRITE(66) CDTM58 WRITE(66) CO2M58 WRITE(66) C2DM58 WRITE(66) CDT51 WRITE(66) CO251 WRITE(66) C2D51 WRITE(66) CDT58 WRITE(66) CO258 WRITE(66) C2D58 CNOV89 WRITE(66) CDT31 WRITE(66) CO231 WRITE(66) C2D31 WRITE(66) CDT38 WRITE(66) CO238 WRITE(66) C2D38 WRITE(66) CDT71 WRITE(66) CO271 WRITE(66) C2D71 WRITE(66) CDT78 WRITE(66) CO278 WRITE(66) C2D78 C NOT USED YET WRITE(66) CO211 C NOT USED YET WRITE(66) CO218 CNOV89 REWIND 66 STOP END SUBROUTINE SETSIG(CI, SI, DEL, SL, CL, RPI) INCLUDE "parmco2" PARAMETER (LEVS=LERIC,LEVP1=LEVS+1,LEVM1=LEVS-1) CCC REAL RK,RK1,RKR DIMENSION CI(LEVP1), SI(LEVP1), 1 DEL(LEVS), SL(LEVS), CL(LEVS), RPI(LEVM1) DIMENSION DELSIG(LEVS) C DATA DELSIG/ .01,.017,0.025,0.055,0.073,0.085,0.093, * 0.096,0.096,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05/ C PRINT 98 98 FORMAT (1H0, 'BEGIN SETSIG') DO 9889 K=1,LEVS 9889 DEL(K)=DELSIG(K) CI(1) = 0.E0 SUMDEL=0.E0 DO 54 K=1,LEVS SUMDEL=SUMDEL+DEL(K) CI(K+1)=CI(K)+DEL(K) 54 CONTINUE CI(LEVP1)=1.E0 RK = 287.05/1004.64 RK1 = RK + 1.E0 RKR = 1.0/RK DO 1 LI=1,LEVP1 1 SI(LI) = 1.E0 - CI(LI) DO 3 LE=1,LEVS DIF = SI(LE)**RK1 - SI(LE+1)**RK1 DIF = DIF / (RK1*(SI(LE)-SI(LE+1))) SL(LE) = DIF**RKR CL(LE) = 1.E0 - SL(LE) 3 CONTINUE C COMPUTE PI RATIOS FOR TEMP. MATRIX. DO 4 LE=1,LEVM1 BASE = SL(LE+1)/SL(LE) 4 RPI(LE) = BASE**RK DO 5 LE=1,LEVP1 PRINT 100, LE, CI(LE), SI(LE) 100 FORMAT (1H , 'LEVEL=', I2, 2X, 'CI=', F6.3, 2X, 'SI=', F6.3) 5 CONTINUE PRINT 97 97 FORMAT (1H0) DO 6 LE=1,LEVS PRINT 101, LE, CL(LE), SL(LE), DEL(LE) 101 FORMAT (1H , 'LAYER=', I2, 2X, 'CL=', F6.3, 2X, 'SL=', F6.3, 2X, 1 'DEL=', F6.3) 6 CONTINUE PRINT 102, (RPI(LE), LE=1,LEVM1) 102 FORMAT (1H0, 'RPI=', (18(1X,F6.3)) ) PRINT 99,LEVS,SUMDEL 99 FORMAT (1H0,'SHALOM FROM MRF',I4,' LEVEL SETSIG--SUMDEL=',E13.4) RETURN END