PROGRAM DELOAD C C WRITTEN BY LIDA A.M. NEJAD, C DEPARTMENT OF MATHEMATICS, C UMIST, C MANCHESTER M60 1QD, C UK. C C C C----------------------------------------------------------------------- C C THE FOLLOWING IS A LISTING OF THE PROGRAM DELOAD WHICH WRITES C CHEMICAL RATE EQUATIONS AND APPROPRIATE CONSERVATION EQUATIONS FOR C A GIVEN REACTION SET. IN THE FOLLOWING EXAMPLE, IT SETS UP THE C RELEVANT EQUATIONS FOR AN INTERSTELLAR CLOUD MODEL CONSISTING OF C 2411 REACTIONS AND 247 SPECIES OF WHICH THE ABUNDANCES OF 13 C CHEMICAL ELEMENTS, INCLUDING 1 FOR ELECTRONS, ARE CALCULATED BY C THE CONSERVATION EQUATIONS. (N.B. THIS MODEL INCLUDES BOTH C NEGATIVE AND POSITIVE IONS). C C----------------------------------------------------------------------- C C INPUT DATA REQUIREMENT C ---------------------- C C (1) A FILE CONTAINING THE REACTION SET. C (2) A LIST CONSISTING OF SPECIES WHOSE CONSERVATION EQUATIONS ARE C REQUIRED TOGETHER WITH PARTICLES SUCH AS PHOTONS. THE CHARACTER C ARRAY 'SPSET' IS USED TO HOLD THIS LIST. A SEPARATE LIST OF C CHEMICAL ELEMENTS, INCLUDING A '+' AND/OR A '-' FOR IONS; THE C CHARACTER ARRAY 'ELSET' IS SPECIFIED TO INCORPORATE THIS LIST. C ===================================== C INPUT PROGRAM FLAGS C ------------------- C ISORT SORTING SPECIFICATIONS; A SPECIES FILE OR A VALID REACTION C FILE CAN BE PROCESSED BY SPECIFYING ONE OF THE FOLLOWING C VALUES: C =0 IN ADDITION TO (1) AND (2) ABOVE, A DATA FILE CONTAINING THE C LIST OF 'N' TIME-DEPENDENT SPECIES IS ALSO GIVEN. NO SORTING C IS PERFORMED. C =1 (1),(2) ARE PROVIDED, THE PROGRAM SELECTS THE EXPLICITLY C TIME-DEPENDENT VARIABLE SPECIES AND WRITES THEM TO A FILE. C =2 THE REACTION SET MAY INVOLVE REDUNDANT SPECIES AND REACTIONS C OR IT HAS BEEN MODIFIED AND REQUIRES RENUMBERING. THE C RELEVENT REACTIONS CORRESPONDING TO PRESPECIFIED SPECIES ARE C SELECTED AND A REDUCED, NUMBERED REACTION SET IS CREATED. C C IMASS THE MOLECULAR MASS CALCULATION SPECIFICATION: C =1 CALCULATION IS REQUIRED C =0 CALCULATION IS NOT REQUIRED C C ISUBS C =1 THE CHEMISTRY INCLUDES TWO DIGIT SUBSCRIPTS C =0 THE CHEMISTRY DOES NOT INCLUDE TWO DIGIT SUBSCRIPTS C C INEGE C =1 NEGATIVE IONS ARE INCLUDED IN THE CHEMISTRY C =0 NEGATIVE IONS ARE NOT INCLUDED IN THE CHEMISTRY C ===================================== C INPUT PARAMETERS C ---------------- C NSPEC NUMBER OF VARIABLES / TIME-DEPENDENT EQUATIONS (WHEN THE EXACT C NUMBER OF SPECIES IS KNOWN) C MAXSP > NSPEC IF THE NUMBER OF SPECIES IS NOT KNOWN (I.E. WHEN THE C EXACT NUMBER OF SPECIES IS NOT KNOWN) C NREAC NUMBER OF REACTIONS C MAXRE > NREAC IF THE NUMBER OF REACTIONS IS NOT KNOWN (I.E. WHEN C ISORT=2 IS CHOSEN) C AMASS REAL ARRAY OF 'NCONS' ELEMENTS, EACH ELEMENT CONTAINING THE C ATOMIC MASS OF THE CORRESPONDING ELEMENTS IN 'ELSET' C NCONS NUMBER OF SPECIES CONSERVED / CONSERVATION EQUATIONS C NELEM NUMBER OF CHEMICAL ELEMENTS AND COSMIC PARTICLES (GIVEN IN C 'SPSET') C NONE NUMBER OF ONE CHARACTER LONG CHEMICAL ELEMENTS INCLUDING A '-' C AND/OR '+' TO ACCOUNT FOR THE IONS (E.G. NONE=8 FROM 'ELSET' C FOR THE FOLLOWING EXAMPLE PROGRAM) C NCHAR MAXIMUM LENGTH OF A SPECIES (IN CHARACTERS) C NHIGH THE HIGHEST TWO DIGIT SUBSCRIPT PRESENT IN THE CHEMISTRY C C A,ALFA, REAL ARRAYS USED FOR COPYING THE CORRESPONDING VALUES TO THE C BETA REDUCED RENUMBERED REACTION FILE C C THE FOLLOWING ARRAYS ARE USED FOR WORKING STORAGE C L,LT,LL,LTL,LP,LTP,SUMM C ===================================== C CHARACTER TYPE VARIABLES C ------------------------ C SPSET ARRAY OF 'MAXSP+NELEM' ELEMENTS EACH OF 'NCHAR' CHARACTERS C HOLDING THE SYMBOLIC NAMES OF THE SPECIES C ELSET ARRAY OF 'NCONS' ELEMENTS EACH OF 2 CHARACTERS HOLDING THE C SYMBOLIC NAMES OF THE CHEMICAL ELEMENTS C RE ARRAY HOLDING THE SYMBOLIC NAMES OF THE REACTANTS C PR ARRAY HOLDING THE SYMBOLIC NAMES OF THE PRODUCTS C INDXR ARRAY; INDEX OF REACTIONS C INDXS ARRAY; INDEX OF SPECIES C INDXC ARRAY; INDEX OF CONSERVED SPECIES C SPTMP ARRAY; TEMPORARY STORAGE USED WHEN ISORT=2, OF SIZE=5 ALLOWING C FOR TWO REACTANTS AND THREE PRODUCTS C X WORKING STORAGE HOLDING TERMS IN A CONSERVATION EQUATION C X1 WORKING STORAGE HOLDING NEGATIVE ION TERMS C LOSS WORKING STORAGE HOLDING THE LOSS TERMS C PROD WORKING STORAGE INITIALLY HOLDING THE PRODUCTION TERMS, THEN C THE CONCATENATED PRODUCTION AND LOSS TERMS C SPT WORKING STORAGE HOLDING A SPECIES FROM THE ARRAY 'SPSET' C THE FOLLOWING VARIABLES REPRESENT A REACTANT:- RE1,RE2 C THE FOLLOWING VARIABLES REPRESENT A PRODUCT:- PR1,PR2,PR3 C XR,IXR WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXR' C XS,XN, WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXS' C XJ C XC,IXC WORKING STORAGE HOLDING AN ELEMENT FROM THE ARRAY 'INDXC' C C THE FOLLOWING CHARACTER VARIABLES ARE USED AS WORKING STORAGE C AND THEIR SIZES DO NOT NEED TO BE CHANGED:- C RENUM,MC,MG,MS,NEX1,NEX2,MULT1,MULT2,NPERL C C THE FOLLOWING CHARACTER VARIABLES REPRESENT THE VARIOUS TYPES OF TERMS C APPEARING IN THE CONSERVATION, LOSS AND PRODUCTION SUMMATIONS: C C TEMP0 = '+Y(I)' C TEMP1 = '+M*Y(I)' RELEVANT TERMS INVOLVED IN CONSERVATION C TEMP2 = '+M*Y(I)' EQUATIONS C C LTEM1 = '-K(R)*Y(I)' C LTEM2 = '-K(R)*Y(I)*D' VARIOUS FORMS OF POSSIBLE DESTRUCTION TERMS C LTEM3 = '-K(R)*Y(I)*X(L)*D' C LTEM4 = '-K(R)*Y(I)*Y(J)*D' C C PTEM1 = '+K(R)*Y(I)' C PTEM2 = '+K(R)*Y(I)*D' C PTEM3 = '+K(R)*Y(I)*X(L)*D' C PTEM4 = '+K(R)*Y(I)*Y(J)*D' VARIOUS FORMS OF POSSIBLE FORMATION TERMS C PTEM5 = '+K(R)*X(L)' C PTEM6 = '+K(R)*X(L)*D' C PTEM7 = '+K(R)*X(L)*X(L)*D' C C WHERE I=J= 3-DIGIT SPECIES INDEX, L=2-DIGIT CONSERVED SPECIES INDEX, C R=4-DIGIT REACTIONS NUMBER INDEX, AND M=1 OR 2-DIGIT FOR COUNTING THE C CHEMICAL ELEMENT'S SUBSCRIPT IN A MOLECULE. C E.G. PTEM3*21, TEMP1*8 AND TEMP2*9 C THE LENGTH OF THE ABOVE PARAMETERS CAN BE MODIFIED IN CHARACTER TYPE C STATEMENTS. C C N.B. THE PRESENCE OF '*D' IN CERTAIN PRODUCTION AND LOSS TERMS ABOVE C IS TO ALLOW FOR THE CALCULATION OF FRACTIONAL ABUNDANCES IN MODELS IN C WHICH DENSITY VARIES (E.G. IN A CIRCUMSTELLAR ENVELOPE OR A C POST-SHOCK REGION). IN SUCH CASES THE VALUE OF THE PARAMETER 'D' MUST C BE SET TO THE 'DENSITY' AT ANY TIME. HOWEVER, THE VALUE OF 'D' MUST C BE SET EQUAL TO 1.0 WHEN THE DENSITY IS ASSUMED TO BE CONSTANT IN C CERTAIN MODELS. (NOTE THAT THE VALUE OF 'D' (I.E. THE DENSITY) NEEDS C TO BE ASSIGNED IN THE SUBROUTINE DIFFUN, WHEN THE GEAR PACKAGE IS C USED; FOR FURTHER DETAILS SEE REFERENCES BELOW). C C----------------------------------------------------------------------- C C GENERAL NOTES AND COMMENTS C -------------------------- C C ALTHOUGH THE PROGRAM ALLOWS FOR UP TO 9999 TWO BODY REACTIONS AND 999 C SPECIES OF UP TO 7 CHARACTERS EACH, ALL THE ASSOCIATED VARIABLES CAN C BE EASILY CHANGED TO THE DESIRED VALUES. IT SHOULD BE NOTED THAT C THERE ARE LIMITATIONS IN TERMS OF THE FORTRAN COMPILER AVAILABLE. FOR C EXAMPLE, STANDARD FORTRAN DOES NOT PERMIT A CHARACTER VARIABLE C EXCEEDING 32767 CHARACTERS. HOWEVER, THIS LIMIT IS SUFFICIENT TO C SATISFY THE REQUIREMENTS FOR A SET OF APPROXIMATELY 10,000 REACTIONS C IN WHICH A SPECIES COULD HAVE UP TO ABOUT 2000 PRODUCTION AND LOSS C TERMS. ALSO NOTE THAT FOR LARGE REACTION SETS, THE EQUATIONS OF C CERTAIN SPECIES ( LIKE H AND C+, FOR EXAMPLE) MAY CONTAIN MORE C THAN 19 CONTINUATION LINES, THE UPPER LIMIT ALLOWED BY STANDARD C FORTRAN. IN SUCH CIRCUMSTANCES, THE APPROPRIATE EQUATIONS NEED TO C BE SPLIT INTO SUB-EQUATIONS. DELOAD COPES WITH THE SPLITTING OF THE C VARIABLE SPECIES. EACH SPLIT EQUATION IS ASSIGNED TO THE VARIABLE C 'YD(J)', WHERE J=1,..NO. OF SPLITS, SAY 5 FOR EXAMPLE, AND FINALLY C YDOT(I)=YD(1)+...YD(5). HOWEVER, SUCH SPLITTING FOR CONSERVATION C EQUATIONS LONGER THAN 20 LINES IS NOT YET CARRIED OUT BY DELOAD AND C SO NEEDS TO BE DONE MANUALLY USING THE SYSTEM EDITOR. C C ===================================== C C AN ABRIDGED SAMPLE RECTION SET USED BY DELOAD IS AS FOLLOWS:- C C 1 H CRP H+ ELECTR 4.60E-01 0.00 0.0 C 2 HE CRP HE+ ELECTR 5.00E-01 0.00 0.0 C ... C 57 H+ NH3 NH3+ H 5.20E-09 0.00 0.0 C 58 H+ C2H2 C2H+ H2 4.30E-09 0.00 0.0 C 59 H+ H2CO H2CO+ H 2.00E-09 0.00 0.0 C ... C 285 C+ C2H4 C3H3+ H 1.02E-09 0.00 0.0 C 286 C+ CH3OH CH3+ HCO 2.08E-09 0.00 0.0 C ... C 999 NH2+ H2O NH3+ OH 1.00E-10 0.00 0.0 C 1000 NH2+ H2O H3O+ NH 2.76E-09 0.00 0.0 C 1001 NH2+ H2O NH4+ O 1.45E-10 0.00 0.0 C ... C 1625 H- HCN CN- H2 3.80E-09 0.00 0.0 C 1626 C- O2 O- CO 4.00E-10 0.00 0.0 C ... C 1718 HEH+ ELECTR H HE 2.00E-07 -0.50 0.0 C 1719 CH+ ELECTR C H 3.00E-07 -0.40 0.0 C ... C 1944 C PHOTON C+ ELECTR 2.16E-10 0.00 2.6 C 1945 NA PHOTON NA+ ELECTR 5.60E-12 0.00 2.0 C ... C 2078 H+ C- H C 2.30E-07 -0.50 0.0 C ... C 2113 H CH C H2 1.73E-11 0.50 2200.0 C 2114 H NH N H2 1.73E-11 0.50 2400.0 C ... C 2411 O CH HCO+ ELECTR 2.00E-11 0.44 0.0 C 9999 C C ===================================== C C AN ABRIDGED SAMPLE SPECIES SET C C Y(1 )=H Y(2 )=CO Y(3 )=H+ Y(4 )=CH ... C Y(6 )=OH Y(7 )=C2 Y(8 )=NO Y(9 )=O2 ... C ... C ... C Y(226)=NA+ Y(227)=HNO Y(228)=CCN Y(229)=C5H ... C Y(231)=S2 Y(232)=O2H Y(233)=H2O2 Y(234)=CH3O ... C C ===================================== C C AN ABRIDGED SAMPLE OUTPUT OF DELOAD IS AS FOLLOWS:- C C X( 1) = TOTAL( 1)-(Y(2 )+Y(4 )+2*Y(7 )+Y(12 )+Y(17 )+Y C * (22 )+Y(25 )+Y(26 )+Y(27 )+2*Y(28 )+3*Y(29 )+Y(30 )+Y( C ... C * 30)+Y(234)+Y(235)) C X( 2) = TOTAL( 2)-(Y(2 )+Y(6 )+Y(8 )+2*Y(9 )+Y(13 )+Y C * (19 )+Y(21 )+Y(24 )+Y(25 )+2*Y(30 )+Y(31 )+Y(33 )+Y(37 C ... C * (234)) C X( 7) = TOTAL( 7)+(Y(3 )+Y(51 )+Y(64 )+Y(70 )+Y(71 )+Y(7 C * 3 )+Y(74 )+Y(75 )+Y(76 )+Y(77 )+Y(78 )+Y(79 )+Y(80 )+Y C ... C * 222)+Y(223)+Y(224)+Y(225)+Y(226))-(+Y(162)+Y(163)+Y(16 C * 4)+Y(165)+Y(166)+Y(167)) C ... C X(13) = TOTAL(13)-(Y(16 )+Y(54 )+Y(55 )+Y(77 )+Y(84 )+Y(1 C * 20)+Y(168)+Y(193)) C YD( 1) = C * +K(7 )*X(6 )*D+K(9 )*X(6 )*D+K(9 )*X(6 )*D+K(12 C ... C * )*Y(35 )*D+K(59 )*Y(3 )*Y(37 )*D C YD( 2) = YD( 1) C * +K(62 )*Y(3 )*Y(38 )*D+K(63 )*Y(3 )*Y(39 )*D+K(64 C ... C * *Y(5 )*Y(64 )*D+K(220 )*Y(6 )*Y(64 )*D C YD( 3) = YD( 2) C * +K(225 )*Y(16 )*Y(64 )*D+K(226 )*Y(11 )*Y(64 )*D+K(235 C ... C ... C YD(15) = YD(14) C * +K(2352)*Y(2 )*Y(15 )*D+K(2399)*Y(28 )*Y(36 )*D+K(240 C ... C * D-K(1686)*Y(163)*D-K(1697)*Y(164)*D)) C YDOT( 1) = +YD(15) C * +(Y( 1)*( C * -K(1705)*Y(166)*D-K(1711)*Y(165)*D-K(1715)*Y(167)*D-K( C ... C * )*Y(6 )*D-K(2403)*Y(9 )*D-K(2404)*Y(24 )*D)) C ... C YDOT( 10) = +K(1772)*Y(175)*X(7 )*D+K(1774)*Y(122)*X C * (7 )*D+K(1809)*Y(184)*X(7 )*D+(Y(10 )*(-K(26 )*Y(3 ) C * *D-K(101 )*Y(51 )*D-K(224 )*Y(64 )*D-K(2181)*X(3 )*D-K C * (2221)*X(2 )*D)) C YDOT( 11) = +K(1757)*Y(119)*X(7 )*D+(Y(11 )*(-K(27 C * )*Y(3 )*D-K(31 )*Y(3 )*D-K(100 )*Y(51 )*D-K(226 )*Y C * (64 )*D-K(427 )*Y(76 )*D-K(428 )*Y(76 )*D-K(619 )*Y(83 C * )*D-K(906 )*Y(100)*D-K(1124)*Y(109)*D-K(1301)*Y(125)* C * D-K(1348)*Y(127)*D-K(1970)-K(2216)*X(2 )*D)) C ... C ... C ... C YDOT(234) = +K(2335)*Y(6 )*Y(34 )*D+K(2373)*Y(9 )* C * Y(34 )*D+K(2386)*Y(34 )*Y(63 )*D+(Y(234)*(-K(2147)*Y(1 C * )*D-K(2204)*X(3 )*D-K(2268)*X(2 )*D-K(2342)*Y(6 )*D C * -K(2377)*Y(9 )*D)) C YDOT(235) = +K(1835)*Y(195)*X(7 )*D+K(2286)*Y(22 )*X C * (4 )*D+(Y(235)*(-K(2201)*X(3 )*D-K(2255)*X(2 )*D)) C C ===================================== C C REFERENCES C C BENNETT, A. 1988. IN 'RATE COEFFICIENTS IN ASTROCHEMISTRY' EDS. T.J. C MILLAR AND D.A.WILLIAMS, KLUWER ACADEMIC PUBL. (DORDRECHT), P339. C NEJAD, L.A.M. AND MILLAR, T.J. 1987, ASTRON. ASTROPHYS.,183, 279. C NEJAD, L.A.M. 1986. PH.D. THESIS, UNIV. OF MANCHESTER. C FARQUHAR, P.R.A. AND MILLAR, T.J., 1993, CCP7 NEWSLETTER NO. 18, 6 C----------------------------------------------------------------------- C PARAMETER(MAXRE=4000,MAXSP=394,NCONS=2,NELEM=5) PARAMETER(ISORT=2,IMASS=0,ISUBS=1,INEGE=1) PARAMETER(NPERL=54,IPL=NPERL-1,LCHAR=18*NPERL,LCHAR1=LCHAR+40, * INUM1=NELEM-1,NHIGH=12,NTWO=NHIGH-9,NTOT=MAXSP+NELEM) PARAMETER(NONE=2,NCHAR=7) CHARACTER*1 NEX1,MULT1,RENUM(0:9) CHARACTER*3 XC1,ELSET(INUM1),SM(INUM1),MS(NTWO),INDXC(INUM1),XC, * NEX2,MULT2 CHARACTER*3 XS,XN,XJ,INDXS(MAXSP),PR3,PR4 CHARACTER*4 XR,IXR,INDXR(MAXRE) CHARACTER*7 RE(2,MAXRE),SPECI(MAXSP),PR(4,MAXRE),RE1,RE2,PR1, * PR2,SPSET(NTOT),SPTMP(6),SPT,DEL*9 CHARACTER MG*10,MC*9,WS*200000, * TEMP1*9,TEMP0*7,TEMP2*10, * LTEM1*9,LTEM2*18,LTEM3*20,LTEM4*18, * PTEM1*16,PTEM2*25,PTEM3*26,PTEM4*15,PTEM5*23, * X1*2000,X*9000,LOSS*200000,PROD*200000 REAL AMASS(NCONS),SUMM(MAXSP),A(MAXRE),ALFA(MAXRE),BETA(MAXRE) INTEGER L(NCONS),LT(NCONS),LL(MAXSP),LTL(MAXSP),LTP(MAXSP), * LP(MAXSP) DATA (SPSET(I),I=1,NELEM)/'ELECTR', * 'CRP','PHOTON','CRPHOT',' '/, * (RENUM(I),I=0,9)/'0','1','2','3','4','5','6','7','8','9'/, * MC/'23456789'/, MG/'123456789'/, * (ELSET(I),I=1,NCONS)/'-','+'/ C C OPEN APPROPRIATE DATA AND OUTPUT FILES C SPECIES FILE IS ON UNIT 10 C REACTION FILE IS ON UNIT 11 C OUTPUT FILE IS ON UNIT 12 C NEW RATE FILE IS ON UNIT 13 C N.B. THE CONTENTS OF THE OPEN STATEMENTS CAN BE MACHINE DEPENDENT OPEN(UNIT=10,FILE="species94.data") OPEN(UNIT=11,FILE="rate95.data") OPEN(UNIT=12,FILE="odes.f") OPEN(UNIT=13,FILE="newrate.data") C******************************* SECTION 1 ***************************** NSPEC = MAXSP IF (ISORT.EQ.0) THEN C ----------------------------- (ISORT=0) ------------------------------ READ(10,1002)(INDXS(J),SPECI(J),J=1,NSPEC) WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC) DO 15 J=1,NSPEC SPSET(NELEM+J) = SPECI(J) 15 CONTINUE READ(11,1001) DO 35 J=1,MAXRE READ(11,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ,ALFAJ, * BETAJ,DEL IF(INDXR(J).EQ.'9999') GO TO 40 35 CONTINUE 40 NREAC=MAXRE ELSE IF (ISORT.EQ.1) THEN C ---------------------------------------------------------------------- C ----------------------------- (ISORT=1) ------------------------------ READ(11,1001) DO 125 J=1,MAXRE READ(11,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ,ALFAJ, * BETAJ,DEL WRITE(12,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4),AJ, * ALFAJ,BETAJ,DEL IF (INDXR(J).EQ.'9999') GO TO 130 125 CONTINUE 130 NREAC = J-1 DO 140 I=1,MAXSP 140 INDXS(I) = INDXR(I) NCODE=NELEM DO 180 IR=1,NREAC DO 170 M=1,2 DO 150 N1=1,NCODE+1 IF (RE(M,IR).EQ.SPSET(N1)) GO TO 160 IF (N1.GT.NCODE) THEN NCODE = NCODE+1 NSPEC = NCODE-NELEM SPECI(NSPEC) = RE(M,IR) SPSET(NCODE) = RE(M,IR) GO TO 160 END IF 150 CONTINUE 160 GO TO (170,180)M 170 CONTINUE 180 CONTINUE WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC) ELSE IF (ISORT.EQ.2) THEN C ---------------------------------------------------------------------- C ----------------------------- (ISORT=2) ------------------------------ READ(10,1002)(INDXS(J),SPECI(J),J=1,NSPEC) WRITE(12,1003)(INDXS(J),SPECI(J),J=1,NSPEC) DO 210 J=1,NSPEC SPSET(NELEM+J) = SPECI(J) 210 CONTINUE IRUSED = 0 READ(11,1001) DO 275 J=1,MAXRE READ(11,1004)IXR,RE1,RE2,PR1,PR2,PR3,PR4,AJ,ALFAJ,BETAJ,DEL SPTMP(1) = RE1 SPTMP(2) = RE2 SPTMP(3) = PR1 SPTMP(4) = PR2 SPTMP(5) = PR3 SPTMP(6) = PR4 DO 255 M=1,6 DO 245 K1=1,NTOT IF (SPTMP(M).EQ.SPSET(K1)) THEN GO TO (255,255,255,255,255,260)M END IF 245 CONTINUE GO TO 275 255 CONTINUE 260 IRUSED = IRUSED+1 RE(1,IRUSED) = RE1 RE(2,IRUSED) = RE2 PR(1,IRUSED) = PR1 PR(2,IRUSED) = PR2 PR(3,IRUSED) = PR3 PR(4,IRUSED) = PR4 A(IRUSED) = AJ ALFA(IRUSED) = ALFAJ BETA(IRUSED) = BETAJ C NUMBER VALID REACTIONS IF (IRUSED.LT.10) THEN INDXR(IRUSED) = RENUM(IRUSED) ELSE IF (IRUSED.LT.100) THEN J1 = INT(IRUSED/10) J2 = IRUSED-J1*10 INDXR(IRUSED) = RENUM(J1)//RENUM(J2) ELSE IF (IRUSED.LT.1000) THEN J1 = INT(IRUSED/100) IR1 = IRUSED-J1*100 J2 = INT(IR1/10) J3 = IR1-J2*10 INDXR(IRUSED) = RENUM(J1)//RENUM(J2)//RENUM(J3) ELSE J1 = INT(IRUSED/1000) IR1 = IRUSED-J1*1000 J2 = INT(IR1/100) IR2 = IR1-J2*100 J3 = INT(IR2/10) IR3 = IR2-J3*10 J4 = IR3 INDXR(IRUSED)= RENUM(J1)//RENUM(J2)//RENUM(J3)//RENUM(J4) END IF IF (IXR.EQ.'9999') GO TO 280 275 CONTINUE 280 NREAC = IRUSED C WRITE VALID REACTION FILE TO UNIT 13 WRITE(13,1001) DO 285 J=1,IRUSED WRITE(13,1004)INDXR(J),(RE(I,J),I=1,2),(PR(I,J),I=1,4), * A(J),ALFA(J),BETA(J) 285 CONTINUE CLOSE(UNIT=13) END IF C ---------------------------------------------------------------------- C LOAD ARRAY MS DO 290 I=1,NTWO I2 = 9+I MS(I) = INDXS(I2) 290 CONTINUE C ------------------- CALCULATION OF MOLECULAR MASSES ------------------ C IF THIS CALCULATION IS NOT REQUIRED, THE FOLLOWING LINES UP TO C STATEMENT 400 AND,ARRAYS SUMM,AMASS,MG AND FLAG 'IMASS' MAY BE C REMOVED FROM THE CODE. C ---------------------------------------------------------------------- IF (IMASS.EQ.1) THEN DO 305 I=1,INUM1 SM(I) = ELSET(I) 305 CONTINUE NMASS = NCONS DO 380 NM = 1,NSPEC SPT = SPECI(NM) I = 0 I1 = 0 J1 = 0 J2 = 0 310 I = I+1+J1+J2 J1 = 0 J2 = 0 I1 = 0 IF (SPT(I:I).NE.' ') THEN DO 360 K=1,NMASS DO 315 M=NONE+1,NCONS IF (I.GT.NCHAR-1) GO TO 320 IF (SPT(I:I+1).EQ.SM(M)) THEN I1 = 1 J1 = 1 GO TO 320 END IF 315 CONTINUE 320 IF ((I+I1).LE.NCHAR) THEN IF (K.GT.NONE) I1=1 IF (SPT(I:I+I1).EQ.SM(K)(1:1+I1)) THEN C C CHECK FOR TWO DIGIT ELEMENTAL SUBSCRIPTS AND TAKE THEM INTO ACCOUNT C IF (ISUBS.NE.0) THEN NEX2 = SPT(I+1+I1:I+2+I1) DO 325 J=1,NTWO J3 = 9+J IF (NEX2.EQ.MS(J)) THEN MJ = J3 J2 = 2 GO TO 350 END IF 325 CONTINUE END IF C C CHECK FOR ONE DIGIT ELEMENTAL SUBSCRIPTS C NEX1 = SPT(I+1+I1:I+1+I1) IF(.NOT.(NEX1.EQ.' '.OR.NEX1.EQ.'+'.OR.NEX1.EQ.'-')) THEN DO 335 J=1,9 IF (NEX1.EQ.MG(J:J)) THEN MJ =J J2 = 1 GO TO 350 END IF 335 CONTINUE END IF MJ = 1 350 SUMM(NM) = SUMM(NM)+AMASS(K)*MJ GO TO 310 END IF END IF 360 CONTINUE GO TO 310 END IF 380 CONTINUE C WRITE(12,1001) WRITE(12,390)(SUMM(J),J=1,NSPEC) 390 FORMAT(5(3X,F6.2,1X)) END IF C----------------------------------------------------------------------- 400 CONTINUE C****************************** SECTION 2 ****************************** C------------------- SET UP CONSERVATION EQUATIONS --------------------- IFLAG = 0 LTMC0 = LEN(TEMP0) LTMC1 = LEN(TEMP1) LTMC2 = LEN(TEMP2) DO 500 K=1,NCONS L(K) = LT(K) DO 460 N1 = 1,NSPEC SPT = SPECI(N1) XS = INDXS(N1) I = 0 I1 = 0 J2 = 0 J1 = 0 410 I = I+1+J1+J2 I1 = 0 J2 = 0 J1 = 0 IF (SPT(I:I).NE.' ') THEN DO 415 M=NONE+1,NCONS IF (I.GT.NCHAR-1) GO TO 420 IF (SPT(I:I+1).EQ.ELSET(M)) THEN I1 = 1 J1 = 1 GO TO 420 END IF 415 CONTINUE 420 IF ((I+I1).LE.NCHAR) THEN IF (K.GT.NONE) I1=1 IF (SPT(I:I+I1).EQ.ELSET(K)(1:1+I1)) THEN C C CHECK FOR TWO DIGIT ELEMENTAL SUBSCRIPTS AND TAKE THEM INTO ACCOUNT C IF (ISUBS.NE.0) THEN NEX2 = SPT(I+1+I1:I+2+I1) DO 430 J=1,NTWO IF (NEX2.EQ.MS(J)) THEN MULT2 = NEX2 TEMP2 = '+'//MULT2//'*Y('//XS//')' L2 = LTMC2+L(K) LT(K) = L2 L1 = L(K)+1 WS(L1:L2) = TEMP2 J2 = 2 GO TO 450 END IF 430 CONTINUE END IF C C CHECK FOR ONE DIGIT ELEMENTAL SUBSCRIPTS C NEX1 = SPT(I+1+I1:I+1+I1) IF(.NOT.(NEX1.EQ.' '.OR.NEX1.EQ.'+'.OR.NEX1.EQ.'-')) THEN DO 440 J=1,9 IF (NEX1.EQ.MC(J:J)) THEN MULT1 = MC(J:J) TEMP1 = '+'//MULT1//'*Y('//XS//')' L2 = LTMC1+L(K) LT(K) = L2 L1 = L(K)+1 WS(L1:L2) = TEMP1 J2 = 1 GO TO 450 END IF 440 CONTINUE END IF TEMP0 = '+'//'Y('//XS//')' L2 = LTMC0+L(K) LT(K) = L2 L1 = L(K) +1 WS(L1:L2) = TEMP0 END IF END IF 450 L(K) = LT(K) GO TO 410 END IF 460 CONTINUE ID = LT(K)+1 WS(ID:ID+60) = ' ' IF (INEGE.EQ.1) THEN IF (ELSET(K).EQ.'-') THEN LCON1 = LT(K)+3 X1 = '-('//WS(1:LT(K))//')' GO TO 500 END IF IF (ELSET(K).EQ.'+') THEN LCON = LCON1+LT(K)+2 X = '+('//WS(2:LT(K))//')'//X1 IFLAG = 1 IB = K-1 IC = IB GO TO 480 END IF IF (IFLAG.EQ.1) THEN IB = K-1 IC = IB GO TO 475 END IF END IF IB = K IC = IB 475 X = '-('//WS(2:LT(K))//')' LCON = LT(K)+2 480 IF (LCON.GT.40) THEN WRITE(12,1005)IC,IB,X(:40) NLINE = (LCON-40)/NPERL IF ((LCON-40).NE.NLINE*NPERL) NLINE=NLINE+1 DO 495 J=1,NLINE J1 = J*40+(4*(J-1)+1)+(J-1)*10 J2 = J1+IPL WRITE(12,1006)X(J1:J2) 495 CONTINUE ELSE WRITE(12,1005)IC,IB,X(:40) END IF 500 CONTINUE C****************************** SECTION 3 ****************************** DO 505 J=1,INUM1 INDXC(J) = INDXS(J) 505 CONTINUE C LTML1 = LEN(LTEM1) LTML2 = LEN(LTEM2) LTML3 = LEN(LTEM3) LTML4 = LEN(LTEM4) C LTML5 = LEN(LTEM5) LTMP1 = LEN(PTEM1) LTMP2 = LEN(PTEM2) LTMP3 = LEN(PTEM3) LTMP4 = LEN(PTEM4) LTMP5 = LEN(PTEM5) C LTMP6 = LEN(PTEM6) C LTMP7 = LEN(PTEM7) IBCK = 5+LEN(XS) C-------------------------- SET UP LOSS TERMS -------------------------- C DO 1000 NS=1,NSPEC XS = INDXS(NS) DO 550 NR= 1,NREAC XR = INDXR(NR) LL(NS) = LTL(NS) DO 510 M=1,2 IF (SPECI(NS).EQ.RE(M,NR)) THEN IF (M.EQ.1) THEN M1=M+1 ELSE M1=M-1 END IF GOTO 530 END IF 510 CONTINUE GO TO 550 530 DO 540 IC=1,INUM1 XC = INDXC(IC) IF (RE(M1,NR).EQ.SPSET(IC)) THEN IF (SPSET(IC).EQ.'PHOTON'.OR.SPSET(IC).EQ.'CRP' * .OR.SPSET(IC).EQ.'CRPHOT') THEN LTEM1 = '-K('//XR//')' L2 = LTML1+LL(NS) LTL(NS) =L2 L1 = LL(NS)+1 LOSS(L1:L2) = LTEM1 C ELSE IF (SPSET(IC).EQ.' ') THEN C LTEM2 = '-K('//XR//')*D' C L2 = LTML2+LL(NS) C LTL(NS) =L2 C L1 = LL(NS)+1 C LOSS(L1:L2) = LTEM2 ELSE LTEM2 = '-K('//XR//')*X('//XC//')*D' L2 = LTML2+LL(NS) LTL(NS) = L2 L1 = LL(NS)+1 LOSS(L1:L2) = LTEM2 END IF GOTO 550 END IF 540 CONTINUE DO 545 J=1,NSPEC XJ = INDXS(J) IF (RE(M1,NR).EQ.RE(M,NR).AND.RE(M1,NR).EQ.SPECI(J)) THEN LTEM3 = '-2*K('//XR//')*Y('//XJ//')*D' L2 = LTML3+LL(NS) LTL(NS) = L2 L1 = LL(NS)+1 LOSS(L1:L2) = LTEM3 GO TO 550 ELSE IF (RE(M1,NR).EQ.SPECI(J)) THEN LTEM4 = '-K('//XR//')*Y('//XJ//')*D' L2 = LTML4+LL(NS) LTL(NS) = L2 L1 = LL(NS)+1 LOSS(L1:L2) = LTEM4 GO TO 550 END IF 545 CONTINUE 550 CONTINUE C----------------------------------------------------------------------- C WRITE LOSS TERMS TO OUTPUT FILE AS A CHECK ON FAILURE - THIS WILL ONLY C BE PERFORMED IF THE 'C' IN COLUMN 1 IS REPLACED WITH A BLANK IN THE C LINES UP TO STATEMENT NUMBER 590. C C IF (LTL(NS).GT.40) THEN C WRITE(12,570)NS,LOSS(:40) C 570 FORMAT(6X,'LOSS(',I3,')= ',A40) C NLINE = (LTL(NS)-40)/NPERL C IF ((LTL(NS)-40).GT.NLINE*NPERL) NLINE = NLINE+1 C DO 580 J=1,NLINE C J1 = J*40+(4*(J-1)+1)+(J-1)*10 C J2 = J1+NPERL-1 C WRITE(12,1006)LOSS(J1:J2) C 580 CONTINUE C ELSE C WRITE(12,570) NS,LOSS(:40) C 590 END IF C----------------------------------------------------------------------- C----------------------- SET UP PRODUCTION TERMS ----------------------- C DO 820 NR=1,NREAC XR = INDXR(NR) LP(NS) = LTP(NS) DO 815 MP=1,3 IF (SPECI(NS).EQ.PR(MP,NR)) THEN DO 750 N1=1,NSPEC XN = INDXS(N1) DO 640 M=1,2 IF (RE(M,NR).EQ.SPECI(N1)) THEN IF (M.EQ.1) THEN M1=M+1 ELSE M1=M-1 END IF GOTO 680 END IF 640 CONTINUE GOTO 750 680 DO 690 IC=1,INUM1 XC = INDXC(IC) IF (RE(M1,NR).EQ.SPSET(IC)) THEN IF (SPSET(IC).EQ.'PHOTON'.OR.SPSET(IC).EQ.'CRP' * .OR.SPSET(IC).EQ.'CRPHOT') THEN PTEM1 = '+K('//XR//')*Y('//XN//')' L2 = LTMP1+LTP(NS) LTP(NS) = L2 L1 = LP(NS)+1 PROD(L1:L2) = PTEM1 C ELSE IF (SPSET(IC).EQ.' ') THEN C PTEM2 = '+K('//XR//')*Y('//XN//')*D' C L2 = LTMP2+LTP(NS) C LTP(NS) = L2 C L1 = LP(NS)+1 C PROD(L1:L2) = PTEM2 ELSE PTEM2 = '+K('//XR//')*Y('//XN//')*X('//XC//')*D' L2 = LTMP2+LTP(NS) LTP(NS) = L2 L1 = LP(NS)+1 PROD(L1:L2) = PTEM2 END IF GO TO 810 END IF 690 CONTINUE C DO 740 J=1,NSPEC XJ =INDXS(J) IF (RE(M1,NR).EQ.SPECI(J)) THEN PTEM3 = '+K('//XR//')*Y('//XN//')*Y('//XJ//')*D' L2 = LTMP3+LTP(NS) LTP(NS) = L2 L1 = LP(NS)+1 PROD(L1:L2) = PTEM3 GO TO 810 END IF 740 CONTINUE 750 CONTINUE C DO 800 JC=1,INUM1 XC = INDXC(JC) DO 765 M=1,2 IF (RE(M,NR).EQ.SPSET(JC)) THEN IF (M.EQ.1) THEN M1=M+1 ELSE M1=M-1 END IF GO TO 785 END IF 765 CONTINUE GO TO 800 785 DO 790 KC=1,INUM1 XC1 = INDXC(KC) IF (RE(M1,NR).EQ.SPSET(KC)) THEN IF (RE(M1,NR).EQ.'PHOTON'.OR.RE(M1,NR).EQ.'CRP' * .OR.SPSET(IC).EQ.'CRPHOT') THEN PTEM4 = '+K('//XR//')*X('//XC//')' L2 = LTMP4+LTP(NS) LTP(NS) = L2 L1 = LP(NS)+1 PROD(L1:L2) = PTEM4 C ELSE IF (RE(M1,NR).EQ.' ') THEN C PTEM6 = '+K('//XR//')*X('//XC//')*D' C L2 = LTMP6+LTP(NS) C LTP(NS) = L2 C L1 = LP(NS)+1 C PROD(L1:L2) = PTEM6 ELSE PTEM5 = '+K('//XR//')*X('//XC//')*X('//XC1//')*D' L2 = LTMP5+LTP(NS) LTP(NS) = L2 L1 = LP(NS)+1 PROD(L1:L2) = PTEM5 END IF GO TO 810 END IF 790 CONTINUE 800 CONTINUE 810 LP(NS) = LTP(NS) END IF 815 CONTINUE 820 CONTINUE C----------------------------------------------------------------------- C WRITE PRODUCTION TERMS TO OUTPUT FILE AS A CHECK ON FAILURE - THIS C WILL ONLY BE PERFORMED IF THE 'C' IN COLUMN 1 IS REPLACED WITH A BLANK C IN THE LINES UP TO STATEMENT NUMBER 870. C C IF (LTP.GT.40) THEN C WRITE(12,840)NS,PROD(:40) C 840 FORMAT(6X,'PROD(',I3,')= ',A40) C NLINE = (LTP(NS)-40)/NPERL C IF ((LTP(NS)-40).GT.NLINE*NPERL) NLINE = NLINE+1 C DO 860 J=1,NLINE C J1 = J*40+(4*(J-1)+1)+(J-1)*10 C J2 = J1+NPERL-1 C WRITE(12,1006)PROD(J1:J2) C 860 CONTINUE C ELSE C WRITE(12,840)NS,PROD(:40) C 870 END IF C----------------------------------------------------------------------- IF (LTL(NS).EQ.0) LTL(NS) = 1 IF (LTP(NS).EQ.0) LTP(NS) = 1 L3 = LTL(NS) L4 = LTP(NS) IF (LTL(NS).EQ.1) THEN PROD(1:1) = ' ' L5 = L4 ELSE IF (LTP(NS).EQ.1) THEN L5 = L3+13 PROD(1:L5) = '+(Y('//XS//')*('//LOSS(1:L3)//'))' ELSE L5 = L3+L4+13 PROD(L4+1:L5) = '+(Y('//XS//')*('//LOSS(1:L3)//'))' END IF C C OUTPUT SECTION C L6 = L5+1 PROD(L6:L6+NPERL) = ' ' IF (L5.LE.LCHAR1) THEN IF (L5.GT.40) THEN WRITE(12,1007)NS,PROD(:40) NLINE = (L5-40)/NPERL IF ((L5-40).GT.NLINE*NPERL) NLINE=NLINE+1 DO 925 J = 1,NLINE J1 = J*40+(4*(J-1)+1)+(J-1)*10 J2 = J1+NPERL-1 WRITE(12,1006)PROD(J1:J2) 925 CONTINUE GO TO 1000 ELSE WRITE(12,1007)NS,PROD(:40) GO TO 1000 END IF ELSE C C THIS ELSE-END IF CLAUSE IS ONLY FOR SUBDIVISION OF EQUATIONS WHOSE C CONTINUATION LINES EXCEED 18 BEFORE OUTPUTTING C I1 = 1 I2 = LCHAR IGROUP = L5/LCHAR IF ((IGROUP*LCHAR).NE.L5) IGROUP = IGROUP+1 DO 990 I=1,IGROUP IF (I.GT.1) I2 = I1+LCHAR LASTL = I2-NPERL+25 DO 935 J=1,26 J1 = LASTL+J IF (PROD(J1:J1+2).EQ.'*(-') THEN J2 = J1-I1-IBCK WS(1:J2) = PROD(I1:J1-IBCK) J1 = J1-IBCK GO TO 940 END IF IF (PROD(J1:J1).EQ.'+') THEN J2 = J1-I1 WS(1:J2) = PROD(I1:J1-1) GO TO 940 END IF IF (PROD(I1:I1).EQ.'+'.AND.PROD(J1:J1).EQ.'-') THEN J2 = J1-I1+2 WS(1:J2) = PROD(I1:J1-1)//'))' GO TO 940 END IF IF (PROD(I1:I1).EQ.'-'.AND.PROD(J1:J1).EQ.'-') THEN J2 = J1-I1+10 WS(1:J2) = '+Y('//XS//')*('//PROD(I1:J1-1)//')' GO TO 940 END IF 935 CONTINUE 940 J3 = J2+1 WS(J3:) = ' ' I1 = J1 LASTG = I-1 IF (I.GT.1) THEN WRITE(12,945) I,LASTG 945 FORMAT(6X,'YD(',I2,') = YD(',I2,')') ELSE WRITE(12,950) I 950 FORMAT(6X,'YD(',I2,') = ') END IF NLINE = J2/NPERL IF ((NLINE*NPERL).LT.J2) NLINE =NLINE+1 DO 955 K =1,NLINE K1 =(K-1)*IPL+K K2 =K1+IPL WRITE(12,1006) WS(K1:K2) 955 CONTINUE IF ((L5-I2).LE.LCHAR) THEN WRITE(12,970) NS,I 970 FORMAT(6X,'YDOT(',I3,') = +YD(',I2,')') L8 =L5-J1+1 NLINE = L8/NPERL IF ((NLINE*NPERL).LT.L8) NLINE =NLINE+1 K3 = J1-1 IF (PROD(J1:J1).EQ.'-') THEN WRITE(12,980)NS 980 FORMAT(5X,'*',8X,'+(Y(',I3,')*(') END IF DO 985 K =1,NLINE K1 = (K-1)*IPL+K K2 =K1+IPL WRITE(12,1006) PROD(K3+K1:K3+K2) 985 CONTINUE GO TO 1000 END IF 990 CONTINUE END IF 1000 CONTINUE C CLOSE(UNIT=10) C CLOSE(UNIT=11) C CLOSE(UNIT=12) STOP 1001 FORMAT(//) 1002 FORMAT(5(3X,A3,2X,A7,:)) 1003 FORMAT(5(1X,'Y(',A3,')=',A,1X,:)) 1004 FORMAT(1X,A4,1X,4(1A7,1X),A3,1X,A3,1PE8.2,1X,0PF5.2,1X, * F8.1,A9) 1005 FORMAT(6X,'X(',I2,') = TOTAL(',I2,')',A40) 1006 FORMAT(5X,'*',8X,A54) 1007 FORMAT(6X,'YDOT(',I3,') = ',5X,A40) END