FUNCTION ANTEMP(L,Z) REAL :: ZB(10,7),C(11,7),DELTA(10,7),TSTAR(7) ! ************** TROPICAL SOUNDING ************************** DATA (ZB(N,1),N=1,10)/ 2.0, 3.0, 16.5, 21.5, 45.0, & 51.0, 70.0, 100., 200., 300./ DATA (C(N,1),N=1,11)/ -6.0, -4.0, -6.7, 4.0, 2.2, & 1.0, -2.8, -.27, 0.0, 0.0, 0.0/ DATA (DELTA(N,1),N=1,10)/.5, .5, .3, .5, 1.0, & 1.0, 1.0, 1.0, 1.0, 1.0/ ! ************** SUB-TROPICAL SUMMER ************************ DATA (ZB(N,2),N=1,10)/ 1.5, 6.5, 13.0, 18.0, 26.0, & 36.0, 48.0, 50.0, 70.0, 100./ DATA (C(N,2),N=1,11)/ -4.0, -6.0, -6.5, 0.0, 1.2, & 2.2, 2.5, 0.0, -3.0, -0.25, 0.0/ DATA (DELTA(N,2),N=1,10)/ .5, 1.0, .5, .5, 1.0, & 1.0, 2.5, .5, 1.0, 1.0/ ! ************** SUB-TROPICAL WINTER ************************ DATA (ZB(N,3),N=1,10)/ 3.0, 10.0, 19.0, 25.0, 32.0, & 44.5, 50.0, 71.0, 98.0, 200.0/ DATA (C(N,3),N=1,11)/ -3.5, -6.0, -0.5, 0.0, 0.4, & 3.2, 1.6, -1.8, -0.7, 0.0, 0.0/ DATA (DELTA(N,3),N=1,10)/ .5, .5, 1.0, 1.0, 1.0, & 1.0, 1.0, 1.0, 1.0, 1.0/ ! ************* SUB-ARCTIC SUMMER ************************* DATA (ZB(N,4),N=1,10)/ 4.7, 10.0, 23.0, 31.8, 44.0, & 50.2, 69.2, 100.0, 102.0, 103.0/ DATA (C(N,4),N=1,11)/ -5.3, -7.0, 0.0, 1.4, 3.0, & 0.7, -3.3, -0.2, 0.0, 0.0, 0.0/ DATA (DELTA(N,4),N=1,10)/ .5, .3, 1.0, 1.0, 2.0, & 1.0, 1.5, 1.0, 1.0, 1.0/ ! ************ SUB-ARCTIC WINTER ***************************** DATA (ZB(N,5),N=1,10)/ 1.0, 3.2, 8.5, 15.5, 25.0, & 30.0, 35.0, 50.0, 70.0, 100.0/ DATA (C(N,5),N=1,11)/ 3.0, -3.2, -6.8, 0.0, -0.6, & 1.0, 1.2, 2.5, -0.7, -1.2, 0.0/ DATA (DELTA(N,5),N=1,10)/ .4, 1.5, .3 , .5, 1.0, & 1.0, 1.0, 1.0, 1.0, 1.0/ ! ************ US STANDARD 1976 ****************************** DATA (ZB(N,6),N=1,10)/ 11.0, 20.0, 32.0, 47.0, 51.0, & 71.0, 84.8520, 90.0, 91.0, 92.0/ DATA (C(N,6),N=1,11)/ -6.5, 0.0, 1.0, 2.80, 0.0, & -2.80, -2.00, 0.0, 0.0, 0.0, 0.0/ DATA (DELTA(N,6),N=1,10)/ 0.3, 1.0, 1.0, 1.0, 1.0, & 1.0, 1.0, 1.0, 1.0, 1.0/ ! ! ************ ENLARGED US STANDARD 1976 ********************** DATA (ZB(N,7),N=1,10)/ 11.0, 20.0, 32.0, 47.0, 51.0, & 71.0, 84.8520, 90.0, 91.0, 92.0/ DATA (C(N,7),N=1,11)/ -6.5, 0.0, 1.0, 2.80, 0.0, & -2.80, -2.00, 0.0, 0.0, 0.0, 0.0/ DATA (DELTA(N,7),N=1,10)/ 0.3, 1.0, 1.0, 1.0, 1.0, & 1.0, 1.0, 1.0, 1.0, 1.0/ ! DATA TSTAR/ 300.0, 294.0, 272.2, 287.0, 257.1, 2*288.15/ ! 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) !JD single-precision change ! FAC=EXP(EXPP)+EXP(-EXPP) !mp write(6,*) '.........................................' !mp what in the hell does the next line do? !mp !mp apparently if statement <0 or =0 then 23, else 24 !mp IF(ABS(EXPO)-100.0) 23,23,24 ! ! changed to a more reasonable value for the workstation ! IF(ABS(EXPO)-50.0) 23,23,24 23 X=EXP(EXPO) Y=X+1.0/X ZLOG=ALOG(Y) GO TO 25 24 ZLOG=ABS(EXPO) !mp 25 IF(EXPP-100.0) 27,27,28 25 IF(EXPP-50.0) 27,27,28 !JD single-precision change 27 FAC=EXP(EXPP)+EXP(-EXPP) FACLOG=ALOG(FAC) GO TO 29 28 FACLOG=EXPP ! TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)* ! 1 ALOG((EXP(EXPO)+EXP(-EXPO))/FAC)) 29 TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)* & (ZLOG-FACLOG)) !mp write(6,*) 'ANTEMP pieces (C,C,ZLOG,FACLOG)', C(N+1,L),C(N,L), !mp + ZLOG,FACLOG 20 CONTINUE ANTEMP=TEMP END FUNCTION ANTEMP !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc SUBROUTINE COEINT(RAT,IR) ! ********************************************************************** ! ! ! THE TRANSMISSION FUNCTION BETWEEN P1 AND P2 IS ASSUMED TO ! THE FUNCTIONAL FORM ! TAU(P1,P2)= 1.0-SQRT(C*LOG(1.0+X*PATH)), ! WHERE ! PATH(P1,P2)=((P1-P2)**2)*(P1+P2+CORE)/ ! (ETA*(P1+P2+CORE)+(P1-P2)) ! ! ! THE PARAMETERS C AND X ARE FUNCTIONS OF P2, AND ARE TO BE DETER ! WHILE CORE IS A PRESPECIFIED NUMBER.ETA IS A FUNCTION OF THE TH ! PRODUCT (CX);IT IS OBTAITED ITERATIVELY. THE DERIVATION OF ALL ! VALUES WILL BE EXPLAINED IN A FORTHCOMING PAPER. ! SUBROUTINE COEINT DETERMINES C(I) AND X(I) BY USING THE ACT ! VALUES OF TAU(P(I-2),P(I)) AND TAU(P(I-1),P(I)) AND THE PREVIOU ! ITERATION VALUE OF ETA. ! DEFINE: ! PATHA=PATH(P(I),P(I-2),CORE,ETA) ! PATHB=PATH(P(I),P(I-1),CORE,ETA); ! THEN ! R=(1-TAU(P(I),P(I-2)))/(1-TAU(P(I),P(I-1))) ! = SQRT(LOG(1+X*PATHA)/LOG(1+X*PATHB)), ! SO THAT ! R**2= LOG(1+X*PATHA)/LOG(1+X*PATHB). ! THIS EQUATION CAN BE SOLVED BY NEWTON S METHOD FOR X AND THEN T ! RESULT USED TO FIND C. THIS IS REPEATED FOR EACH VALUE OF I GRE ! THAN 2 TO GIVE THE ARRAYS X(I) AND C(I). ! NEWTON S METHOD FOR SOLVING THE EQUATION ! F(X)=0 ! MAKES USE OF THE LOOP XNEW= XOLD-F(XOLD)/FPRIME(XOLD). ! THIS IS ITERATED 20 TIMES, WHICH IS PROBABLY EXCESSIVE. ! THE FIRST GUESS FOR ETA IS 3.2E-4*EXP(-P(I)/1000),WHICH HAS ! BEEN FOUND TO BE FAIRLY REALISTIC BY EXPERIMENT; WE ITERATE 5 T ! (AGAIN,PROBABLY EXCESSIVELY) TO OBTAIN THE VALUES FOR C,X,ETA T ! USED FOR INTERPOLATION. ! THERE ARE SEVERAL POSSIBLE PITFALLS: ! 1) IN THE COURSE OF ITERATION, X MAY REACH A VALUE WHICH ! 1+X*PATHA NEGATIVE; IN THIS CASE THE ITERATION IS STOP ! AND AN ERROR MESSAGE IS PRINTED OUT. ! 2) EVEN IF (1) DOES NOT OCCUR, IT IS STILL POSSIBLE THAT ! BE NEGATIVE AND LARGE ENOUGH TO MAKE 1+X*PATH(P(I),0,C ! NEGATIVE. THIS IS CHECKED FOR IN A FINAL LOOP, AND IF ! A WARNING IS PRINTED OUT. ! ! ********************************************************************* !.... ! 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) INTEGER :: IERR DATA SINV/2.74992,2.12731,4.38111,0.0832926/ !NOV89 DIMENSION SINV(3) !NOV89 DATA SINV/2.74992,2.12731,4.38111/ !O222 OLD CODE USED 2.7528 RATHER THAN 2.74992 ---K.A.C. OCTOBER 1988 !O222 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) arg1=path(pa(i),pa(i-2),core,eta(i)) arg2=path(pa(i),pa(i-1),core,eta(i)) 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))/ & (F2*F2) XX=XX-F/FPRIME CHECK=1.0D0+XX*PATHA !!!! IF (CHECK) 1020,1020,1025 IF(CHECK.LE.0.)THEN WRITE(errmess,360)I,LL,CHECK WRITE(errmess,*)' xx=',xx,' patha=',patha 360 FORMAT(' ERROR,I=',I3,'LL=',I3,'CHECK=',F20.10) CALL wrf_error_fatal ( errmess ) ENDIF 1010 CONTINUE CA(I)=(1.0D0-TRANSA(I,I-2))**(UEXP/SEXP)/ & (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) !+++ 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)* & (CA(I)*XA(I))**(1./UEXP) 1035 CONTINUE ! ! THE ETA FORMULATION IS DETAILED IN SCHWARZKOPF AND FELS(1985). ! THE QUANTITY SINV=(G*DELTANU)/(RCO2*D*S) ! IN CGS UNITS,WITH D,THE DIFFUSICITY FACTOR=2, AND ! S,THE SUM OF CO2 LINE STRENGTHS OVER THE 15UM CO2 BAND ! ALSO,THE DENOMINATOR IS MULTIPLIED BY ! 1000 TO PERMIT USE OF MB UNITS FOR PRESSURE. ! S IS ACTUALLY WEIGHTED BY B(250) AT 10 CM-1 WIDE INTERVALS,IN ! ORDER TO BE CONSISTENT WITH THE METHODS USED TO OBTAIN THE LBL ! 1-BAND CONSOLIDATED TRANCMISSION FUNCTIONS. ! FOR THE 490-850 INTERVAL (DELTANU=360,IR=1) SINV=2.74992. ! (SLIGHTLY DIFFERENT FROM 2.7528 USED IN EARLIER VERSIONS) ! FOR THE 490-670 INTERVAL (IR=2) SINV=2.12731 ! FOR THE 670-850 INTERVAL (IR=3) SINV=4.38111 ! FOR THE 2270-2380 INTERVAL (IR=4) SINV=0.0832926 ! SINV HAS BEEN OBTAINED USING THE 1982 AFGL CATALOG FOR CO2 ! RAT IS THE ACTUAL CO2 MIXING RATIO IN UNITS OF 330 PPMV, ! LETTING USE OF THIS FORMULATION FOR ANY CO2 CONCENTRATION. ! ! WRITE (6,366) (NP,I,CA(I),XA(I),ETA(I),SEXPV(I),I=1,109) !366 FORMAT (2I4,4E20.12) 1200 CONTINUE 361 FORMAT (' **WARNING:** 1+XA*PATH(PA(I),0) IS NEGATIVE,I= ',I3,/ & 20X,'PATH0(I)=',F16.6,' XA(I)=',F16.6) RETURN END SUBROUTINE COEINT !-------------- !CCC PROGRAM CO2INS SUBROUTINE CO2INS(T22,T23,T66,IQ,L,LP1,iflag) ! ********************************************************* ! SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 ****** ! ..... K.CAMPANA MARCH 1988,OCTOBER 1988... ! ..... K.CAMPANA DECEMBER 1988-CLEANED UP FOR LAUNCHER ! ..... K.CAMPANA NOVEMBER 1989-ALTERED FOR NEW RADIATION ! ********************************************************* 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) !CC ITIN=22 !CC ITIN1=23 !O222 LATEST CODE HAD IQ=1 !CC IQ=4 1011 FORMAT (4F20.14) !CC READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1) !CC 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) !NOV89 IF (IQ.EQ.5) GO TO 300 !NOV89 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) !NOV89 IF (IQ.EQ.5) GO TO 301 !NOV89 CO2801(I,J) = T23(I,J,2) CO2802(I,J) = T23(I,J,3) 301 CONTINUE !***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS ! ARE: ! IQ=1 560-800 (CONSOL.=490-850) ! IQ=2 560-670 (CONSOL.=490-670) ! IQ=3 670-800 (CONSOL.=670-850) ! IQ=4 560-760 (ORIGINAL CODE) (CONSOL.=490-850) !NOV89 ! IQ=5 2270-2380 (CONSOL.=2270-2380) !NOV89 ! THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS ! USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT ! WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S. !NOV89 ! NOTE: ALTHOUGH THE BAND TRANSMISSION FUNCTIONS ARE ! COMPUTED FOR ALL RADIATIVE BANDS, AS OF 9/28/88, THEY ! ARE WRITTEN OUT IN FULL ONLY FOR THE FULL 15 UM BAND CASES ! (IQ=1,4). IN OTHER CASES, THE TRANSMISSIVITIES (1,K) ARE ! WRITTEN OUT, AS THESE ARE THE ONLY ONES NEEDED FOR CTS ! CALCULATIONS. ALSO, FOR THE 4.3 UM BAND (IQ=5) THE TEMP. ! DERIVATIVE TERMS ARE NOT WRITTEN OUT, AS THEY ARE UNUSED. !NOV89 IF (IQ.EQ.1) THEN C1=1.5 C2x=0.5 ENDIF IF (IQ.EQ.2) THEN C1=18./11. C2x=7./11. ENDIF IF (IQ.EQ.3) THEN C1=18./13. C2x=5./13. ENDIF IF (IQ.EQ.4) THEN C1=1.8 C2x=0.8 ENDIF !NOV89 IF (IQ.EQ.5) THEN C1=1.0 C2x=0.0 ENDIF !NOV89 DO 1021 I=1,LP1 DO 1021 J=1,LP1 CO2PO(J,I)=C1*CO2PO(J,I)-C2x CO2800(J,I)=C1*CO2800(J,I)-C2x !NOV89 IF (IQ.EQ.5) GO TO 1021 !NOV89 CO2PO1(J,I)=C1*CO2PO1(J,I)-C2x CO2801(J,I)=C1*CO2801(J,I)-C2x CO2PO2(J,I)=C1*CO2PO2(J,I)-C2x CO2802(J,I)=C1*CO2802(J,I)-C2x 1021 CONTINUE !NOV89 IF (IQ.GE.1.AND.IQ.LE.4) THEN !NOV89 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 !NOV89 ENDIF !NOV89 !O222 ********************************************************* !CC REWIND 66 ! SAVE CDT51,CO251,C2D51,CDT58,CO258,C2D58..ON TEMPO FILE !CC WRITE (66) DCDT10 !CC WRITE (66) CO2PO !CC WRITE (66) D2CT10 !CC WRITE (66) DCDT8 !CC WRITE (66) CO2800 !CC WRITE (66) D2CT8 !CC REWIND 66 !NOV89 IF (IQ.EQ.1.OR.IQ.EQ.4) THEN !NOV89 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 !NOV89 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 !NOV89 !O222 ********************************************************* RETURN END SUBROUTINE CO2INS !O222 PROGRAM CO2INT(INPUT,TAPE5=INPUT) !NOV89 SUBROUTINE CO2INT(ITAPE,T15A,T15B,T22,RATIO,IR,NMETHD,NLEVLS,NLP1,NLP2) !NOV89 ! ********************************************************* ! CHANGES TO DATA READ AND FORMAT SEE CO222 *** ! ..... K.CAMPANA MARCH 1988,OCTOBER 1988 ! CHANGES TO PASS ITAPE,AND IF IR=4,READ 1 CO2 REC..KAC NOV89 ! ********************************************************* ! CO2INT INTERPOLATES CARBON DIOXIDE TRANSMISSION FUNCTIONS ! FROM THE 109 LEVEL GRID,FOR WHICH THE TRANSMISSION FUNCTIONS ! HAVE BEEN PRE-CALCULATED, TO THE GRID STRUCTURE SPECIFIED BY THE ! USER. ! ! METHOD: ! ! CO2INT IS EMPLOYABLE FOR TWO PURPOSES: 1) TO OBTAIN TRANSMIS- ! SIVITIES BETWEEN ANY 2 OF AN ARRAY OF USER-DEFINED PRESSURES; AND ! 2) TO OBTAIN LAYER-MEAN TRANSMISSIVITIES BETWEEN ANY 2 OF AN ARRAY ! OF USER-DEFINED PRESSURE LAYERS.TO CLARIFY THESE TWO PURPOSES,SEE ! THE DIAGRAM AND DISCUSSION BELOW. ! CO2INT MAY BE USED TO EXECUTE ONLY ONE PURPOSE AT ONE TIME. ! ! LET P BE AN ARRAY OF USER-DEFINED PRESSURES ! AND PD BE USER-DEFINED PRESSURE LAYERS. ! ! - - - - - - - - - PD(I-1) --- ! ^ ! ----------------- P(I) ^ PRESSURE LAYER I (PLM(I)) ! ^ ! - - - - - - - - - PD(I) --- ! ^ ! ----------------- P(I+1) ^ PRESSURE LAYER I+1 (PLM(I+1)) ! ^ ! - - - - - - - - - PD(I+1)--- ! ... (THE NOTATION USED IS ! ... CONSISTENT WITH THE CODE) ! ... ! - - - - - - - - - PD(J-1) ! ! ----------------- P(J) ! ! - - - - - - - - - PD(J) ! ! PURPOSE 1: THE TRANSMISSIVITY BETWEEN SPECIFIC PRESSURES ! P(I) AND P(J) ,TAU(P(I),P(J)) IS COMPUTED BY THIS PROGRAM. ! IN THIS MODE,THERE IS NO REFERENCE TO LAYER PRESSURES PD ! (PD,PLM ARE NOT INPUTTED). ! ! PURPOSE 2: THE LAYER-MEAN TRANSMISSIVITY BETWEEN A LAYER- ! MEAN PRESSURE PLM(J) AND PRESSURE LAYER I IS GIVEN BY ! TAULM(PLM(I),PLM(J)). IT IS COMPUTED BY THE INTEGRAL ! ! PD(I) ! ---- ! 1 ^ ! ------------- * ^ TAU ( P',PLM(J) ) DP' ! PD(I)-PD(I-1) ^ ! ---- ! PD(I-1) ! ! THE LAYER-MEAN PRESSURE PLM(I) IS SPECIFIED BY THE USER. ! FOR MANY PURPOSES,PLM WILL BE CHOSEN TO BE THE AVERAGE ! PRESSURE IN THE LAYER-IE,PLM(I)=0.5*(PD(I-1)+PD(I)). ! FOR LAYER-MEAN TRANSMISSIVITIES,THE USER THUS INPUTS ! A PRESSURE ARRAY (PD) DEFINING THE PRESSURE LAYERS AND AN ! ARRAY (PLM) DEFINING THE LAYER-MEAN PRESSURES.THE CALCULATION ! DOES NOT DEPEND ON THE P ARRAY USED FOR PURPOSE 1 (P IS NOT ! INPUTTED). ! ! THE FOLLOWING PARAGRAPHS DEPICT THE UTILIZATION OF THIS ! CODE WHEN USED TO COMPUTE TRANSMISSIVITIES BETWEEN SPECIFIC ! PRESSURES. LATER PARAGRAPHS DESCRIBE ADDITIONAL FEATURES NEEDED ! FOR LAYER-MEAN TRANSMISSIVITIES. ! ! FOR A GIVEN CO2 MIXING RATIO AND STANDARD TEMPERATURE ! PROFILE,A TABLE OF TRANSMISSION FUNCTIONS FOR A FIXED GRID ! OF ATMOSPHERIC PRESSURES HAS BEEN PRE-CALCULATED. ! THE STANDARD TEMPERATURE PROFILE IS COMPUTED FROM THE US ! STANDARD ATMOSPHERE (1977) TABLE.ADDITIONALLY, THE ! SAME TRANSMISSION FUNCTIONS HAVE BEEN PRE-CALCULATED FOR A ! TEMPERATURE PROFILE INCREASED AND DECREASED (AT ALL LEVELS) ! BY 25 DEGREES. ! THIS PROGRAM READS IN THE PRESPECIFIED TRANSMISSION FUNCTIONS ! AND A USER-SUPPLIED PRESSURE GRID (P(I)) AND CALCULATES TRANS- ! MISSION FUNCTIONS ,TAU(P(I),P(J)), FOR ALL P(I) S AND P(J) S. ! A LOGARITHMIC INTERPOLATION SCHEME IS USED. ! THIS METHOD IS REPEATED FOR THE THREE TEMPERATURE PROFILES ! GIVEN ABOVE .THEREFORE OUTPUTS FROM THE PROGRAM ARE THREE TABLES ! OF TRANSMISSION FUNCTIONS FOR THE USER-SUPPLIED PRESSURE GRID. ! THE EXISTENCE OF THE THREE TABLES PERMITS SUBSEQUENT INTERPO- ! LATION TO A USER-SUPPLIED TEMPERATURE PROFILE USING THE METHOD ! DESCRIBED IN THE REFERENCE.SEE LIMITATIONS SECTION IF THE ! USER DESIRES TO OBTAIN ONLY 1 TABLE OF TRANSMISSIVITIES. ! ! MODIFICATIONS FOR LAYER-MEAN TRANSMISSIVITIES: ! THE PRESSURES INPUTTED ARE THE LAYER-MEAN PRESSURES,PD, ! AND THE LAYER-MEAN PRESSURES ,PLM. A SERIES OF TRANSMISSIVITIES ! (TAU(P'',PLM(J)) ARE COMPUTED AND THE INTEGRAL GIVEN IN THE ! DISCUSSION OF PURPOSE 2 IS COMPUTED.FOR PLM(I) NOT EQUAL TO ! PLM(J) SIMPSON S RULE IS USED WITH 5 POINTS. IF PLM(I)=PLM(J) ! (THE -NEARBY LAYER- CASE) A 49-POINT QUADRATURE IS USED FOR ! GREATER ACCURACY.THE OUTPUT IS IN TAULM(PLM(I),PLM(J)). ! NOTE: ! TAULM IS NOT A SYMMETRICAL MATRIX. FOR THE ARRAY ELEMENT ! TAULM(PLM(I),PLM(J)),THE INNER(FIRST,MOST RAPIDLY VARYING) ! DIMENSION IS THE VARYING LAYER-MEAN PRESSURE,PLM(I);THE OUTER ! (SECOND) DIMENSION IS THE FIXED LAYER-MEAN PRESSURE PLM(J). ! THUS THE ELEMENT TAULM(2,3) IS THE TRANSMISSION FUNCTION BETWEEN ! THE FIXED PRESSURE PLM(3) AND THE PRESSURE LAYER HAVING AN AVERAG ! PRESSURE OF PLM(2). ! ALSO NOTE THAT NO QUADRATURE IS PERFORMED OVER THE LAYER ! BETWEEN THE SMALLEST NONZERO PRESSURE AND ZERO PRESSURE; ! TAULM IS TAULM(0,PLM(J)) IN THIS CASE,AND TAULM(0,0)=1. ! ! ! REFERENCE: ! S.B.FELS AND M.D.SCHWARZKOPF,-AN EFFICIENT ACCURATE ! ALGORITHM FOR CALCULATING CO2 15 UM BAND COOLING RATES-,JOURNAL ! OF GEOPHYSICAL RESEARCH,VOL.86,NO. C2, PP.1205-1232,1981. ! MODIFICATIONS TO THE ALGORITHM HAVE BEEN MADE BY THE AUTHORS; ! CONTACT S.B.F.OR M.D.S. FOR FURTHER DETAILS.A NOTE TO J.G.R. ! IS PLANNED TO DOCUMENT THESE CHANGES. ! ! AUTHOR: M.DANIEL SCHWARZKOPF ! ! DATE: 14 JULY 1983 ! ! ADDRESS: ! ! G.F.D.L. ! P.O.BOX 308 ! PRINCETON,N.J.08540 ! U.S.A. ! TELEPHONE: (609) 452-6521 ! ! INFORMATION ON TAPE: THIS SOURCE IS THE FIRST FILE ! ON THIS TAPE.THE SIX FILES THAT FOLLOW ARE CO2 TRANS- ! MISSIVITIES FOR THE 500-850 CM-1 INTERVAL FOR CO2 ! CONCENTRATIONS OF 330 PPMV (1X) ,660 PPMV (2X), AND ! 1320 PPMV (4X). THE FILES ARE ARRANGED AS FOLLOWS: ! FILE 2 1X,CONSOLIDATED USING B(250) WEIGHTING FCTN. ! FILE 3 1X,CONSOLIDATED WITH NO WEIGHTING FCTN. ! FILE 4 2X,CONSOLIDATED USING B(250) WEIGHTING FCTN. ! FILE 5 2X,CONSOLIDATED WITH NO WEIGHTING FCTN. ! FILE 6 4X,CONSOLIDATED USING B(250) WEIGHTING FCTN. ! FILE 7 4X,CONSOLIDATED WITH NO WEIGHTING FCTN. ! FILES 2,4,6 ARE RECOMMENDED FOR USE IN OBTAINING ! TRANSMISSION FUNCTIONS FOR USE IN HEATING RATE ! COMPUTATIONS;THEY CORRESPOND TO THE TRANSMISSIVITIES ! DISCUSSED IN THE 1980 PAPER.FILES 3,5,7 ARE PROVIDED ! TO FACILITATE COMPARISON WITH OBSERVATION AND WITH OTHER ! CALCULATIONS. ! ! PROGRAM LANGUAGE: FORTRAN 1977,INCLUDING PARAMETER ! AND PROGRAM STATEMENTS.THE PROGRAM IS WRITTEN ON A ! CYBER 170-730.SEE THE SECTION ON LIMITATIONS FOR ! ADAPTATIONS TO OTHER MACHINES. ! ! INPUT UNITS,FORMATS AND FORMAT STATEMENT NOS: ! ! UNIT NO VARIABLES FORMAT STATEMENT NO. TYPE ! 5 P (PURPOSE 1) (5E16.9) 201 CARDS ! 5 PD (PURPOSE 2) (5E16.9) 201 CARDS ! 5 PLM(PURPOSE 2) (5E16.9) 201 CARDS ! 5 NMETHD (I3) 202 CARDS ! 20 TRANSA (4F20.14) 102 TAPE !NOV89 ! ITAPE TRANSA (4F20.14) 102 TAPE !NOV89 ! ! OUTPUT UNITS,FORMATS AND FORMAT STATEMENT NOS: ! ! UNIT NO VARIABLES FORMAT STATEMENT NO. ! 6 TRNFCT (1X,8F15.8) 301 PRINT ! 22 TRNFCT (4F20.14) 102 TAPE ! ! PARAMETER INPUTS: ! A) NLEVLS : NLEVLS IS AN (INTEGER) PARAMETER DENOTING ! THE NUMBER OF NONZERO PRESSURE LEVELS FOR PURPOSE 1 ! OR THE NUMBER OF NONZERO LAYER PRESSURES NEEDED TO ! SPECIFY THE PRESSURE LAYERS(PURPOSE 2) IN THE OUTPUT ! GRID. FOR EXAMPLE,IN PURPOSE 1,IF P=0,100,1000,NLEVLS=2. ! IF,IN PURPOSE 2,PD=0,100,500,1000,THE NUMBER OF NONZERO ! PRESSURE LAYERS=2,SO NLEVLS=2 ! IN THE CODE AS WRITTEN,NLEVLS=40; THE USER SHOULD ! CHANGE THIS VALUE TO A USER-SPECIFIED VALUE. ! B) NLP1,NLP2 : INTEGER PARAMETERS DEFINED AS: NLP1=NLEVLS+1; ! NLP2=NLEVLS+2. ! SEE LIMITATIONS FOR CODE MODIFICATIONS IF PARAMETER ! STATEMENTS ARE NOT ALLOWED ON YOUR MACHINE. ! ! INPUTS: ! ! A) TRANSA : THE 109X109 GRID OF TRANSMISSION FUNCTIONS ! TRANSA IS A DOUBLE PRECISION REAL ARRAY. ! ! TRANSA IS READ FROM FILE 20. THIS FILE CONTAINS 3 ! RECORDS,AS FOLLOWS: ! 1) TRANSA, STANDARD TEMPERATURE PROFILE ! 3) TRANSA, STANDARD TEMPERATURES + 25 DEG ! 5) TRANSA, STANDARD TEMPERATURES - 25 DEG ! ! B) NMETHD: AN INTEGER WHOSE VALUE IS EITHER 1 (IF CO2INT IS ! TO BE USED FOR PURPOSE 1) OR 2 (IF CO2INT IS TO BE USED FOR ! PURPOSE 2). ! ! C) P,PD,PLM : ! P IS A REAL ARRAY (LENGTH NLP1) SPECIFYING THE PRESSURE ! GRID AT WHICH TRANSMISSION FUNCTIONS ARE TO BE COMPUTED FOR ! PURPOSE 1.THE DIMENSION OF P IS IN MILLIBARS.THE ! FOLLOWING LIMITATIONS WILL BE EXPLAINED MORE ! IN THE SECTION ON LIMITATIONS: P(1) MUST BE ZERO; P(NLP1),THE ! LARGEST PRESSURE, MUST NOT EXCEED 1165 MILLIBARS. ! PD IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE PRESSURE ! LAYERS FOR WHICH LAYER-AVERAGED TRANSMISSION FUNCTIONS ARE ! TO BE COMPUTED.THE DIMENSION OF PD IS MILLIBARS.THE LIMITATIONS ! FOR PD ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON ! LIMITATIONS. ! PLM IS A REAL ARRAY (LENGTH NLP2) SPECIFYING THE LAYER-MEAN ! PRESSURES. THE DIMENSION OF PLM IS MILLIBARS. THE LIMITATIONS ! FOR PLM ARE THE SAME AS FOR P,AND ARE GIVEN IN THE SECTION ON ! LIMITATIONS.PD IS READ IN BEFORE PLM. ! ! NOTE: AGAIN,WE NOTE THAT THE USER WILL INPUT EITHER P (FOR ! PURPOSE 1) OR PD AND PLM(FOR PURPOSE 2) BUT NOT BOTH. ! ! ! ! ! LIMITATIONS: ! 1) P(1)=0.,PD(1)=0.,PLM(1)=0. THE TOP PRESSURE LEVEL ! MUST BE ZERO,OR THE TOP PRESSURE LAYER MUST BE BOUNDED BY ZERO. ! THE TOP LAYER-MEAN PRESSURE (PLM(1)) MUST BE ZERO; NO ! QUADRATURE IS DONE ON THE TOP PRESSURE LAYER.EVEN IF ONE IS ! NOT INTERESTED IN THE TRANSMISSION FUNCTION BETWEEN 0 AND P(J), ! ONE MUST INCLUDE SUCH A LEVEL. ! 2) PD(NLP2)=P(NLP1) IS LESS THAN OR EQUAL TO 1165 MB. ! EXTRAPOLATION TO HIGHER PRESSURES IS NOT POSSIBLE. ! 3) IF PROGRAM IS NOT PERMITTED ON YOUR COMPILER, ! SIMPLY DELETE THE LINE. ! 4) IF PARAMETER IS NOT PERMITTED,DO THE FOLLOWING: ! 1) DELETE ALL PARAMETER STATEMENTS IN CO2INT ! 2) AT THE POINT WHERE NMETHOD IS READ IN,ADD: ! READ (5,202) NLEVLS ! NLP1=NLEVLS+1 ! NLP2=NLEVLS+2 ! 3) CHANGE DIMENSION AND/OR COMMON STATEMENTS DEFINING ! ARRAYS TRNS,DELTA,P,PD,TRNFCT,PS,PDS,PLM IN CO2INT. ! THE NUMERICAL VALUE OF (NLEVLS+1) SHOULD BE INSERTED ! IN DIMENSION OR COMMON STATEMENTS FOR TRNS,DELTA, ! P,TRNFCT,PS,PLM; THE NUMERICAL VALUE OF (NLEVLS+2) ! IN DIMENSION OR COMMON STATEMENTS FOR PD,PDS. ! 5) PARAMETER (NLEVLS=40) AND THE OTHER PARAMETER ! STATEMENTS ARE WRITTEN IN CDC FORTRAN; ON OTHER MACHINES THE ! SAME STATEMENT MAY BE WRITTEN DIFFERENTLY,FOR EXAMPLE AS ! PARAMETER NLEVLS=40 ! 6) -DOUBLE PRECISION- IS USED INSTEAD OF -REAL*8- ,DUE TO ! REQUIREMENTS OF CDC FORTAN. ! 7) THE STATEMENT -DO 400 KKK=1,3- CONTROLS THE NUMBER OF ! TRANSMISSIVITY OUTPUT MATRICES PORDUCED BY THE PROGRAM.TO ! PRODUCE 1 OUTPUT MATRIX,DELETE THIS STATEMENT. ! ! OUTPUT: ! A) TRNFCT IS AN (NLP1,NLP1) REAL ARRAY OF THE TRANSMISSION ! FUNCTIONS APPROPRIATE TO YOUR ARRAY. IT IS TO BE SAVED ON FILE 22. ! THE PROCEDURE FOR SAVING MAY BE MODIFIED; AS GIVEN HERE,THE ! OUTPUT IS IN CARD IMAGE FORM WITH A FORMAT OF (4F20.14). ! ! B) PRINTED OUTPUT IS A LISTING OF TRNFCT ON UNIT 6, IN ! THE FORMAT (1X,8F15.8) (FORMAT STATEMENT 301). THE USER MAY ! MODIFY OR ELIMINATE THIS AT WILL. ! ! ************ FUNCTION INTERPOLATER ROUTINE ***************** ! ! ! ****** THE FOLLOWING PARAMETER GIVES THE NUMBER OF ******* ! ****** DATA LEVELS IN THE MODEL ******* ! **************************************************************** ! **************************************************************** 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 TRNS(NLP1,NLP1) DIMENSION 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) LOGICAL , EXTERNAL :: wrf_dm_on_monitor DATA NRTAB/1,2,4/ !*********************************** ! THE FOLLOWING ARE THE INPUT FORMATS 100 FORMAT (4F20.14) 743 FORMAT (F20.14) 201 FORMAT (5E16.9) 202 FORMAT (I3) !O222 203 FORMAT (F12.6,I2) 203 FORMAT (F12.6) ! THE FOLLOWING ARE THE OUTPUT FORMATS 102 FORMAT (4F20.14) 301 FORMAT (1X,8F15.8) ! !CC REWIND 15 !CC REWIND 20 !NOV89 REWIND ITAPE !NOV89 !CC REWIND 22 ! ! CALCULATION OF PA -THE -TABLE- OF 109 GRID PRESSURES ! 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 ! N=25 NLV=NLEVLS NLP1V=NLP1 NLP2V=NLP2 ! READ IN THE CO2 MIXING RATIO(IN UNITS OF 330 PPMV),AND AN INDEX ! GIVING THE FREQUENCY RANGE OF THE LBL DATA !O222 READ (5,203) RATIO,IR !CC IR = 1 !CC READ (5,203) RATIO !O222 *********************************** !***VALUES FOR IR***** ! IR=1 CONSOL. LBL TRANS. =490-850 ! IR=2 CONSOL. LBL TRANS. =490-670 ! IR=3 CONSOL. LBL TRANS. =670-850 ! IR=4 CONSOL. LBL TRANS. =2270-2380 !*** IR MUST BE 1,2,3 OR 4 FOR THE PGM. TO WORK ! ALSO READ IN THE METHOD NO.(1 OR 2) !CC READ (5,202) NMETHD GO TO 621 !NOV89 621 ITAP1=20 621 ITAP1=ITAPE !NOV89 NTAP=1 IF (NMETHD.EQ.2) GO TO 502 ! *****CARDS FOR PURPOSE 1(NMETHD=1) !CC 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 ! *****CARDS FOR PURPOSE 2(NMETHD=2) !CC READ (15,201) (PD(I),I=1,NLP2) !CC 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 ! 503 CONTINUE ! *****DO LOOP CONTROLLING NUMBER OF OUTPUT MATRICES !NOV89 !NOV89 DO 400 KKK=1,3 ICLOOP = 3 IF (IR.EQ.4) ICLOOP = 1 DO 400 KKK=1,ICLOOP !NOV89 ! ********************** IF (NMETHD.EQ.2) GO TO 505 ! *****CARDS FOR PURPOSE 1(NMETHD=1) DO 803 I=1,NLP1 P(I)=PS(I) 803 CONTINUE GO TO 506 505 CONTINUE ! *****CARDS FOR PURPOSE 2(NMETHD=2) DO 804 I=1,NLP1 PD(I)=PDS(I) P(I)=PS(I) 804 CONTINUE ! 506 CONTINUE IA=108 IAP=IA+1 !NOV89 IF (NTAP.EQ.1) READ (20,100) ((TRANSA(I,J),I=1,109),J=1,109) !mp IF (NTAP.EQ.1) READ (ITAPE,100) ((TRANSA(I,J),I=1,109),J=1,109) IF (NTAP.EQ.1) THEN IF ( wrf_dm_on_monitor() ) READ (ITAPE,743) ((TRANSA(I,J),I=1,109),J=1,109) CALL wrf_dm_bcast_bytes ( TRANSA , size ( TRANSA ) * RWORDSIZE ) ENDIF !mp IF (NTAP.EQ.1) READ (ITAPE,100) (tmp(I),I=1,11881 !mp do J=109,1,-6 !mp write(6,697)(TRANSA(I,J),I=5,105,10) enddo ! 697 format(11(f5.3,1x)) !mp !NOV89 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 ! *****THIS IS THE END OF PURPOSE 1 CALCULATIONS IF (NMETHD.EQ.1) GO TO 2872 ! 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 ! *****THIS IS THE END OF PURPOSE 2 CALCULATIONS 2872 CONTINUE ! !+++ WRITE (6,301) ((TRNS(I,J),I=1,NLP1),J=1,NLP1) !CC 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 SUBROUTINE CO2INT !CCC PROGRAM CO2IN1 SUBROUTINE CO2IN1(T20,T21,T66,IQ,L,LP1) ! CO2IN1=CO2INS FOR METHOD 1 ! ********************************************************* ! SAVE DATA ON PERMANENT DATA SET DENOTED BY CO222 *** ! ..... K.CAMPANA MARCH 1988,OCTOBER 1988 ! ..... K.CAMPANA DECEMBER 88 CLEANED UP FOR LAUNCHER ! ********************************************************* 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 !O222 LATEST CODE HAS IQ=1 !CC IQ=4 1011 FORMAT (4F20.14) !CC READ (ITIN,1011) ((CO2PO(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN1,1011) ((CO2800(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN,1011) ((CO2PO1(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN1,1011) ((CO2801(I,J),I=1,LP1),J=1,LP1) !CC READ (ITIN,1011) ((CO2PO2(I,J),I=1,LP1),J=1,LP1) !CC 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) !NOV89 IF (IQ.EQ.5) GO TO 300 !NOV89 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) !NOV89 IF (IQ.EQ.5) GO TO 301 !NOV89 CO2801(I,J) = T21(I,J,2) CO2802(I,J) = T21(I,J,3) 301 CONTINUE !***THE FOLLOWING CODE IS REWRITTEN SO THAT THE RADIATIVE BANDS ! ARE: ! IQ=1 560-800 (CONSOL.=490-850) ! IQ=2 560-670 (CONSOL.=490-670) ! IQ=3 670-800 (CONSOL.=670-850) ! IQ=4 560-760 (ORIGINAL CODE) (CONSOL.=490-850) !NOV89 ! IQ=5 2270-2380 (CONSOL.=2270-2380) !NOV89 ! THE FOLLOWING LOOP OBTAINS TRANSMISSION FUNCTIONS FOR BANDS ! USED IN RADIATIVE MODEL CALCULATIONS,WITH THE EQUIVALENT ! WIDTHS KEPT FROM THE ORIGINAL CONSOLIDATED CO2 TF S. IF (IQ.EQ.1) THEN C1=1.5 C2x=0.5 ENDIF IF (IQ.EQ.2) THEN C1=18./11. C2x=7./11. ENDIF IF (IQ.EQ.3) THEN C1=18./13. C2x=5./13. ENDIF IF (IQ.EQ.4) THEN C1=1.8 C2x=0.8 ENDIF !NOV89 IF (IQ.EQ.5) THEN C1=1.0 C2x=0.0 ENDIF !NOV89 DO 1021 I=1,LP1 DO 1021 J=1,LP1 CO2PO(J,I)=C1*CO2PO(J,I)-C2x CO2800(J,I)=C1*CO2800(J,I)-C2x !NOV89 IF (IQ.EQ.5) GO TO 1021 !NOV89 CO2PO1(J,I)=C1*CO2PO1(J,I)-C2x CO2801(J,I)=C1*CO2801(J,I)-C2x CO2PO2(J,I)=C1*CO2PO2(J,I)-C2x CO2802(J,I)=C1*CO2802(J,I)-C2x 1021 CONTINUE !NOV89 IF (IQ.GE.1.AND.IQ.LE.4) THEN !NOV89 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 !NOV89 ENDIF !NOV89 !O222 ********************************************************* !CC REWIND 66 ! SAVE CDTM51,CO2M51,C2DM51,CDTM58,CO2M58,C2DM58..ON TEMPO FILE !CC WRITE (66) (DCDT10(I,I+1),I=1,L) !CC WRITE (66) (CO2PO(I,I+1),I=1,L) !CC WRITE (66) (D2CT10(I,I+1),I=1,L) !CC WRITE (66) (DCDT8(I,I+1),I=1,L) !CC WRITE (66) (CO2800(I,I+1),I=1,L) !CC WRITE (66) (D2CT8(I,I+1),I=1,L) !CC REWIND 66 !O222 ********************************************************* DO 400 I=1,L T66(I,2) = CO2PO(I,I+1) T66(I,5) = CO2800(I,I+1) !NOV89 IF (IQ.EQ.5) GO TO 400 !NOV89 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 SUBROUTINE CO2IN1 !CCC PROGRAM PTZ - COURTESY OF DAN SCHWARZKOPF,GFDL DEC 1987.... SUBROUTINE CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, & SFULL,SHALF,PPTOP,LREAD,NL,NLP,NLP2) ! ! ** THIS PROGRAM CALCULATES TEMPERATURES ,H2O MIXING RATIOS ! ** AND O3 MIXING RATIOS BY USING AN ANALYTICAL ! ** FUNCTION WHICH APPROXIMATES ! ** THE US STANDARD (1976). THIS IS ! ** CALCULATED IN FUNCTION 'ANTEMP', WHICH IS CALLED BY THE ! ** MAIN PROGRAM. THE FORM OF THE ANALYTICAL FUNCTION WAS ! ** SUGGESTED TO ME IN 1971 BY RICHARD S. LINDZEN. ! ****************************************************************** ! CODE TO SAVE STEMP,GTEMP ON DATA SET,BRACKETED BY CO222 ** ! ....K. CAMPANA MARCH 88,OCTOBER 88 DIMENSION SGTEMP(NLP,2),T41(NLP2,2),T42(NLP), & T43(NLP2,2),T44(NLP) DIMENSION SGLVNU(NLP),SIGLNU(NL) DIMENSION SFULL(NLP),SHALF(NL) ! ****************************************************************** ! !*****THIS VERSION IS ONLY USABLE FOR 1976 US STD ATM AND OBTAINS ! QUANTITIES FOR CO2 INTERPOLATION AND INSERTION INTO OPERA- ! TIONAL RADIATION CODES ! 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) ! ! DATA PROFIL/ & 'US STANDARD 1976'/ DATA PSMAX/1013.250/ ! ! ** NTYPE IS AN INTEGER VARIABLE WHICH HAS THE FOLLOWING ! ** VALUES: 0 =SIGMA LEVELS ARE USED; 1= SKYHI L40 LEVELS ! ** ARE USED; 2 = SKYHI L80 LEVELS ARE USED. DEFAULT: 0 ! NTYPE=0 !O222 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) !*******DETERMINE THE PRESSURES (PRESS) PSTAR=PSMAX ! !*** LTOP COMPUTATION MOVED FROM MODEL INITIALIZATION ! LTOP(1)=0 LTOP(2)=0 LTOP(3)=0 DO 30 N=1,NL PCLD=(PSTAR-PPTOP*10.)*SHALF(N)+PPTOP*10. IF(PCLD.GE.642.)LTOP(1)=N IF(PCLD.GE.350.)LTOP(2)=N IF(PCLD.GE.150.)LTOP(3)=N ! PRINT *,N,PCLD,SHALF(N),PSTAR,PPTOP 30 CONTINUE ! !O222 IF (NTYPE.EQ.1) CALL SKYP(PSTAR,PD,GTEMP) !O222 IF (NTYPE.EQ.2) CALL SKY80P(PSTAR,PD,GTEMP) !O222 IF (NTYPE.EQ.0) CALL SIGP(PSTAR,PD,GTEMP) !CC---- CALL SIGP(PSTAR,PD,GTEMP) NLM=NL-1 CALL SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, & SFULL,SHALF,PPTOP,LREAD,NL,NLP,NLM,NLP2) PD(NLP2)=PSTAR DO 40 N=1,NLP PRSINT(N)=PD(NLP2+1-N) 40 CONTINUE ! *** 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) !********************* DO 100 N=1,NLEV ! ! ** ESTABLISH COMPUTATATIONAL LEVELS BETWEEN USER LEVELS AT ! ** INTERVALS OF APPROXIMATELY 'DELZAP' KM. ! 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*G0*ZNINT) HT=ALT(N) ! ! ** CALCULATE HEIGHT AT NEXT USER LEVEL BY MEANS OF ! ** RUNGE-KUTTA INTEGRATION. ! 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 !mp write(6,*) 'RK values,DZ ', RK1,RK2,RK3,RK4,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 504 CONTINUE !O222 ***************************************************** !***OUTPUT TEMPERATURES !O222 ***************************************************** DO 901 N=1,NLP SGTEMP(N,1) = TMPINT(NLP2-N,1) 901 CONTINUE !O222 ***************************************************** !***OUTPUT GTEMP !O222 ***************************************************** DO 902 N=1,NLP SGTEMP(N,2) = GTEMP(N) 902 CONTINUE !O222 ***************************************************** RETURN END SUBROUTINE CO2PTZ FUNCTION PATH(A,B,C,E) !.... ! 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 FUNCTION PATH !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE QINTRP(XM,X0,XP,FM,F0,FP,X,F) !.... ! 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 QINTRP 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 ! *****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 QUADSR !--------------------------------------------------------------------- SUBROUTINE SIGP(PSTAR,PD,GTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, & SIGLV,SIGLY,PPTOP,LREAD,KD,KP,KM,KP2) 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), & T43(KP2,2),T44(KP) ! integer :: retval ! character(50) :: prsmid='prsmid' !CC 18 LEVEL SIGMAS FOR NMC MRF(NEW) MODEL !CC DATA Q/.021,.074,.124,.175,.225,.275,.325,.375,.425,.497, & !CC .594,.688,.777,.856,.920,.960,.981,.995/ ! FOR SIGMA MODELS,Q=SIGMA,QMH=0.5(Q(I)+Q(I+1), ! PD=Q*PSS,PLM=QMH*PSS.PSS=SURFACE PRESSURE(SPEC.) ! !..... GET NMC SIGMA STRUCTURE !CC IF (LREAD.GT.0) GO TO 914 !--- PPTOP IS MODEL TOP PRESSURE IN CB.... ! SIGMA DATA IS BOTTOM OF ATMOSPHERE TO T.O.A..... !cccc PPTOP=5.0 ! READ(11,PPTOP,END=12321) 12321 CONTINUE ! WRITE(6,88221)PPTOP,KD,KP !88221 FORMAT(' ENTER SIGP PPTOP=',E12.5,' KD=',I2,' KP=',I2) ! open(unit=23,file='fort.23',form='unformatted' & ! , access='sequential') ! REWIND 23 ! READ(23)SIGLY ! DO KKK=1,KD ! SIGLY(KKK)=1.-(FLOAT(KKK)-0.5)/KD ! END DO ! 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 ! DO KKK=1,KP ! SIGLV(KKK)=1.-(FLOAT(KKK-1))/KD ! END DO ! WRITE(6,88223) !88223 FORMAT(' READ ETA') ! PRINT 704,(SIGLY(K),K=1,KD) ! PRINT 704,(SIGLV(K),K=1,KP) ! 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. !--- IF PTOP NOT EQUAL TO ZERO ADJUST SIGMA SO AS TO GET PROPER STD ATM ! 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 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 ! call int_get_fresh_handle(retval) ! close(retval) ! write(0,*)' before open in CO2O3' ! open(unit=retval,file=prsmid,form='UNFORMATTED',iostat=ier) ! write(0,*)' after open1' ! do k=1,62 ! write(retval)pd(k) ! enddo ! close(retval) 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. !+++ WRITE (6,100) (GTEMP(K),K=1,KD) !+++ WRITE (6,100) (PD(K),K=1,KP2) !+++ WRITE (6,100) (PLM(K),K=1,KP) !***TAPES 41,42 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM (PS=1013MB) ! 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 !CC WRITE (41,101) (PD(K),K=1,KP2) !CC WRITE (41,101) (PLM(K),K=1,KP) !CC 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 !***STORE AS PDT,SO THAT RIGHT PD IS RETURNED TO PTZ DO 12 I=1,KP2 PDT(I)=PD(I) 12 CONTINUE !***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 !+++ WRITE (6,100) (PD(K),K=1,KP2) !+++ WRITE (6,100) (PLM(K),K=1,KP) !***TAPES 43,44 ARE OUTPUT TO THE CO2 INTERPOLATION PROGRAM(PS=810 MB) ! 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 !CC WRITE (43,101) (PD(K),K=1,KP2) !CC WRITE (43,101) (PLM(K),K=1,KP) !CC 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 !***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 SIGP !--------------------------------------------------------------------- SUBROUTINE SINTR2 !.... ! 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 ! DETERMINE ETAP,THE VALUE OF ETA TO USE BY LINEAR INTERPOLATION ! 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))/ & (PA(IETAP1)-PA(IETA)) SEXP=SEXPV(IETA)+(PETA-PA(IETA))*(SEXPV(IETAP1)- & 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 SUBROUTINE SINTR2 SUBROUTINE CO2O3(SFULL,SHALF,PPTOP,L,LP1,LP2,RATIN,PSC) !CCC PROGRAM CO2O3 = CONSOLIDATION OF A NUMBER OF DAN SCHWARZKOPF,GFDL ! CODES TO PRODUCE A FILE OF CO2 HGT DATA ! FOR ANY VERTICAL COORDINATE (READ BY SUBROUTINE ! CONRAD IN THE GFDL RADIATION CODES)-K.A.C. JUN89. !NOV89--UPDATED (NOV 89) FOR LATEST GFDL LW RADIATION.....K.A.C. LOGICAL :: opened LOGICAL , EXTERNAL :: wrf_dm_on_monitor CHARACTER*80 errmess ! integer :: retval,kk,ka,kb ! character(50) :: co2='co2' INTEGER etarad_unit61, etarad_unit62, etarad_unit63,IERROR DIMENSION SGTEMP(LP1,2),CO2D1D(L,6),CO2D2D(LP1,LP1,6) !NOV89 DIMENSION CO2IQ2(LP1,LP1,6),CO2IQ3(LP1,LP1,6),CO2IQ5(LP1,LP1,6) !NOV89 DIMENSION T41(LP2,2),T42(LP1), & 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 SFULL(LP1),SHALF(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) !NOV89 ! 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)) !NOV89 ! 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)) ! ! Deallocate before reading. This is required for nested domain init. ! IF(ALLOCATED (CO251))DEALLOCATE(CO251) IF(ALLOCATED (CDT51))DEALLOCATE(CDT51) IF(ALLOCATED (C2D51))DEALLOCATE(C2D51) IF(ALLOCATED (CO258))DEALLOCATE(CO258) IF(ALLOCATED (CDT58))DEALLOCATE(CDT58) IF(ALLOCATED (C2D58))DEALLOCATE(C2D58) IF(ALLOCATED (STEMP))DEALLOCATE(STEMP) IF(ALLOCATED (GTEMP))DEALLOCATE(GTEMP) IF(ALLOCATED (CO231))DEALLOCATE(CO231) IF(ALLOCATED (CDT31))DEALLOCATE(CDT31) IF(ALLOCATED (C2D31))DEALLOCATE(C2D31) IF(ALLOCATED (CO238))DEALLOCATE(CO238) IF(ALLOCATED (CDT38))DEALLOCATE(CDT38) IF(ALLOCATED (C2D38))DEALLOCATE(C2D38) IF(ALLOCATED (CO271))DEALLOCATE(CO271) IF(ALLOCATED (CDT71))DEALLOCATE(CDT71) IF(ALLOCATED (C2D71))DEALLOCATE(C2D71) IF(ALLOCATED (CO278))DEALLOCATE(CO278) IF(ALLOCATED (CDT78))DEALLOCATE(CDT78) IF(ALLOCATED (C2D78))DEALLOCATE(C2D78) IF(ALLOCATED (CO2M51))DEALLOCATE(CO2M51) IF(ALLOCATED (CDTM51))DEALLOCATE(CDTM51) IF(ALLOCATED (C2DM51))DEALLOCATE(C2DM51) IF(ALLOCATED (CO2M58))DEALLOCATE(CO2M58) IF(ALLOCATED (CDTM58))DEALLOCATE(CDTM58) IF(ALLOCATED (C2DM58))DEALLOCATE(C2DM58) ! ALLOCATE(CO251(LP1,LP1)) ALLOCATE(CDT51(LP1,LP1)) ALLOCATE(C2D51(LP1,LP1)) ALLOCATE(CO258(LP1,LP1)) ALLOCATE(CDT58(LP1,LP1)) ALLOCATE(C2D58(LP1,LP1)) ALLOCATE(STEMP(LP1)) ALLOCATE(GTEMP(LP1)) ALLOCATE(CO231(LP1)) ALLOCATE(CDT31(LP1)) ALLOCATE(C2D31(LP1)) ALLOCATE(CO238(LP1)) ALLOCATE(CDT38(LP1)) ALLOCATE(C2D38(LP1)) ALLOCATE(CO271(LP1)) ALLOCATE(CDT71(LP1)) ALLOCATE(C2D71(LP1)) ALLOCATE(CO278(LP1)) ALLOCATE(CDT78(LP1)) ALLOCATE(C2D78(LP1)) ALLOCATE(CO2M51(L)) ALLOCATE(CDTM51(L)) ALLOCATE(C2DM51(L)) ALLOCATE(CO2M58(L)) ALLOCATE(CDTM58(L)) ALLOCATE(C2DM58(L)) IF ( wrf_dm_on_monitor() ) THEN DO i = 61,99 INQUIRE ( i , OPENED = opened ) IF ( .NOT. opened ) THEN etarad_unit61 = i GOTO 2061 ENDIF ENDDO etarad_unit61 = -1 2061 CONTINUE DO i = 62,99 INQUIRE ( i , OPENED = opened ) IF ( .NOT. opened ) THEN etarad_unit62 = i GOTO 2062 ENDIF ENDDO etarad_unit62 = -1 2062 CONTINUE DO i = 63,99 INQUIRE ( i , OPENED = opened ) IF ( .NOT. opened ) THEN etarad_unit63 = i GOTO 2063 ENDIF ENDDO etarad_unit63 = -1 2063 CONTINUE ENDIF CALL wrf_dm_bcast_bytes ( etarad_unit61 , IWORDSIZE ) IF ( etarad_unit61 < 0 ) THEN CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' ) ENDIF CALL wrf_dm_bcast_bytes ( etarad_unit62 , IWORDSIZE ) IF ( etarad_unit62 < 0 ) THEN CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' ) ENDIF CALL wrf_dm_bcast_bytes ( etarad_unit63 , IWORDSIZE ) IF ( etarad_unit63 < 0 ) THEN CALL wrf_error_fatal ( 'module_ra_gfdleta: co2o3: Can not find unused fortran unit to read in lookup table.' ) ENDIF IF ( wrf_dm_on_monitor() ) THEN OPEN(etarad_unit61,FILE='tr49t85', & FORM='FORMATTED',STATUS='OLD',ERR=9061,IOSTAT=IERROR) ENDIF IF ( wrf_dm_on_monitor() ) THEN OPEN(etarad_unit62,FILE='tr49t67', & FORM='FORMATTED',STATUS='OLD',ERR=9062,IOSTAT=IERROR) ENDIF IF ( wrf_dm_on_monitor() ) THEN OPEN(etarad_unit63,FILE='tr67t85', & FORM='FORMATTED',STATUS='OLD',ERR=9063,IOSTAT=IERROR) ENDIF !===> GET SGTEMP AND OUTPUT WHICH USED TO BE ON UNITS 41,42,43,44.... LREAD = 0 ! DO KKK=1,L !JD READ(23)SIGLNU(KKK) ! SIGLNU(KKK)=1.-FLOAT(KKK)/LP1 ! END DO CALL CO2PTZ(SGTEMP,T41,T42,T43,T44,SGLVNU,SIGLNU, & SFULL,SHALF,PPTOP,LREAD,L,LP1,LP2) !--- ETA_PORT: escala de concentracao por pressao (Curtis-Godson): amount x ! broadening ~ r <=> avaliar a tabela 330ppm em pressoes r*p. PSC=ppm/330. T41 = T41*PSC T42 = T42*PSC T43 = T43*PSC T44 = T44*PSC ! call int_get_fresh_handle(retval) ! close(retval) ! open(unit=retval,file=co2,form='UNFORMATTED',iostat=ier) ! do kk=1,2 ! write(retval)(sgtemp(k,kk),k=1,61) ! enddo DO K=1,LP1 STEMP(K)=SGTEMP(K,1) GTEMP(K)=SGTEMP(K,2) ENDDO !===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. ! IR=1,IQ=1 IS FOR COMMON /CO2BD3/ IN RADIATION CODE... ! FOR THE CONSOLIDATED 490-850 CM-1 BAND... !NOV89 ! ICO2TP=61 ICO2TP=etarad_unit61 !NOV89 IR = 1 RATIO = RATIN NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2) IR = 1 RATIO = RATIN NMETHD = 1 CALL CO2INT(ICO2TP,T41,T42,T20,RATIO,IR,NMETHD,L,LP1,LP2) IR = 1 RATIO = RATIN NMETHD = 2 CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2) IR = 1 RATIO = RATIN NMETHD = 1 CALL CO2INT(ICO2TP,T43,T44,T21,RATIO,IR,NMETHD,L,LP1,LP2) !===> FILL UP THE CO2D1D ARRAY ! THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND ! THEIR DERIVATIVES FOR TAU(I,I+1),I=1,LEVS, ! WHERE THE VALUES ARE NOT OBTAINED BY QUADRATURE BUT ARE THE ! ACTUAL TRANSMISSIVITIES,ETC,BETWEEN A PAIR OF PRESSURES. THESE ! ARE USED ONLY FOR NEARBY LAYER CALCULATIONS INCLUDING H2O.. ! IQ = 1 CALL CO2IN1(T20,T21,CO2D1D,IQ,L,LP1) ! do kk=1,6 ! write(retval)(co2d1d(k,kk),k=1,60) ! enddo DO K=1,L CDTM51(K)=CO2D1D(K,1) CO2M51(K)=CO2D1D(K,2) C2DM51(K)=CO2D1D(K,3) CDTM58(K)=CO2D1D(K,4) CO2M58(K)=CO2D1D(K,5) C2DM58(K)=CO2D1D(K,6) ENDDO ! !===> FILL UP THE CO2D2D ARRAY ! THE FOLLOWING GETS CO2 TRANSMISSION FUNCTIONS AND THEIR DERIVATIVES ! FROM 109-LEVEL LINE-BY-LINE CALCULATIONS MADE USING THE 1982 ! MCCLATCHY TAPE (12511 LINES),CONSOLIDATED,INTERPOLATED ! TO THE MRF VERTICAL COORDINATE,AND RE-CONSOLIDATED TO A ! 200 CM-1 BANDWIDTH. THE INTERPOLATION METHOD IS DESCRIBED IN ! SCHWARZKOPF AND FELS (J.G.R.,1985). ! CALL CO2INS(T22,T23,CO2D2D,IQ,L,LP1,1) ! do kk=1,6 ! write(retval)((co2d2d(ka,kb,kk),ka=1,61),kb=1,61) ! enddo DO K1=1,LP1 DO K2=1,LP1 CDT51(K1,K2)=CO2D2D(K1,K2,1) CO251(K1,K2)=CO2D2D(K1,K2,2) C2D51(K1,K2)=CO2D2D(K1,K2,3) CDT58(K1,K2)=CO2D2D(K1,K2,4) CO258(K1,K2)=CO2D2D(K1,K2,5) C2D58(K1,K2)=CO2D2D(K1,K2,6) ENDDO ENDDO ! !NOV89 !===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. ! IR=2,IQ=2 IS FOR COMMON /CO2BD2/ IN RADIATION CODE... ! FOR THE CONSOLIDATED 490-670 CM-1 BAND... ! ICO2TP=62 ICO2TP=etarad_unit62 IR = 2 RATIO = RATIN NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2) CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2) IQ = 2 CALL CO2INS(T22,T23,CO2IQ2,IQ,L,LP1,2) ! do kk=1,6 ! write(retval)(co2iq2(k,1,kk),k=1,61) ! enddo DO K=1,LP1 CDT31(K)=CO2IQ2(K,1,1) CO231(K)=CO2IQ2(K,1,2) C2D31(K)=CO2IQ2(K,1,3) CDT38(K)=CO2IQ2(K,1,4) CO238(K)=CO2IQ2(K,1,5) C2D38(K)=CO2IQ2(K,1,6) ENDDO !===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. ! IR=3,IQ=3 IS FOR COMMON /CO2BD4/ IN RADIATION CODE... ! FOR THE CONSOLIDATED 670-850 CM-1 BAND... ! ICO2TP=63 ICO2TP=etarad_unit63 IR = 3 RATIO = RATIN NMETHD = 2 CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD,L,LP1,LP2) CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD,L,LP1,LP2) IQ = 3 CALL CO2INS(T22,T23,CO2IQ3,IQ,L,LP1,3) ! do kk=1,6 ! write(retval)(co2iq3(k,1,kk),k=1,61) ! enddo ! close(retval) DO K=1,LP1 CDT71(K)=CO2IQ3(K,1,1) CO271(K)=CO2IQ3(K,1,2) C2D71(K)=CO2IQ3(K,1,3) CDT78(K)=CO2IQ3(K,1,4) CO278(K)=CO2IQ3(K,1,5) C2D78(K)=CO2IQ3(K,1,6) ENDDO !--- FOLLOWING CODE NOT WORKING AND NOT NEEDED YET !===> INTERPOLATE DESIRED CO2 DATA FROM THE DETAILED(109,109) GRID.. ! IR=4,IQ=5 IS FOR COMMON /CO2BD5/ IN RADIATION CODE... ! FOR THE 4.3 MICRON BAND... ! NOT USED YET ICO2TP=65 ! NOT USED YET IR = 4 ! NOT USED YET RATIO = 1.0 ! DAN SCHWARZ --- USE 300PPMV RATIO = 0.9091 (NOT TESTED YET)..... ! NOT USED YET NMETHD = 2 ! NOT USED YET CALL CO2INT(ICO2TP,T41,T42,T22,RATIO,IR,NMETHD) ! NOT USED YET CALL CO2INT(ICO2TP,T43,T44,T23,RATIO,IR,NMETHD) ! NOT USED YET IQ = 5 ! NOT USED YET CALL CO2INS(T22,T23,CO2IQ5,IQ) !NOV89 !... WRITE DATA TO DISK.. ! ...SINCE THESE CODES ARE COMPILED WITH AUTODBL,THE CO2 DATA ! IS CONVERTED TO SINGLE PRECISION IN A LATER JOB STEP.. ! NOT USED YET WRITE(66) CO211 ! NOT USED YET WRITE(66) CO218 !NOV89 IF ( wrf_dm_on_monitor() ) THEN CLOSE (etarad_unit61) CLOSE (etarad_unit62) CLOSE (etarad_unit63) ENDIF RETURN 9061 CONTINUE WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr49t85 on unit ',etarad_unit61 write(0,*)' IERROR=',IERROR CALL wrf_error_fatal(errmess) 9062 CONTINUE WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr49t67 on unit ',etarad_unit62 write(0,*)' IERROR=',IERROR CALL wrf_error_fatal(errmess) 9063 CONTINUE WRITE( errmess , '(A49,I4)' ) 'module_ra_gfdleta: error reading tr67t85 on unit ',etarad_unit63 write(0,*)' IERROR=',IERROR CALL wrf_error_fatal(errmess) END SUBROUTINE CO2O3