ESSA TR ERL 177-SDL 15 A UNITED STATES DEPARTMENT OF COMMERCE PUBLICATION ,-»"*■' °'^o^ V ESSA Technical Report ERL 177-SDL 15 U.S. DEPARTMENT OF COMMERCE Environmental Science Services Administration Research Laboratories SIXBIT A Generalized Reaction Kinetics Program G. W. ADAMS L. R. MEGILL c».v oc NiKUlA S/, "^r^ 1 ' ^/: cN^ C =LL^^ ESSA RESEARCH LABORATORIES The mission of the Research Laboratories is to study the oceans, inland waters, the lower and upper atmosphere, the space environment, and the earth, in search of the under- standing needed to provide more useful services in improving man's prospects for survivad as influenced by the physical environment. Laboratories contributing to these studies are: Earth Sciences Laboratories: Geomagnetism, seismology, geodesy, and related earth sciences; earthquake processes, internal structure and accurate figure of the Ecirth, and distribution of the Earth's mass. Atlantic Oceanographic and Meteorological Laboratories: Oceanography, with emphasis on the geology and geophysics of ocean basins, oceanic processes, sea- air interactions, hurricane research, and weather modification (Miami, Florida). Pacific Oceanographic Laboratories: Oceanography; geology and geophysics of the Pacific Basin and margins; oceanic processes and dynamics; tsunami generation, propaga- tion, modification, detection, and monitoring (Seattle, Washington). Atmospheric Physics and Chemistry Laboratory: Cloud physics and precipitation; chem- ical composition and nucleating substances in the lower atmosphere; and laboratory and field experiments toward developing feasible methods of weather modification. Air Resources Laboratories: Diffusion, transport, and dissipation of atmospheric con- taminants; development of methods for prediction and control of atmospheric pollution (Silver Spring, Maryland). Geophysical Fluid Dynamics Laboratory: Dynamics and physics of geophysical fluid systems; development of a theoretical basis, through mathematical modeUng and computer simulation, for the behavior and properties of the atmosphere and the oceans (Princeton, New Jersey). National Severe Storms Laboratory: Tornadoes, squall lines, thunderstorms, and other severe local convective phenomena toward achieving improved methods of forecasting, detecting, and providing advance warnings (Norman, Oklahoma). Space Disturbances Laboratory: Nature, behavior, and mechanisms of space disturb- ances; development and use of techniques for continuous monitoring and early detection and reporting of important disturbances. Aeronomy Laboratory: Theoretical, laboratory, rocket, and satellite studies of the physical and chemical processes controlling the ionosphere and exosphere of the earth and other planets. Wave Propagation Laboratory: Development of new methods for remote sensing of the geophysical environment; special emphasis on propagation of sound waves, and electro- magnetic waves at millimeter, infrared, and optical frequencies. Institute for Telecommunication Sciences: Central federal agency for research and services in propagation of radio waves, radio properties of the earth and its atmosphere, nature of radio noise and interference, information transmission and antennas, and meth- ods for the more effective use of the radio spectrum for telecommunications. Research Flight Facility: Outfits and operates aircraft specially instrumented for re- search; and meets needs of ESSA and other groups for environmental measurements for aircraft (Miami, Florida). ENVIRONMENTAL SCIENCE SERVICES ADMINISTRATION BOULDER, COLORADO 80302 -vW^eNTQ^ '■'(NCi %19^'iv U.S. DEPARTMENT OF COMMERCE Maurice H. Stans, Secretary ENVIRONMENTAL SCIENCE SERVICES ADMINISTRATION Robert M. White, Administrator RESEARCH LABORATORIES Wilmot N. Hess, Director ESSA TECHNICAL REPORT ERL 177-SDL 15 SIXBIT A Generalized Reaction Kinetics Program G. W. ADAMS L. R. MEGILL a 6 SPACE DISTURBANCES LABORATORY BOULDER, COLORADO August 1970 For sale by the Superintendent of Documents, U. S. Government Printing Office, Washington, D. C. 20402 Price $1.00 TABLE OF CONTENTS ABSTRACT 1 1. INTRODUCTION 1 2. USING SIXBIT 3 2.1 Data Input 3 2.2 Program Output 10 2.3 Other Things 13 3. THE INTERNAL WORKINGS OF SIXBIT l4 3.1 Variable Definitions 1^ 3.2 How SIXBIT Sets up the Differential Equations 22 3.3 How SIXBIT Handles Sunlight and Photochemistry 2U 3.k Methods of Solution of the Time- Dependent Differential Equations 2? 3.5 Through the Program in Detail 29 h. CONCLUSION 51 5 . REFERENCES 63 APPENDIX A - SAMPI^ INPUT DECK 6^ APPENDIX B - SAMPLE OUTPUT 66 APPENDIX C - SIXBIT LISTING 75 Digitized by the Internet Archive in 2012 with funding from LYRASIS IVIembers and Sloan Foundation i^^p http://archive.org/details/sixbitgeneralizeOOadam SrXBIT - A GENERALIZED REACTION KINETICS PROGRAM G. W. Adams and L. R. Megill A FORTRAN-IV computer program, SIXBIT, has been written as a general reaction-handling tool. This report describes the use of SIXBIT as a "black-box" program, and describes the program in detail. Key Words: atmosphere, chemistry, computer, differential equations, program, reactions. 1 . BITRODUCTION In the past few years a number of computer programs designed to solve sets of coupled, time-dependent differential equations, describing the behavior of a set of reactants, have been written. See, for example, Keneshea [I962, I963, I967] and Hunt [I966]. The justification for writing yet another of these programs is as follows. Past programs have been relatively inflexible, requiring rather major reprogramming to vary the set of reactants. In addition, it is often difficult to determine the major process operating in a complicated system. SIXBIT attempts to resolve these two difficulties. The first objective is answered by subroutine ALGEBRA, which takes a series of cards describing the reactions and writes the appropriate set of differential equations. Once a basic set of reactions is given to the program, adding or deleting one or two reactions is a trivial task. The second objective is more difficult to answer. An attempt to get out the pertinent information is made by printing out a table giving the major production and destruction mechanisms for each species. These two features make the program a useful tool for ionospheric research and perhaps for some laboratory investigations. Of necessity, a number of approximations have been made in handling the solar flux in the atmosphere. The succeeding pages specify those approximations. In addition, no diffusive processes have been included. There are, of course, a significant number of problems in which this omission will make the results invalid. Extensive measures have been taken both to maintain the accuracy and to keep track of a parameter describing the accuracy of the calcu- lations. The program automatically adjusts the step size as needed to keep the accuracy up. This allows for a great increase in speed, since we find that the step size may change from as short as a few microsec- onds during twilight to steps of the order of an hour in the middle of the day or night. The resulting program has turned out to be quite fast. Another feature of the program is that the sun can be "held" in a single position (e.g., noon) while the program is run to steady state. This allows starting values to be obtained so that diurnal variations may be calculated without a steady drift in value of species that have time constants long with respect to a day. The program is written in FORTRAN IV as used at the Boulder Labora- tories CDC 3800. This language allows extensive indexing of indices which is not allowed on many FORTRAN IV systems. However, at least some versions of FORTRAN V allow this flexibility. No extensive measures have been taken to trap for overflows and underflows, and attempts to use this program on machines allowing exponent ranges of ±38 may result in difficulties The program as presently dimensioned will take 50 reactions, 20 reactants ,and 20 heights. These limitations are not fundamental and can easily be extended by modifying dimension statements. The emphasis in this work, however, is to concentrate on the physical processes operat- ing in a system and it is felt that more reactions than this may be ex- ceedingly difficult to interpret. This report is in two sections. Section 2 describes the use of the program as a "black box." This involves the use of the program without reference to its internal workings. Section 3 describes the program in detail sufficient for the person who may want to modify it. 2. USING SIXBIT 2.1 Data Input (Appendix A is the listing of one of the data decks actually used for the calculations described in Section 3.) A. The first input card contains a single integer; or 1. A tells the program to take all the input data from the data cards. A 1 picks up the species concentrations from a magnetic tape, where they were written out in a previous run, but gets the rest of the data from cards. This restart capability is particularly useful on computer systems with restric- tive running times. B. The second card also contains a to 1. A 1 tells the program to make microfilm plots of the profiles; tells it not to. C. The next three cards contain an alpha-numeric note to yourself so you caxL remember what you are doing. These three cards will be printed out as the first output of the program. D. Reactions . Write down on paper your complete reaction set, and make a list of all the species involved. Number the species consecutively. 1 and 2 are not used as species numbers; these are preempted by the program. 3 is Op and h is Oo. From 5 on, there is no restriction on the species identities. If your reaction set uses Op and/or Oo, use these numbers for them. Otherwise, start your numbering with "5". In either case. Op and Oo profiles will be printed out. Now, rewrite your reaction set in terms of the species numbers, making one list for the regular reactions and a second list for the photoreactions. Each reaction must appear in the form of 3 body _♦ 3 "body, so anywhere you have blanks, enter a "1". As an example, consider the reaction set (1) A+B-^C + D k = 10-10 cm3/sec (2) C+D^A+B k= 10-15 cm3/sec . Assuming that neither Op nor On is part of this set, we would number the species A = 5jB=6, = 7,0=^8. Then the reactions would be -3- (1) 05 + 06 + 01 _ 07 + 08 + 01 (2) 07 + 08 + 01 -> 05 + 06 + 01 . Now express each rate coefficient in the form , / T v^ b/T k = k^ (-) e / where k , A, n, and B are constants, and T is temperature. (The tem- perature iR entered later.) If the rate coefficients are just constants, as in the example, you would have for the first reaction k^ = 10" , A=1.0, n=0,B=0. Note that A must not = 0, since it is used as a divisor. However, with n = 0, any positive value of A will do equally well. The (non-photo-) reactions can now he card -punched, one reaction (with rate coefficient) to a card. There can be any number of reactions up to a maximum of 50. (This number is limited only by the dimension size in the program or the machine memory size and is not fundamental.) In order for the program to know when all the reactions have been read in, the last card in this set must have a "99" punched in columns 1 and 2, with the rest of the card blank. The reaction set will be printed out when the program is run. FORMAT (3(l2,rK),lX,3(l2,lX),6x,2(E8.2,2X),2(E9.2,2X)) SAMPLE ( / IS EDGE OF CARD) /05-06_01__07-08_01 1.00E-10__1.00E+00__+0.00E+00__+0.00E+00 /07-08_01__05_06_01 1.00E-15__1.00E+00_-+0.00E+00__+0.00E+00 E. Phot ore act ions . The photoreactions now exist in terms of the species numbers. For these, the reaction cross sections as a function of wave- length must be specified. SIXBIT works with the wavelength region ^ X <■ 12,000A in lOOA steps, numbered as 1 - - lOOA 2 = 100 - 200A 3 = 200 - 3OOA 120 - 11,900 - 12, 000 A -k- For each reaction, you must specify (a) the number of the first wave- length region with non-zero cross-section, (b) the number of the last region with non-zero cross-section, and (c) the cross-section for each region. To our example reaction set, add the photoreaction C + hv - A + A (07_01_01__05-05-Ol) with a cross-section i^ich is 10" ^cm in the region I85O-2UOOA and is zero elsewhere. Then we would specify I9 and 2k as the first and last -20 -19 wavelength regions, followed by 5 x 10 and five 1 x 10 . The first cross-section has been adjusted to give the average cross-section over the region I8OO-I9OOA. The photoreaction list can contain any number of reactions up to a maximum of 20. It must be terminated with a "99" just like the reaction list was. FORMAT (3(12, IK), IX, 3(12, ]:K)) FORMAT (I3,2X,I3) FORMAT (7(e8.2,2X)) SAMPLE /07-01-01--05-05-01 /019__02l+ /5.OOE-2O--I.OOE-I9--I.OOE-I9--I.OOE-I9--I.OOE-I9--I.OOE-I9 Notice that both here and in the calculation of the photon atten- uation, which is handled separately, lOOA bands are used. This will lead to problems any time the cross -sections have structure, since there is no way to average across such structure and still obey Beer's Law. There are several ways this problem can be handled, but it is left to the user to find a reasonable way to treat any particular case. F. Altitude Information . Altitudes are in kilometers, with the highest one first. There are two ways to get this information in. 1. Constant height interval. This is denoted by the first number (an integer) being 0. The next three numbers are the highest altitude, the lowest altitude, and the height interval. For example, if you want 30-100 km in 5 -km intervals, you punch (to tell it it's constant height interval), then l.OOOE+02, 3.000E+01, 5-OOOE+OO. 2. Non-constant height interval. The initial integer is the number of heights to be given, followed by the altitude values in descending order. For either form of input, a maximum of 20 altitude values is allowed by the present dimension sizes. FORMAT (I2,6(2X,E9.3)) - FIRST CARD FORMAT (6(E9.3,2X)) - ALL OTHER CARDS SAMPLE NO. 1 (CONSTANT STEP SIZE) /00__1 . OOOE+02__3 . 000E+01__5 . OOOE+OO SAMPLE NO. 2 (VARIABLE STEP SIZE) /0U__1 . 000E+02_ _9 . 500E+01_ _9 . 000E+01_ _8 . OOOE+01 G. Temperature Profile . One value of temperature must "be given for each altitude, starting with the highest altitude. If none of the rate coeffi- cients has any temperature dependence, then any non-zero values can be entered here. FORMAT (7(e8.2,2X)) SAMPLE /1.23E+02__1.37E+02__1.^6e+02__1.70E+02 H. Number of Atoms per Species . The easiest way to maintain accuracy in a calculation like this is to calculate all the densities except one, then determine what that one must be to keep the total number of atoms/cm^ a constant. Unless the errors get so large that one of the densities goes negative, you will never know anything about the accuracy this way. A better way Is to calculate all the densities, then multiply them all by a factor which will keep the total n\imber of atoms/cm-^ a constant. As long as this factor is reasonably close to unity at each step, this approach appears to do well at keeping errors from growing. However, this only controls the total number of atoms, so that oxygen can get converted to nitrogen (for example) and still keep the total number of atoms con- stant. Therefore, the procedure in SIXBIT is to determine all the dif- ferent types of atoms which should keep their total numbers/cm-^ constant, such as hydrogen, oxygen, etc. Order these types so that the least abun- dant one is first and the most abundant one is last. For each species. -6- the nimber of each different type of atom is given. Every print-out time, SIXBIT will (l) determine the total number of atoms/cm^ of the least abun- dant type, (2) determine the multiplicative factor required to get the total back to the total that was input, (3) multiply those densities which contain that type of atom by this factor, {k) determine the total number of atoms/cm-^ of the next least abundant type, etc. Normally the most abundant type of atom will be the total of all types of atoms. Each time a new factor is determined and used to multiply one of the sets, species #2 will also be multiplied by the factor. Species #2 starts with a value of 1 at all altitudes, and nothing else happens to it. It gets printed out with the other densities, and thus reflects the total amount of "squeezing" that has occurred at each altitude. FORMAT (12, 2X, 10(11, IK)) SAMPLE /07— 2_1_0_0_1_0_^ I. Time and Miscellaneous Flags . This consists of one card containing 9 numbers . 1. T, the reaction time relative to the time at which the calcu- lations were started. This is normally input as zero, although it is sometimes convenient to change it when a reaction set is being restarted. Its only essential use for you (although not for the program) is to keep track of time for reaction sets with character- istic times of less than, a second. 2. TPRINT, the time interval for which you want the data printed out. This number will automatically be adjusted by the program. If you tell it to print out every second (TPRINT = l.O), then it will print out the results of the calculations every second for 10 seconds. At this point it will examine the average time steps that the program used in the last print-out intervals. If the -7- minimum step-size is greater than TPRINT * lO"^, then TPRINT will be multiplied by 10. TPRINT will continue to be increased by factors of 10 at 100 seconds, 1000 seconds, etc., up to a maximum (see below), as long as the step-size continues to grow (the normal circumstance for most reaction sets) . 3. The solar zenith angle and the photon fluxes are recalculated only at print-out time. For this reason, you do not want TPRINT to get too large. TMAXl is the maximum value that TPRIMT is allowed to have when the sun is above the horizon or when the region of the atmosphere for which you are doing calculations is in total darkness. A reasonable value here is 1 hour (3.6 X lo3) for calculations of diurnal variations. k. TMAX2 is the corresponding upper limit for TPRINT during twilight. One minute is reasonable here, unless you are rimning to steady state. 5. TQUIT is the stopping time for the program (in seconds). Typically seconds) . Typically this might be several days (1 day = S.Gk x 10 6. IFLAG8 chooses between two methods of solution. A chooses a standard fourth-order Runge-Kutta routine with variable step-size. A 1 chooses a partial integration routine. The choice depends on the type of reaction set you are dealing with. Basically, if the rapidly varying species are minor (relatively low concentration) species in a "sea" of relatively stable major species, then the partial integration routine is usually faster. If the major species are also the rapidly changing species, then the Runge-Kutta routine is usualfy fester. Most atmospheric problems will therefore run best with the partial integration routine, while the two-reaction set we have been using for an example will be fastest with Riinge-Kutta. 7. IFLAGl is an integer that controls whether or not species can go into steady state. A 1 here lets species go into steady state when they satisfy the criteria (input later on this card). A does not let species go into steady state. 8. EQUIL is the criterion for letting species go into steady state. If X is the species concentration, then X will go into steady state (assuming it's allowed to go by IFLAGl) if I i < wiL _o 10 is a reasonable value for EQUIL. 9. JOSHUA controls the motion of the sun. If JOSHUA = 0, the sun will go through its normal diurnal variation. If JOSHUA = 1. the sun will remain stationary in the sky. This is used to get steady-state starting values for a reaction set by holding the sun at noon and running to steady state. FORMAT (5(e8.2,2X),2(I1,2X),E8.2,2X,I1) SAMPLE /O.OOE+00__1.00E+00__3.60E+03-_6.00E+01__8.6J+E+05__0_„l__1.00E-08__0 J. Local Time Information . This card contains 1. LAT = latitude in degrees 2. LONG = longitude in degrees 3. TMONTH ^ number of the month h. TDAY = day of the month -9- 5. THOUR = local hour 6. TMIN = local minute 7. TSEC ^ local seconds FORMAT (2(E9.2,2X),5(I2,2X)) SAMPLE /+3.72E+0l__+9.25E+0l__03__05__23__00__00 K. Starting Species Concentration Profiles . The initial value of each species concentration at each altitude can be specified, hut not all of them have to be. If O2 and Oo (#3 and #^) are not input, internal pro- files will be used. However, all altitudes must be in the range 30-200 km to use this feature. Any other species profile not specified will be taken as 5 x 10 times the Op concentration. (This is somewhat faster than starting the calculations with zeroes.) The species member is entered on a card by itself; the altitude profile is entered on the following card. This set must be terminated with a "99" in the first two columns, as before. FORMAT (12) - FIRST CARD - SPECIES NUMBER FORMAT (7(E8.2,2X)) - FOLLOWING CARDS - ALTITUDE PROFILE SAMPLE /OU /9.37E+08__i+.25E+09__8.3^E+09 2.2 Program Output (See Appendix B for a sample) A. Information Output Only at the Beginning 1. The three alphanumeric cards you input will be printed out. 2. A statement telling how many reactions (excluding photo- reactions) and how many species you input is printed, (if -10- your highest niombered species is 13, then you put in 11 species, since species numbering starts with #3.) 3. A list of the reactions and rate coefficients that you input, numbered sequentially in the order you put them in. k. The number of photo-reactions input. 5. A list of the photo-reactions, numbered sequentially from the reaction set. 6. A label which says "Starting Values". 7. The local time as you input it. 8. The amount of computer time used up to this point. 9. The reaction time (t) as you input it. 10. The zenith angle of the svoa. 11. A table labeled "Reaction Terms" containing, at this point, the rate coefficients for the reactions as a function of altitude. 12. A table labeled "Species Densities" containing the profiles you input, plus the density for species #2, which is 1.0 at this point. Information Output Every Print -Out Time 1. Local time, as above. 2. Total computer time used up to this point. •11- 3. T, the reaction time in seconds relative to the starting point . k. The solar zenith angle which was used during the last step, therefore the zenith angle appropriate to the previous time, not to the present time. 5. A table labeled "Reaction Terms." For the reaction A+B + C -. D+E+F with rate coefficient K, the reaction term is K[A] [B] [C] . In other words, the reaction term is the absolute value of the contribution to the time derivative of species A, B, C, D, E, or F due to this particular reaction. Another interpretation is that it is the number of times per second that this parti- cular reaction is occurring. Note that the first entry at each altitude in this table is the average step-size in seconds used by the program over the previous print -out interval. 6. A table labeled "Species Densities". Self-explanatory. A minus sign in front of a density indicates that the species is in steady-state. A minus sign in front of an altitude means that all species at that altitude are in steady state. 7. A table labeled "Major Destruction Reactions". For each species at each altitude, a survey is made of all the reactions to see which is the major destroyer. This reaction number is printed out, followed by a slash and then the numbers of the next three destroying reactions, provided that each of the next three is within a factor of iO of the biggest destroyer. A -12- "-1" indicates that there are no non-zero destroying reactions for that species. A "-0" indicates that there is no such species. 8. A table labeled "Major Production Reactions". As above. 2.3 Other Things The program can maXe microfiljn plots of the altitude profiles at each print-out time. Figure 1 shows a sample plot. For local users, the card in subroutine MICROPLT which gives user's name and extension will have to be changed. For other users, subroutine MICROPLT and the CALL MICROPLT statement in the main program (SIXBIT) will have to be removed. The profiles are written out on tape unit 13 every print-out and then 13 rewinds, so that the last profile is always available. Having the first card input value =1 lets the program read starting profiles, but not the rest of the data, from tape. A dummy read statement is the very first operation in SIXBIT. Therefore, to get started, re- move that initial read statement for the first run so that the tape gets labeled, then replace the read statement so that you are always initially reading a labeled tape. If you find any bugs in the program, please let us know. For people on other machines — We would appreciate learning what machine you are running on and what changes you find need to be made for your machines. This will be helpful to other people using similar machines. -13- 3. THE INTERNAL WORKINGS OF SIXBIT 3.1 Variable Definitions (Appendix C is a complete listing of the program.) A. Indices I - reaction nurriber, goes from 1 to NREAC . II - photoreaction number, goes from 1 to NFREAC . IJZ - reaction rate constant index for reaction I at altitude JZ, goes from 1 to (l)(jz) . ILOOK - searching index used in PHYSICS, goes from 1 to i|. ITCH - type-of-atom index, goes from 1 to ITCHXX. JEK - Runge-Kutta index, goes from 1 to ^. JZ - altitude number, goes from 1 to J ZEND . JZTAB - altitude index used for tabulated values of _0 and . Runs from 1 to 87, corresponding to altitudes of 200-28 km in 2-km steps. LAMDA - wavelength number, goes from 1 to 120 corresponding to wavelength bands 0-100A, 100-200A, etc., up to 11,900-12000A. M - species number, goes from 1 to NSPEC . B. Dimensioned Variables (All contained in COMMON /a/) A(l) - Input, read and used in ALGEBRA, part of the temperature-dependent rate coefficient, K = K0^(TEMP/a)^^N*EXPF(b/tEMP) . Should never be zero, since it's used as a divisor. ARK(JRK) - one of the sets of constants (set in START) used in the Runge- Kutta technique. B(I) - See A(I). BRK(JTiK) - See ARK. CRK(JEK) - See ARK. -lU- F(m) - The siom of all gain terms for Y(m), so that the differential equation is written %^ = F(m) - G(m)^Y(m) where F(m) > 90 ? the absorption is controlled primarily by the density at the point of closest approach to the earth of the sun's rays, which is tJie point at which a = 90 . Hqp and Hq are calculated from the existing profiles of and . Fqp axid Fq are determined from equations given by Swider [I96U]. The flux at a given altitude is determined entirely from absorption coefficients carried by the program. The photo-production (or loss) rates of particular species are determined from the calculated fluxes and from absorption coefficients which must be read into the program. Absorption coefficients, like the flux values, are handled in 100-A bands. Values must therefore be entered so that the average over the lOOA band is correct. A special case exists for Lya (lAMDA = 13). Since most of the energy is at 1215A,the cross section of O2 suitable to 1215A is used in the cross- section value. -25- All the photo calculations are repeated every print-out time, so print-out times must be controlled so that the change in solar zenith angle is not excessive. (See TMAXl and TMAX2 in the input data.) Solar spectra at a variety of altitudes, calculated by SIXBIT, are shown in Figures 5, 6, 7, and 8. -26- 3-^ Methods of Solution of the Time -Dependent Differential Equations There are two methods of solution within SIXBIT, with selection made via a O/l flag in the input data. is a Runge-Kutta solution; 1 is a partial integration technique. In general the partial integration tech- nique works best if the rapidly changing species are minor species which are imbedded in a relatively stable sea of major species (which covers most atmospheric problems) - Runge-Kutta would be best for the simple 2-reaction set we have been using as an example. A. The Runge-Kutta Routine, RIMKUT This nas been taken with no changes from Romanelli [I96U]. B. The Partial Integration Routine, WEIRD The equation for any species can be written as K. Y, . Y^. Y^. - Y^ > K. Y, . Y^. 1 li 2i 31 ra 7^ J Ij 2j whicn says that X^ can be factored out of every loss term and never appears in the gain terms. (This is an approximation if species m appears twice or more on the left-hand side of a reaction and not at all on the right- hand side.) Now rewrite this as dY m dt f - g Y m °m m where f and g__ can be functions of all the Y's except Y m nn ^ n -27- NoWj taking f and g to be constant over a step in time, this can be integrated directly to give YJH) = I YJO) - I ( e-8»" m^ ' / m^ ' g„. I g^ where f and g are evaluated at t = 0. The accuracy of a step is checked by starting with Y (h) and calculating backwards in time (t=H) to get Y^(0) , which should be the same as Y (O). Agreement to 1 part in lO-^ is satisfactory here in most cases. If the agreement is not this good, the step-size is halved and the calculation repeated. If the agreement is satisfactory, then the step forward is recalculated using the average f and g over the interval. Any typical method for solving differential equations (such as Runge- Kutta) makes two approximations; that most quantities stay constant over a step and that the integral can be approximated by a Taylor's series. The partial integration method merely recognizes that there is no reason to approximate the integral since it can be evaluated as it stands. ■28- SIXBIT REACTS 3.5 Through the Program in Detail In this section we will go through the program in the logical (work- ing) order, with special attention to areas which may have to be changed to run on other machines. There are no separate dimension statements. All dimensioned variables are in COMMON/a/, which is carried by the main program and all the sub- routines. COMMON/b/ carries all the non-indexed variables needed in more than 1 subroutine. Type statements (REAL and INTEGER) have been used to allow more reasonable variable names. The comment cards carry enough information to get the data in properly. READ (13,997) ICRUD - This is a dummy read wnich Just insures that the program will always read from unit 13 (which holds the profiles for restart purposes) before it writes on I3. This is necessary on our system to avoid changing control cards between runs. CALL REACTS - Subroutine REACTS reads in the reactions, counts the reactions and the species numbers, and reads and counts the photoreactions. The following describes REACTS. DO 10 11=1,51 10 CGLITINUE - This loop reads in the reactions until it finds a 99 in columns 1 and 2, at which point it jumps to statement 20. If it does not find a 99 in the first 51 reactions, it will print an error message, set IFLAG9 = 1, return to the main program, and terminate. This loop also searches out the largest species number existing in the reactions. 20 KREAC = Il-l - The member of reactions is Il-l since it was on II when it found the 99? and that is not a reaction. -29- REACTS SIXBIT HEIGHTS SIXBIT ALGEBRA SIXBIT START NSPEC = NSPEC-2 - This gets the actual number of species input, not the highest numbered species. (This assumes that you have used all the numbers between 3 and NSPEC. If NSPEC = 10 but 5 is never used as a species number, no harm is done, but the program thinks 5 is a real species.) DO 30 12-1.21 30 CONTINUE - This reads in the phot ore act ions; same game as the regular reactions. RETURN to SIXBIT CALL HEIGHTS - This is a trivial subroutine which sets up the alti- tudes. IF(JELAG6.EQ.0)GO to 100 - JEIAG6=0 says that you have given the program the highest altitude, the lowest altitude, and the step size (in kilometers). JELAG6^0 says you have given it JFLAG6 alti- tudes, the highest one first. RETURN to SIXBIT CALL ALGEBRA - ALBEGRA calculates the rate coefficient array K and sets up the arrays NASRAS, NUSRUS, NTERM, and NETERM. DO 260 I2=1,NREAC 260 CONTINUE - This loop sets up NASRAS and NTERM. The first part, down to statement 200, just eliminates species which appear on both sides of the reaction. RETURN to SIXBIT CALL START - This does most of the initializing needed. The first DATA statement contains the three arrays of constants needed by RUNKUT (ARK, BRK, and CRK) and initializes several arrays. HEAR is set to 1 since WEIRD uses HBAR to get a beginning step-size. The first two species in YHOLD are set to 1, since these are the unit density species that the program uses. -30- START The second and third DATA statements contain FLXiriF, the solar photon flux at infinity in photons/cm^-sec-lOOA, starting with the 0-lOOA band (LAMDA - l) . The fourth and fifth DATA statements contain KABSN2, the absorption coefficient for Ng in cm 2^ starting with the value for the 0-lOOA band. This absorption coefficient is not used at present. If you ever want to do F-region calculations, UV absorption by K2 may become important . The sixth and seventh DATA statements contain KABS03, the absorption p coefficient for Oo in cm . The eighth and ninth contain KABS02, the absorption coefficient for 2 Op in cm The tenth contains TAB02, the tabulated Op densities in molecules/cm-^, given from 200 km to 28 km in 2 km steps. The eleventh DATA statement contains TAB03, the tabulated Oo densities, just like TAB02. DO 10 M1=3,NSPEC 10 CONTINUE - reads in species numbers and the number of atoms /molecule of the different types. Note that one card must be present for every species, but that they do not have to be in order. READ (60,99!+) . . . READ (60,993) ..." the latitude and calendar information gets read in here even though it is not needed until subroutine ZENANG. This is done so that the starting concentration profiles, which do get read in START, can be the last cards in the data deck. By having the data arranged in this way, JFLAG2 can be set to pick up the starting pro- files from cards or magnetic tape without having to rearrange the data deck. ■31- START DO 20 M3-2,NSPEC 20 CONTINUE - This loop reads in the starting profiles. It keeps reading until it finds a "99" for a species number, so the loop is set to go from 2 to NSPEC . This way, even if you give it a profile for every species, it will try to read one more and thus find the "99"- Each time it reads a species number, it looks to see if that species is #3 or #h (O2 or Oo). That way, if Op and/or O3 profiles are not input, it can generate its own from TAB02 and/ or TABOS, since it has to have Op and Oo profiles to do the sun- light calculations. It also set YNEW(M^) =2.0 for each species for which it has profiles, so that any species not specified can be given a profile later of 5 x 10" times the Op profile. This is just to avoid having zeroes for profiles, since the solutions are often quite slow when they have to work with true zeroes rather than just small numbers . PRINT 996 IFLAG9=1 GO TO 9000 - if the previous loop never found a "99" to indicate the end of the profiles, it will print an error message, return to SIXBIT, and terminate the r\in. 30 IF (lUNIT.EQ. 13) REWIND 13 - if the loop found a "99", it will come here and rewind tape unit 13 if it read 13 to get the profiles, IF(IFLAG6.EQ.1 etc. - if both and profiles were input, it's happy and will go to ^0. If the Op and/or Oo profiles were not input, then IF(Z(l) .LE. 200.0 etc. - the altitude range must lie within 30-200 km or else the program cannot construct Op and/or Oo profiles internally. If it is not, it will terminate. i+0 DO 50 JZ3=1,JZEND 50 CONTINUE - This loop sets up the Op profiles to be used for the sunlight calculations ( Y3SUN) from the tabulated values, ■32- START SIXBIT regardless of whether or not a profile was input. If the calcu- lations are "being started in daytime, then this profile will be replaced after the first step by the calculated profile. If the calculations are being started at night, then this avoids having the sunrise period being controlled by a nighttime profile. 12 and II are the index values of the TAB02 values which are at alti- tudes that are required for the calculations. This loop also puts the O2 profile derived from TAB02 into YHOLD if an 0^ profile was not input (i.e., if IFLAG7 = O) . DO 70 JZU=1,JZEND 70 CONTINUE - This does the same thing for the Oo profiles that the previous loop did for the O2 profiles. 80 DO 85 m6=5jNSPEC 85 CONTINUE - This loop makes use of the flag values contained in TNE¥ (set in the "DO 20" read loop) to put in profile = 5 x 10"° times the O2 profile if the profile was not input. DO 110 ITCH=1,10 110 CONTINUE - This loop calculates the total number of atoms/cm-^ of each different type in the input data. CHECK is used to find out how many different types of atoms have been declared. ITCHXX is set to the total number of different types of atoms (actually to the first value of ITCH which is followed by CHECK < 1.0). RETURN TO SIXBIT Everything down to this point will be passed through by the program once per run. The statements after this are in a big loop which will be passed through many times. •33- START ZENANG HANDLE 10 CALL ZENANG - This subroutine calculates the zenith angle of the sun (alpha), prints out the starting values that were input, and calculates the scale heights of O2 and Oo for all the altitudes. IF(IFLAG3.NE.0) go to 10 - IFLAG3 was set to zero in the first DATA statement in START. The first time through ZENANG, ALPHA is calcu- lated, the starting values are printed out, and IFLAG3 is set to 1. The rest of the time, since IFLAG3 is not = 0, printing the starting values is omitted, and the calculation of ALPHA is skipped if the sun is being held fixed in the sky. DEC = . . . - This calculates the declination of the sun in degrees. This equation is a sine fit to the annual variation of the solar declination, and carries the approximation that every month has 30.^ days. PRINT 998 - This puts the "STARTING VALUES" label on the print-out. JCNTRL ^0 - This is a carriage control character that will keep the print-out from skipping to the top of the next page, so that the "STARTING VALUES" label will be in the proper place. THOLD = TPRINT TPRINT =0 - This will leep time from being advanced by the print- out routine as is normally done. THOLD will let TPRINT be recovered after the print-out. TZEN = . . . - This calculates the hour angle of the sun, which is measured in degrees (l hour - 15 degrees) relative to noon. ALPHA = . . . - This calculates the zenith angle of the sun, in degrees away from vertical. CALL HANDLE - This is the subroutine that does all the print-out. TX = . . . •3^- HANDLE ZENMG TSEC = TX - This section adds TPRINT onto the previous calendar time to get the new calendar time. Note that TX, TDAY, THOUR, TMIN, TSEC are all integers by virtue of the type statement. ITSEC = TSEC + 100 ITDAY = TDAY + 100 - on the CDC 38OO, an I format which is too small drops digits from the left. Therefore, 10^ printed in an 12 format comes out an Ok. This lets 8 AM on June 1 be printed out as 06/OI 08/00/00 rather than as _6/_l _8/_0/_0. ITIME = KLOCK(l) TIME = ITIME/IOOO - KLOCK is a system subroutine that reads a clock in the CDC 3800 in milliseconds, relative to the start of the program. The rest of this subroutine is a straightforward print -out of the reaction terms, the species densities, the major destruction reactions, and the major production reactions. Right after statement #1+0 (CONTINUE) is IF( JC]\1TRL.EQ.0) GO TO 9OOO. Since the major destruction and production reactions are not deter- mined until after the calculation has been advanced through a step, their print -out is skipped the first time through. (Note that JCNTRL is used as a carriage control character in FORMAT 999, used in the first print statement.) RETUEN to ZENANG JCNTRL =1 - This carriage control character is now set to 1 and will stay that way throughout the rest of the run. TPRINT = THOLD - The true value of TPRINT is recovered. GO TO 100 - This will skip by the section which starts with statement #10. This section, reached by the first statement in the program {W IFLAG3 • • • ) is J^st the calculation of ALPHA without all the extra statements required the first time through. ■35- ZENMG IF(ALPHA.GT-90.0+SQPTF(z(1))) go to 9000 - The equation ^s ~ 90 + -v/Zg J where a^ is the zenith angle of the sun in degrees and Zg is the height of the shadow of the solid earth in kilometers, is a useful approximation. This statement says to skip all the calculations about solar fluxes if the height of the shadow is above the highest altitude of interest. IF(ALPHA.GT.90.0) go to 300 - if a < 90°, the following loop will take the O2 and Oo profiles to be used for the solar fl\ax attenuation calculaticns (Y3SIM and yUsUN) as the current calculated profiles (YHOID) If a > 90°, it wants to use the profiles that existed at 90° (or the closest it got to 90°)? so it skips the replacement loop. DO 200 JZ1=1JZEWD COWTOTUE - This just puts the Og and 0^ profile into Y3SUN and Z^SUN. The rest of ZENMG calculates the scale heights of O2 and 0^ at the altitudes of interest, since this is what PRODUC needs to calculate the attenuation of the solar photon flux. It does this by fitting the top two points in each profile with an exponential, and deter- mining from this the total amounts of On and 0^ that lie above the top altitude. At each altitude, it then determines the total amount of O2 and O3 from there to the top by doing a simple sum, adds on the amount determined to lie above the top altitude, and then deter- mines what the scale heights have to be at that point for the inte- grated densities above that point to come out right. The exponential form is Y(Z^) = YCZg) e-(Zl-Z2)/H so that, if Z-, is the top altitude and Zp is the next -to -the -top, then H = (Z2-Zi)/ln(Y(Zi)/Y(Z2)) . •36- ZENMG The total number of molecules /cm above Z-j^ is Z =0= ■(Z - Zi)/H X = J Y(Z^) Z^Zi which gives X = HY(Z-j_) Now back to the program. 300 IF(Y3Sm(2).GT.Y3SlOT(l)) GO TO 310 - If the top two points on the profile indicate that the density is increasing with altitude, then the exponential extrapolation of the profile to infinity will give an infinite number of molecules. Therefore, if the density is in- creasing, the program will print out a message (PRINT 997), set IFIAG = 1, return to SIXBIT, and terminate. If the density is de- creasing, it goes to statement 310. 310 IF(Y3SUN(1)-GT-0.0) go to 320 - This avoids a division by zero in determining the scale height for O2. If the top Og density is zero, then H2H0n) =1.0 - The scale height is taken as 1 and it jumps to 330. Otherwise it goes to 320 H2H0II>(1) = (z(2)-Z(11)/LOGF(Y3SU¥(1)/Y3SUW(2)) which is just the equation for the scale height given earlier. 330 IF(Y^SUN(2).GT.YUSIM(1)) GO TO 3^0. 350 H3H0LD(1) = . . . - This is the same as the preceding block, except for 0-3 instead of Og. i+00 imD3 = Y3SUN(l)-)fH2HOUD(l) ADDU = Y^SUN(1)^H3H0LD(1) - ADD3 and ADD1+ are the integrated den- sities of Og and Oo. These statements, using the equation for X given before, initialize the sums. DO Ul0 JZ2=2,JZEND - This starts the integration. •37- ZENANG SIXBIT PRODUC ADD3-ADD3+(Y3SUW(JZ2-l)+y3SUN(jZ2))^(Z(JZ2-l)-Z(JZ2))/2.0 pj)Dk = . . . - These take the total numbers of molecules/cm^ between two altitudes as the linear average of the end-points times the distance between the two altitudes, and adds them onto the sum. H2H0LD( JZ2 ) =ADD3/Y3SUW( JZ2 ) H3H0LD(JZ2) = . . . - This uses the equation X = HYCZ-l) to deter- mine the effective scale height at JZ2. ^10 CONTINUE RETUM to SIXBIT DO 100 JZMAIN=1,JZEND - Here begins the calculation proper. JZMAIN is in COMMON/b/ and will be used in several places to know what altitude is being treated. IJZX ^ JZEND - JZMA.IN - This is part of a compound index that will be needed in RUMUT or WEIRD. DO 30 M1=1,NSPEC 30 CONTimiE - This transfers the species densities at this altitude to a singly-dimensioned array, since it is faster to use a single index than to use a double index. CALL PRODUC - This is the subroutine that uses the scale heights calculated in ZENANG to determine the solar flux at Z( JZMAIN), and then determines the FTERM's for the photoreactions from this flux. CALL Q9EXUN - This is a system subroutine which suspends the under- flow detection, and simply replaces an underflow (a number less than ~10~3^°) by zero without writing an error message. IF(NEREAC.EQ.o) go to 9000 - If there are no photoreactions, this subroutine will be skipped. ALPHAR=ALPHA^Pl/l80.0 - ALPHAR is ALPHA in radians. REARTH = 6371.0 - The radius of the earth in kilometers. •38- PRODUC SIXBIT IF ( ALPHA. GT. 90.0) GO TO 100 - This statement and statement 100 divide the calculations into 3 regions; a < 90°, 90 < a < 90 +-^z", and a > 90 + ■^^. For a < 90°, the calculation of F02 and F03 is done using equation 53 in Swider [196^]. For 90 < ct < 90 + -v/ Z , F02 and F03 are calculated with equation ^7 in Swider. If a > 90 + "VZ J the photo terms are set to zero with the loop at statement 1000. F02 and F03 can "be thought of as sec a modified for the curvature of the earth. Having calculated F02 and F03 for either a < 90° or 90 < a < 90 + '^~z' . you wind up at 300 DO 310 LAMDA1=1,120 310 CONTINUE - This loop calculates the solar photon flux as a function of wavelength, FLUX(LAMDA), using the flujc at infinity (FLXINF), which was given in a DATA statement in START, the absorption coef- ficients for Oo and 0^ (KABS02 and KABS03) which were given in DATA Oo (KABS02 and KABS03) which were given in statements in START, the 0^ and 0, densities (Y3SUW & yUSIM) which were set up in ZENANG, and F02 and F03 which were just calculated. DO 330 II=1,NFREAC 330 CONTINUE - The rate that a photoreaction of the form A + hu -» B + C goes is [A] / :^(X)0(\)C3A . dA dt o In SIXBIT, the integral is called FRATE. This loop just calculates FRATE. GO TO 9000 CALL R9EXUN - This restores the system underflow detection. RETURN to SIXBIT HCOIMC = - This is a counter that will he used later. -39- SIXBIT RUNKUT IF(IFLAG8.EQ.1) go to h^ - IFLAG8 is the input flag (O or l) that picks the method of solution, either RIMKUT (the Runge-Kutta solu- tion) or WEIRD (the partial integration solution) . RIMKUT: CALL Q,9EXOT - This suspends underflow detection, as in ZEMWG. H=HBAR(JZMAIN)^3.0 - HBAR was set to 1.0 in START. Later HEAR will be the average step-size over the previous print interval. This statement assumes that the biggest allowable step is probably bigger than the previous average. H is the step to be used in the calculation. 10 TCIiECK = T-TOID+H IFLAG2 ==1 - If a step of size H will put the time beyond the time for the next print-out, then H is reduced to get time just to the print time, and IFLAG2 is set to 1. This will return the program to SIXBIT at the end of the step. 20 DO 30 Ml=2, NSPEC 30 CONTIMJE - This holds onto the current values of Y and QRK in case the error is too big in the next step and the calculation must be repeated. ^0 DO 90 JRK=l,i^- - This starts the h calculations involved in a fourth-order Runge-Kutta solution. DO 50 I1-1,NREAC 50 CONTIMJE - This calculates the TERM for each reaction (see section 2). DO 52 IT1-1,NFREAC 52 CONTimJE - This calculates the FTERM's for the phot ore act ions. JRKX = (JIlK-l)^NSPEC - Part of a compound index to be used later for KRK. -i^O- RUNKUT DO 70 M2=3,NSPEC - This staxts the loop that calculates the time derivatives for the species. MJRK = JRKX+M2 - The rest of the compound index to be used later for KRK. IIIX=(M2-3)^NREAC - Part of a compound index to be used for NASRA.S. DY=0.0 - DY is the time derivative for species M2. It will be transferred to an indexed variable after all the TERM'S have been added on. MST0P-NTERM(M2) - This is the number of TERM'S in NASRAS which affect species M2. This statement is made because NTERM has to be the end of a DO loop specification, and indexed variables cannot be used for that. DO 60 ]M=1,]MST0P - This begins the loop to pick the TERM'S for M2 out of NASRAS and add them onto DY. IF(MSRAS(lirx:+M).LE.NREAC) GO TO 55 - This picks a reaction number from NASRAS and decides if it is a production reaction or a destruction reaction. If it is >NREAC, it is a production reaction. and it goes to DY=DY4^ERM(NASRAS(IIX+NN)), which adds the TERM for that reaction to DY. If the reaction number is NTX, it is looking at photoreactions. If NX > NEREA.C, then it is looking at a producer of Ml. It thus goes to statement #300- IF(]M.GE.WTX.AM).KX.GT.WREAC) go to 300 - If M ^ mX, it is look- ing at regular reactions. If NX > KREAC, then it is looking at a producer of Ml. It thus goes to statement #300. DO 200 IL00K3=1,IIX 200 CONTINUE - This looks through the list of the destruction reactions already found to see if this reaction has been entered previously. If it has, it jumps to statement #500. Otherwise it winds up at 250 IF(NT^.GT.NTX) go TO 260 - This statement and the nine following (through the third GO TO 500) compare the TERM (or FTERM) to the biggest TERM found so far (lIMAX) . If this TERM (or FTERM) is bigger than TLMAX, then TLMAX is set equal to this TERM, the reaction number (NX) (or NX+NREAC in the case of photoreactions) is stored in ILMAX, and it goes to statement -1+7- PHYSICS SIXBIT HAM)LE 300 IF(m.GT.NTX) GO TO 310 500 CONTIKTUE - This section handles production reactions just as destruc- tion reactions were handled in the preceding section. The biggest TEEM is stored in TGMAX, and the reaction number is IGMAX. IF(IL00K2.EQ.1) GO TO 550 - If this is the first time through, it goes to 550 ILOSS(l)=ILMAX BIGG=TGMAX - Which puts IIMAX into ILOSS(l), IGMAX into IGAIN(i), and holds the TERM values in BIG2 and BIGG for later use. If IL00K2 is greater than 1, it does IF(TLMAX.GE.0.1^BIG2) . . . GO TO 600 - Which enters those reactions into the appropriate ILOSS and IGAIN only if they are within XIO of the maximum values . Having finished the IL00K2 loop, it goes to IF(ILOSS(i).EQ.0) go to 610 - If IK)SS(l)=0), then there were no loss processes for Ml, and -1 will be entered into the final arrays with statement #6l0. Otherwise, the final arrays (ILOSSA, ILOSSB, ILOSSC, and ILOSSD) are filled with the values of ILOSS +100, which gives leading zeroes as in HANDLE. 620 IF(IGAIN(i).EQ.0) . . . 630 ... - Same as above for the gain processes. 700 CONTINUE - Return to get another altitude. 800 CONTINUE - Return to get another species. RETURN TO SIXBIT CALL HANDLE - This subroutine, which prints out the TERM's and RTERM's the species densities, and the major production and loss reactions determined in PHYSICS, was discussed earlier. IF(JJFLAG.NE.I) GO TO 175 - JJFLAG was input as or 1. If it was 0, the next statement is skipped. If it was 1, it goes through -lt8- 1 SIXBIT MICROPLT CALL MICROPLT - This subroutine generates microfilm plots of the species densities. Samples are given in Appendix D. No details of the workings of this subroutine will be given. Anyone not working on the Boulder CDC 38OO will have to delete this sub- routine and the call to it. For people using the CDC 38OO, note that the name and extension in the first DATA statement must be changed so that we do not get your microfilm. RETURN TO SIXBIT 175 TOLD=T - TOLD is advanced in preparation for the next set of calculations. IF(TOLD.GE.TQUIT) go to 800 - TQUIT was the time to quit which was input. If T < TQUIT, it goes to the next section, which in- creases the print interval (TPRINT). This section works well only if T was input as 0.0 initially. In normal operation, the result of this section is to have the print- outs come at T=l, 2, 3, • • •, 9, 10, 20, 30, • • •, 90, 100, 200, 300, • • •, 900, 1000", 2000, etc., which keeps you from being buried in paper while still giving regions of rapid change in sufficient detail to see what is going on. IE'LAG^=IFLAGU+1 - IFLAG^ is a counter which is treated as the second significant digit of T. (This is true if T starts at 0.0 and if TPRINT is 1 X 10^.) IF (IFLAGU.lt. 10) GO TO I88 - The adjustment of TPRINT is done every tenth print-out. When IFLAGU = 10, then it sets it back to and advances IFLAG5 by 1. IFLAG5 is treated as the first significant digit, so it gets reset to zero when it reaches ten. HMIN=TPRINT l80 CONTINUE - This determines HMIN, the minimi;an average step-size (HBAR) over the previous print interval. IF (HMIN/TPRINT.LT. 0.001) GO TO I88 - If more than 1000 steps were required at any altitude during the previous print interval, then -k9- SIXBIT the next section is skipped, i.e., TPRINT is not adjusted. If less than 1000 steps were required, then TPRINT will be adjusted, so it goes to IF(ALPHA.LT. 90.0. OR. ALPHA. GT. 90. 0+SQPTF(Z(l))) GO TO l82 - If a < 90 or a > 90 +'\/Z , then the altitudes of interest are either in total daylight or total darkness. If neither is true, then the region is in twilight. The following section, i.e., IF(TPRINT^10.0.LT.TMAK2) GO TO l8J+ 18^ TPRINT =TPRINT-x-10.0 - Increase TPRINT either by XIO or up to TMAKl (whichever is smaller) for daylight or dark conditions, or increase TPRINT either by XIO or up to TMAX2 (whichever is smaller) for twi- light conditions. (TMAXl and TMAX2 are the input upper limits for TPRINT ) . Then IFLAGi+=IFLAG5 IFLAG5=0 GO TO 190 - which shifts the identities of the first and second significant digits in T, and goes to statement #190. 188 IF ( ALPHA. GT. 90. AND. ALPHA. LT. 90. 0+SQPTF(Z(l)). AND. TPRINT. GT.TMAX:2) TPRINT =TMAX2 - This statement was reached either when it was not yet time to check on the adjustment of TPRINT (3 statements after #175) or when the average step-size was too small to increase TPRINT (one statement past #l80) . In either case, this statement checks to see if the sun has entered the twilight region, and if it has, de- creases TPRINT if need be down to TMAK2. 190 DO 210 JZ^=1,JZEND 210 CONTINUE - YHOLD is shifted into YLAST so that YLAST will be avail- able at the next print-out for checking on steady-state. DO 220 JZ5=1,JZEND IF(JPULL(JZ5).EQ.0) GO TO 10 -50- SIXBIT 220 CONTINUE - JZPULL(JZ5) was set =1 if all the species at that alti- tudes were at steady-state. If any of the JZPULL's are not =1, then it goes back to statement #10 to start another set of calculations. If all the JZPULL's are :=1, then it winds up at 700 PRINT 999 GO TO 9000 - Which prints out a statement that "ALL SPECIES AT ALL ALTITUDES ARE NOW IN EQUILIBRIUM" and quits. h. CONCLUSION In this report we have given the details of the use of a computer program designed to give solutions of sets of reactions. The main virtue of the program is its generality and the ease with which different reactions can be used. ■51- 100 80 X 60 X CO 20 / \ \ l\ \ O' 10^ 10' 10* 10' 10' 10' 10' 10* 10" 10" 10'^ 10" 10" 10" 10'* 10' ALPHA= 22.9 12: OHO MCNTH 6 DAY 1 Figure 1. Sample miarofilm plot. 52 10" : I0l6 1015 10'^ ' / — --— — — — / [ / : 8 - E o r - : ^ ; U- : lO'O m 1 " - n 1 : I09 i ; 10' - in? " 1 1 1 1 1 1 1 1000 2000 3000 4000 5000 6000 7000 8000 9000 lOOOO 11,000 12,000 Wavelength, A Figure 2. Solar photon flux outside the atmosphere data are also tabulated in table 1.) CThei 53 ! ' ' ' ' ^Ionization Continuum 10-17 f\ ^Schunnann-Runge Bands n 1 1 1 J^ 1 Absorpfi )n Cross for Oz -Section 10-18 - : 10-19 : Ly a : 'e 5 in-21 k! '0 i 10-22 : Herzb erg Ban is - 10-23 1 ; 10-2^ ■ : ; r . ^ Absa rption C ross-Se tion for Nz: - : 1 — 1 : ■ ; - ; ; 2000 3000 4000 5000 Wavelength, & Figure 3. Absorption cross sections for N2 and O2. (Thi data also tabulated in table 1.) 54 Ionization Continuum „. riley Be nds ' ' : ] J [ 1 1 Abs 3rption ( ;ross-Se ction fo r O3 ^1 \ 1 / /Chapp \ is Bond s / 1 \ \ . , , 1 10 00 20 00 30 00 4C 00 50 00 60 00 70 00 80 00 9C 00 10^ 00 ll,C 00 I2,C 00 Wavelength, A Figure 4. Absorption aross seotions for 0^. (These data are also tabulated in table 1.) 55 I0 \ I0'5 10'^ 1013 I0l2 8 10" loio c 10^ 10^ —Flux at 30 km — Flux at Infinity \ r-T- Flux at 40 km Flux at Infinity 1000 1500 2000 2500 3000 3500 4000 1000 1500 2000 2500 3000 3500 4000 Figure 5. Solar -photon flux at 20 and 40 km in the Earth'i atmosphere for an overhead sun. Also shown is the flux outside the atmosphere from figure 2. 56 I o o 1012 10". IQlO 10^ : 1 1 1 1 : / / : ~- r " r-- E r 1 ; r-i 1 uJ : 1 1 ' III - 1 LJ 1 1 i J ; ; : : ^ - ^ ; ' F "lux at 50 km " ! F 'lux at Infinity = ' ^ : / i : ^ ^ : r ; \ " r 1 1 J : -1 1 s : - \ : : I : z ^ Fk Flu X at 6( X at In Dkm finity = 1000 1500 2000 2500 3000 3500 4000 1000 1500 2000 2500 3000 3500 4000 Figure 6. Solar photon flux at 50 and 60 km in the Earth 'i atmosphere for an overhead sun. Also shown is the flux outside the atmosphere from figure 2. 57 1016 I0'5 10'^. 1013 1012. o lo'i loio iqs . : / ' : / : : r 1 ^ 1 1 J : ; 1 : 1 1 1 Li : \ \ : : : : F 'lux at ■|ux at 70 km ' nfinity": '- 1 1 III., F 1 1 1 1 : 1^-^ ^ : / -. y" \ r-l 1 : r -1 1 : "l. -J 1 -J : : : : : -. : : ^ \ Flu Flu ,11, X at 8 X at In 1,11 km finity : 1,11 : 000 15 00 >0od"25 00 30 00 35 00 40C Figure 7. Solar photon flux at 70 and 80 km in the Earth'i atmosphere for an overhead sun. Also shown is the flux outside the atmosphere from figure 2. 58 10^ : _^ r ^ / : y : ^ ^ ; : H ^■^ : i r" 1 -J ' - M ' '- - : ^ F 'lux at ■|ux at 90 km ' nfinity j ^ r : / : : y ^ 1 r \ : n 1 -J - n ■-' r- \ ir \ : - : - : - : FIl Flu X at 10 X at In km finity : : - 1500 2000 2500 3000 3500 4000 1000 1500 2000 2500 3000 3500 4000 Figure 8. Solar photon flux at 90 and 100 km in the Earth 'i atmosphere for an overhead sun. Also shown is the flux outside the atmosphere from figure 2. 59 Table 1. Internal Data' Wavelength Solar Photon Flux Absorption Coefficients (cm^/lOOA) Band( KJ Outside the Atmosphere N o_ 0^ (photons/cm2-sec-100A) 2 2 3 0-100 3.79(8) 5.0( -18) 6.o(-i9) 1.5(-20) 100-200 2.65(9) 1.1( -17) 6.o(-l8) 3.0(-20) 200-300 3.80(9) 1.2( -17) i.M-17) 5.0(-20) 300-1^00 8.83(9) iM -17) 2.0(-17) 9.0(-20) UOO-500 3.52(9) 1.5( -17) 2.3(-17) l.7(-19) 500-600 7.71(9) 1.6( -17) 2.3(-17) 3.2(-l9) 600-700 5.08(9) 1.8( -17) 2.1(-17) 6.o(-i9) 700-800 5.88(9) 1.9( -17) 1.5(-17) i.o(-i8) 800-900 1.20(10) 1.0( -18) 9.0(-l8) 2.0(-l8) 900-1000 2.16(10) 1.0( -18) 3.8(-l8) 3.6(-l8) 1000-1100 1.75(10) 1.0( -18) 1.1(-18) 7.0(-l8) 1100-1200 1.03(10) 2.0( -21) 3.M-19) 1.0(-17) 1200-1300 3.9MII) 2.0( -21) 1.0(-20) 1.0(-17) 1300-lUOO 2.72(10) 2.0( -21) 7.2(-l8) 1.0(-17) 1^00-1500 7.30(10) 2.0( -21) 1.3(-17) 3.6(-l8) 1500-1600 2.5M11) 0.0 8.8(-l8) 1.6(-18) 1600-1700 1.03(12) 0.0 3.3(-l8) 8.3(-19) 1700-1800 1.59(12) 0.0 5.3(-19) 6.0(-19) 1800-1900 3.7^12) 0.0 1.8(-20) ^.5(-l9) 1900-2000 8.88(12) 0.0 2.8(-23) 3.2(-19) 2000-2100 I.6M13) 0.0 1.6(-23) 5.0(-l9) 2100-2200 3.58(13) 0.0 l.0(-23) l.0(-l8) 2200-2300 5.68(13) 0.0 6.6(-2U) 2.3(-l8) 23OO-2UOO 7.10(13) 0.0 l+.7(-2i+) 5.0(-l8) 2^+00-2500 9.90(13) 0.0 l.5(-2l+) 9.0(-l8) 2500-2600 1.^5(1^^) 0.0 5.0(-25) 1.0(-17) 2600-2700 2.21(11+) 0.0 0.0 5.0(-l8) 2700-2800 3.12(11+) 0.0 0.0 2.3(-l8) 2800-2900 6.12(11+) 0.0 0.0 1.3(-l8) 2900-3000 9.38(11+) 0.0 0.0 7.0(-19) 3000-3100 1.09(15) 0.0 0.0 2.0(-19) 3100-3200 1.27(15) 0.0 0.0 8.0(-20) 3200-3300 1.1+9(15) 0.0 0.0 2.0(-20) 3300-3'+00 1.83(15) 0.0 0.0 1.0(-21) 3^00-3500 2.05(15) 0.0 0.0 1.0(-21) 3500-3600 2.12(15) 0.0 0.0 1.0(-21) 3600-3700 2.23(15) 0.0 0.0 0.0 3700-3800 2.1+0(15) 0.0 0.0 0.0 3800-3900 2.36(15) 0.0 0.0 0.0 3900-i+OOO 3.06(15) 0.0 0.0 0.0 -60- Table 1. Internal Data (Cont'd ) Wavelength Solar Photon Flux Absorption Coefficients (cm /lOOA) Band(A) Outside the Atmosphere Nq O2 3 (photons/cm2-sec-100A ) 2 _) ilOOO-i^lOO 3.53(15) 0.0 0.0 0.0 !+100-i+200 3.93(15) 0.0 0.0 0.0 i+200-U300 ^.25(15) 0.0 0.0 0.0 U300-i+^00 ^.^5(15) 0.0 0.0 0.0 l4.UOO-i+500 ^.73(15) 0.0 0.0 1.0 [-2k) U5OO-U6OO ^.9M15) 0.0 0.0 1.3 :-22) U60O-U7OO 5.00(15) 0.0 0.0 2.0 :-22) U700-lf800 5.20(15) 0.0 0.0 3.3 '-22) i+800-1^900 5.20(15) 0.0 0.0 5.0 :-22) ^1900-5000 5.20(15) 0.0 0.0 8.0 [-22) 5000-5100 5.20(15) 0.0 0.0 1.1 .-21) 5100-5200 5.20(15) 0.0 0.0 1.3 :-2i) 5200-5300 5.30(15) 0.0 0.0 1.5 -21) 5300-5^00 5.30(15) 0.0 0.0 1.8 :-2i) 5^00-5500 5.^0(15) 0.0 0.0 2.1 -21) 5500-5600 5.^0(15) 0.0 0.0 2.5 -21) 5600-5700 5.^0(15) 0.0 0.0 2.9 -21) 5700-5800 5.^0(15) 0.0 0.0 3.M -21) 5800-5900 5.^0(15) 0.0 0.0 k.0( -21) 5900-6000 5.^0(15) 0.0 0.0 ^.7( -21) 6000-6100 5.^0(15) 0.0 0.0 J+.8( -21) 6100-6200 5.^0(15) 0.0 0.0 k.0{ -21) 6200-6300 5.^0(15) 0.0 0.0 3.3( -21) 6300-6^00 5.35(15) 0.0 0.0 2.8( -21) 6^00-6500 5.35(15) 0.0 0.0 2.M -21) 6500-6600 5.30(15) 0.0 0.0 2.0( -21) 6600-6700 5.25(15) 0.0 0.0 1.7( -21) 6700-6800 5.20(15) 0.0 0.0 1.5( -21) 6800-6900 5.15(15) 0.0 0.0 1.2( -21) 6900-7000 5.10(15) 0.0 0.0 1.1( -21) 7000-7100 5.05(15) 0.0 0.0 7.8( -22) 7100-7200 5.00(15) 0.0 0.0 5.0( -22) 7200-7300 ^.95(15) 0.0 0.0 3.2( -22) 7300-7^00 ^.90(15) 0.0 0.0 2.0( -22) 7^00-7500 ^.85(15) 0.0 0.0 0.0 7500-7600 1+. 80(15) 0.0 0.0 0.0 7600-7700 ^.72(15) 0.0 0.0 0.0 7700-7800 i]-.6ii(l5) 0.0 0.0 0.0 7800-7900 ^.56(15) 0.0 0.0 0.0 7900-8000 U.ii8(l5) 0.0 0.0 0.0 8000-8100 ^Alf(l5) 0.0 0.0 0.0 8100-8200 ^.39(15) 0.0 0.0 0.0 8200-8300 U. 3^(15) 0.0 0.0 0.0 8300-8i+00 ^.30(15) 0.0 0.0 0.0 -61- Table 1. Internal Data (Cont'd' Wavelength Solar Photon Flux Absorption Coefficients (Gill's/loo A) Bandd) Outside the Atmosphere ^2 °2 °3 (photons/cm^-sec-lOOA) c. C- 8ifOO-8500 ^.25(15) 0.0 0.0 0.0 8500-8600 ^.21(15) 0.0 0.0 0.0 8600-8700 ^.17(15) 0.0 0.0 0.0 8700-8800 ^.13(15) 0.0 0.0 0.0 8800-8900 U. 09(15) 0.0 0.0 0.0 8900-9000 ^.05(15) 0.0 0.0 0.0 9000-9100 1+. 00(15) 0.0 0.0 0.0 9100-9200 3.95(15) 0.0 0.0 0.0 9200-9300 3.90(15) 0.0 0.0 0.0 9300-9U00 3.85(15) 0.0 0.0 0.0 9J+00-9500 3.80(15) 0.0 0.0 0.0 9500-9600 3.77(15) 0.0 0.0 0.0 9600-9700 3.7M15) 0.0 0.0 0.0 9700-9800 3.71(15) 0.0 0.0 0.0 9800-9900 3.68(15) 0.0 0.0 0.0 9900-10000 3.65(15) 0.0 0.0 0.0 10000-10100 3.62(15) 0.0 0,0 0.0 10100-10200 3.58(15) 0.0 0.0 0.0 10200-10300 3.5M15) 0.0 0.0 0.0 10300-10^00 3.50(15) 0.0 0.0 0.0 10^00-10500 3.^6(15) 0.0 0.0 0.0 10500-10600 3.^2(15) 0.0 0.0 0.0 10600-10700 3.39(15) 0.0 0.0 0.0 10700-10800 3.36(15) 0.0 0.0 0.0 10800-10900 3.33(15) 0.0 0.0 0.0 10900-11000 3.30(15) 0.0 0.0 0.0 11000-11100 3.26(15) 0.0 0.0 0.0 11100-1 1 200 3.22(15) 0.0 0.0 0.0 11200-11300 3.18(15) 0.0 0.0 0.0 11300-11^00 3.1^(15) 0.0 0.0 0.0 lluoo-11500 3.10(15) 0.0 0.0 0.0 11500-11600 3.06(15) 0.0 0.0 0.0 11600-11700 3.01(15) 0.0 0.0 0.0 11700-11800 2.97(15) 0.0 0.0 0.0 11800-11900 2.93(15) 0.0 0.0 0.0 11900-12000 2.88(15) 0.0 0.0 0.0 after Allen,, I965. -62- 5 . REFERENCES Allen, C. W. (1963), Astrophysical Quantities, 2nd ed., University of London, The Athlone Press. Hunt, B. G. (1966), Photochemistry of ozone in a moist atmosphere, J. Geophys. Res., 71, 1385-1398- Kenesha, T. J. (I962), A computer program for solving the reaction rate equations in the E ionospheric region, AECRL Res. Rpt . 62-828. Kenesha, T. J. (I963), A solution to the reaction rate equations in the atmosphere helow I50 kilometers, AFCRL Res. Rpt. 63-7II. Kenesha, T. J. (I967), A technique for solving the general reaction- rate equations in the atmosphere, AECRL Res. Rpt. 67-0221. Romanelli, J. J. (196^), Runge-Kutta methods for the solution of ordinary differential equations, in Mathematical Methods for Digital Computers, A Ralston and H. S. Wilf, editors, John Wiley and Sons, New York. Swider, W. , Jr. (196^), The determination of the optical depth at large solar zenith distances. Planet. Space Sci., 12, 76I-782. -63- APPENDIX A. SAMPLE INPUT DECK THIS IS THE UV. W RUN FIRST 06 03 05 06 04 01 06 06 05 07 04 01 07 05 01 08 04 01 09 06 01 09 04 01 09 09 01 09 10 01 08 10 01 08 10 01 06 10 01 08 03 05 08 09 05 08 08 05 10 10 01 09 13 01 06 13 01 08 13 01 07 12 01 07 11 01 06 09 05 10 04 01 99 03 01 01 018 025 2.65E-19 3.75E-25 03 01 01 012 018 3.40E-19 03 01 01 001 Oil 6.00E-19 1.50E-17 04 01 01 032 074 8.00E-20 0,00E+00 1.30E-22 1 .50E-21 4.65F-21 1.70E-21 2.05E-22 04 01 01 001 031 1.5nE-20 l.OOE-18 3.60E-18 1 .OOE-18 1.30E-18 11 01 01 002 019 l.OOE-18 l.OOE-17 2. OOE-18 13 01 01 HUNTS ITH THE TO 5TE 04 05 03 03 03 05 03 03 06 05 09 03 08 03 10 03 11 06 11 03 12 03 09 09 09 03 10 05 11 05 12 05 13 03 11 10 09 10 09 09 10 05 09 03 REACT DEAC ADY 5 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 03 1.80E-20 06 07 01 l.OOE-20 07 07 01 6. OOE-18 9. OOE-18 03 06 01 2.00E-20 O.OOE+00 2.00E-22 1.80E-21 4,e0E-21 1.50E-21 07 03 01 3.00E-20 2. OOE-18 1.60E-18 2.30E-18 7.00E-19 09 08 01 1 .OOE-17 l.OOE-17 4. OOE-18 09 09 01 ION SE TIVATI TATE. 8. 5. 2. 1. 1. 2. 5. 5. 2. 1. 2. 1. 1. 7, 2. 2. 3. 4. 1. 1. 1. 1. 1. 1. T, CHANGED ON OF OlD AND THEN R OOE-35 3. 60E-11 70E-33 OOE-11 OOE-10 60E-11 OOE-11 OOE-13 80E-12 OOE-11 OOE-13 OOE-11 OOE-11 40E-32 50E-31 60E-32 OOE-12 OOE-13 OOE-15 OOE-13 OOE-11 OOE-11 40E-31 OCE-14 TO EXTEND THE SOLAR INCREASED TO 1.0E-10» ESTART WITH JOSHUA AT OOE+02 +0.OOE+O0 +4 OOE+02 +0.00E+O0 -2 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0,00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0,00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0,OOE+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +O.OOE+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0,00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0,OOE+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 OOE+02 +0.00E+00 +0 FLUX DOWN THROUGH AND EQUIL l.E-8 AT LAT 90.0 DEG .45E+02 .85E+03 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 ,OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OOE+00 .OCE+00 .OOE+00 .OOE+00 2.80E-23 1.60E-23 l.OOE-23 6.60E-24 4.60E-24 7.20E-18 1.30E-17 8.80E-18 3.30E-18 2.65E-19 1.40E-17 3.80E-18 l.OOE-21 0. OOE+00 3.30E-22 2.10E-21 4.00E-21 1.25E-21 5.00E-20 3.60E-18 8.30E-19 5.00E-18 2.00E-19 2.00E-17 8. OOE-18 2. OOE-18 2. OOE-17 l.lOE-18 l.OOE-21 O.OOE+00 5.00E-22 2.50E-21 3.25E-21 1.08E-21 9.00E-20 7. OOE-18 6.00E-19 9. OOE-18 1.30E-17 5. OOE-18 2.00E-19 2.30E-17 2.30E-17 2.10E-17 l.OOE-21 0. OOE+00 8.00E-22 2.90E-21 2.e0E-21 7.80E-22 1.70E-19 l.OOE-17 4.50E-19 l.OOE-17 1.50E-17 1.40E-17 0. OOE+00 0. OOE+00 1.13E-21 3.35E-21 2.40E-21 5.00E-22 3.20E-19 1. OOE-17 3.20E-19 5. OOE-18 1.80E-17 3. OOE-18 O.OOE+00 l.OOE-24 1.30E-21 4.00E-21 2.00E-21 3.20E-22 6.00E-19 1. OOE-17 5.00E-19 2.30E-18 1.20E-17 7.00E-19 -6k- 002 030 l.OOE-18 l.OOE-17 2.00E-18 2.60E-19 1.60E-20 99 00 1.000 2.12E+02 2.28E+02 03 2 3 1 1 l.OOE-17 2.00E-17 1.30E-17 1.50E-17 l.OOE-17 8.00E-18 5.00E-18 l.'t0E-17 4.00E-18 2.00E-18 1.40E-18 7.00E-19 1.80E-19 1.50E-19 l.OOE-19 5.OOE-20 1.80E-17 1.20E-17 3.00E-18 7.00E-19 4.70E-19 3.50E-19 3.00E-20 2.00E-20 E+02 3.000E+01 l.OOOE+01 1.81E+02 1.85E+02 2.17E+02 1 1 1 2 O.OOE+00 +9,OOE+01 3 3.64+012 7.66+016 2.58E+02 2.70E+02 2.48E+02 1.80E+03 1.30E+06 2.59E+05 3.30E+07 +0.00E+00 09 22 12 00 00 1.32+013 0.00+000 8.07+013 3.68+014 0.00+000 0.00+000 27+015 00+000 l.OOE-08 4.26+015 1.69+016 0.00+000 0.00+000 8.21+006 2.89+012 5 9.98+012 3.83+017 6 3.34+012 7.50+007 7 1 .11+002 4.35+000 3.57+007 1.74+008 5.04+008 1. 0.00+000 0.00+000 0.00+000 0. 6.48+013 4.03+014 1.84+015 6 0.03+000 0.00+000 0.00+000 8.69+011 5.58+010 8.47+009 2 0.00+000 0.00+000 0.00+000 4.12+001 0.00+000 2.92+001 0.00+000 1.71+001 1. 0.00+000 0. 26+010 00+000 .33+015 .00+000 .22+010 .00+000 16+002 00+000 1.17+011 1.41+012 0.00+000 0.00+000 2.13+016 8.46+016 0.00+000 o.oo+ooo 1.10+010 1.62+009 0.00+000 0.00+000 1.70+002 6.37+001 0.00+000 0.00+000 9.95+007 1 .15+000 9 1. 20+002 6.88+005 10 8.00+000 3.35+007 11 1.50-002 1.51+012 12 4.30+005 8.72+005 13 5.50-008 6.91+009 99 4.75+008 0.00+000 1. 00+004 0.00+000 3.49+003 0.00+Oon 1.71+002 0.00+000 1.99+008 0.00+000 3^56-002 0.00+000 9.07+007 0.00+000 2.26+005 0.00+oon 3.92+005 0.00+000 7.35+005 O.OO+OOO 3.11+009 0.00+000 2.41+003 0.00+000 7.92+006 0.00+000 1.19+006 0.00+000 4.71+006 0.00+000 6.67+009 0.00+000 4.45+009 0.00+000 4.94+005 0.00+000 1.58+006 0.00+000 1.32+006 0.00+000 4.29+006 0.00+000 3.07+010 0.00+000 .88+008 .00 + 000 3.87+005 0.00+000 1.83+005 0.00+000 3.24+006 0.00+000 1.27+007 0.00+000 1.52+011 0.00+000 1.88+008 0.00+000 5.67+006 0.00+000 1.07+003 0.00+000 1.89+006 0.00+000 4.72H O.OOh 007 000 8.05+011 0.00+000 5.26+007 0.00+000 5.97+008 0.00+000 -65- APPENDIX B. SDCBIT OUTPUT FOR THE INPUT SHOWN IN APPENDIX A. THIS IS HUNTS REACTION SET» ChANgED TO EXTEND ThE SOlAR FlUX DOWnj THROUGH THE UV. WITH THE DEACTIVATION OF OlQ INCREASED TO l.oE-lOi AND EQUlL l.E-8 RUN First to steady state* and then restart with joshua at o at lat as.o deg THERE HAVE REEN 24 REACTIONS INPUT INVOLVING U SpECjES. THE RFiCTIONS AND THEIR RATE COEFFICIENTS ARE 1 6* 3+ 5 4* 5* 1 K = 8.00-035«(TEMP/3.00*002)«» 0.00»EXpF( 4,45*002/TEMP) 2 f* 4+ 1 3* 3* 1 K = 5.60-011«(TEMP/3.00*002)»» . 00»EXpF ( -2.854.003/TEMP) 3 h* 6* 5 3* 5* 1 K = 2.70-033*(TEMP/3.00*002)»» 0.00»EXpF( . O0»000/TEMP) 4 7* 4+ 1 3+3+1 K = 1.00-011«(TEMP/3.00*002)«* 0.00»EXpF( . 00 + OOO/TEMP) 5 7+ 5* 6 8* 4- 7 q* 6-t 8 9+ 4- 9 9* 9'i 10 9*10^ 11 8*10^ 12 n+10^ 6+ 5+ 1 K = 1.00-010*(TEMP/3.00+002)*» 0.00»EXPF( . OO+OOO/TEMP ) 9* 3* 1 K = 2.60-011*(TEMP/3.00*002)»» 0.0O»EXpF( . 00*000/TEMP) 8* 3* 1 K = 5.00-011*(TEMP/3.00*002)*» 0.00»EXpF( , 00*000/TEMP) — ID* 3* 1 K = 5.00-013*(TEMP/3.00*002)«* 0,00*EXpF< . 00*000/TEMP) --11+ 6+1 K = 2.80-012*{TEMP/3.00*002)«» 0.00»EXpF( . 00 + OOO/TEMP) 11* 3+1 K = 1.00-011«{TEMP/3.00*002)«» 0.00»EXPF( , 00*000/TEMP) 12* 3* 1 K = 2.00-013»(TEMP/3.00*n02)«» O.00»EXpF( . OO + OOO/TEMP) o _. ^n o -H n o ^ n o — Kn o -^fl o -•no -Hn o — r) n ro o ♦ no* n o ■» f>-) O ■)■ mo* no* no* no* in ivj — -< o 1 o ^ O 1 o .-• O 1 o »^ C 1 o ^ O 1 o ^ o 1 o —< O 1 o -^ O 1 o -.-- z »H 1 o o ^ 1 o o ^ 1 oo --« » o o ^ 1 oo ^ 1 o o ^ 1 o o -^ 1 o o ~ ~ ^ a. o o o o o o o o o o o o o o o o o o o o O O O C3 r«- X a lij 1 o^ » 1 Os*^ • 1 o <^ • > o-* • 1 Ov* « 5 0-* • 1 O •* o 1 o -* • — a u. »- O IT . O o in • o o in so O If! • O o in • o o in • o O IT . O z uj »- O • rt o » -^ O o -< O o ^ O • r-1 O • -H a t~ C (VI C (V O CVJ o rvj o rvj O (V c w O (V — O — (\j n -* (Vl — -t - s: r q: « z a ui — CE IjJ I— i: uj f- a. h- lU tvi ^ o no* ^ o I o -H I O O o o o o I o o • o -♦ • o -< I o o o o o o I O O o o -* • o o o o o loo* O 4^ e O ^ S o o o o o o ( O O o O ■* • o CVl rf o no* IT-) O t O ^ I O O O O O O I O O o O <»• a O O e ^ %0 (^ (VJ -I O no* ^ o 1 o ^ I o o o o o o too » o .* • o (VJ •-< O no* ^ O I o .-1 I o o o o o 1 o o • o -* • o o • -^ tVJ '-^ o no* ^ I o o o o o o loo* o ~* • o — a — 1 •-< o -1 <-■ <=> -. -< o rl— 1 O r-. ,-. O ^-1 o .-■ -I o ^ ^ o ^ — (V ^ o * — o * ^ c + — o -^ .- o * »- c * ^ c * ^ o * n w ~ O O 1 o O O 1 o O O o o C 1 o O O I o o o t o o o 1 o o O 1 o -• — s: r- 1 O O ^ 1 o o r-l 1 O O —> 1 o o p-t 1 o o ^ 1 o o ^ 1 o o ~ — 3: QC O O O O o o o o O O O o o o o o O O OO o o o o o o o o o o o o in 3: a iij 1 OO . 1 o o • BOO. 1 OO o 1 O O e 1 o o • loo. loo* — a iij »- O O • o O O • o o o • o O O • O o o • o O Ci • o o o • o s: u h- o • ^ O • r-< o • ^ O • Pi O e —! o . ^ o • -^ O • .-1 a }- O =-t O — O —1 O •-! o ^ c -* o ^ — oc —. — , o ^ -. c ^ .- o .- ^ o — 1 -< O -1 -< c — 1 -1 o .-> ^ o -~ O fVl ^o * ^ o * «-• o * —to* -40 ♦ r-. O * -fO * OJ fU — r-. O 1 O ^ O 1 o .-» O 5 O -< O 1 o r-< O 1 O — 1 O I o -i O 1 o -< O 1 O — - z: -^ 1 O o ^ t o o ^ 1 o o -« 1 o o rH I O O r-l 1 O O -< 1 o o ^ 1 o o - - i: a: o o o <:> o o o o O O O o -« s: q: L. loos 1 OO • 1 O O o 1 OO o 1 o o » 1 OO e loo. 1 oo . " CE liJ »- o o * o O O • O O O o O O O o O o o • o O O o O o o » o o o * o 2; lij h- O • .-4 o • ^ O O -4 O • -1 o • — • o • -• o • -• „ in o in o in o in o in o in o in o in o — t-- n -. o n -< o n — ■ o n -« c n ri o n r-l o n — c ~ £> M ^o * -. o * r-l C * -• O * -1 o * -. o * — c * -< o * nolo nolo n o 6 o nolo nolo nolo nolo nolo ^ — s: n 1 o o n 1 o o n 1 o o n J o o n B o o n » o o n 1 oo n 1 o o s a c o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o n z auj I o o •» 1 O O o [ O O o too. 1 oo . loo. — a u. 1- o o . o o o . o O O o o O O o o c o . o o o » o . ^ o — o 1 ^010 in 1 >o 1 >o 1 -1-2: -.100 -.100 -.100 ^100 -. ) 00 — — s: a 0000 0000 0000 0000 0000 0000 0000 0000 fvj i: a uj 100. 1 00 * 1 • J 00 » too. loo. 100. loo* - a iij ^- -* . in • CVI . C3 r^ « 00 • cr . o^ . f^ . 2: IjJ H- (VI • ^ tVi .sx -* . -* . <* (VJ . St in . sf -. .s* 00 . -< 1 — 1- (VJ (VJ CVI (VJ (VI (VJ CVJ CVI 1 l±J (/) 10 < in (VI -< (VI — 1 (VJ -^ (VJ r^ (VI —1 (VI -^ CVJ ^ CVI ^ 1 K to (Z -J f^ (V -^ c — — . c ♦ — < * -. * ^ * ^ ♦ 1 < z C 2 ,^ r-< <» 1 -* -* 1 -*^ 1 -* 1 sf 1 -* I •* 1 1 CC z 0^ X n 1 n 1 n 1 n 1 n 1 n 1 n 1 n 1 u to „ 2; a 000 C3 0000 0000 c 0000 0000 1 1/1 Ul U UJ Ijj —1 2: ct 1 1 . 1 • 1 » loo* 1 CD » loo. 1 • 1 z (/) llJ to UJ a UJ r- CD • 00 • «D CO 00 a> 00 . c^ CO . 00 CD . (VJ CO . n CD .0 1 —> * 10 q: 2: UJ t- CM . n in » n vo * n -ton 00 on in . n -• • n t^ ' r, in UJ UJ z>z -* 1- >c 0^ CO -£) kn (vj — >*• n n n n n n nm n n n n c5 n n n n ^ 5 n >c (V r-t rt 1 -H 1 -. 1 -• 1 ^ c 1 r-I 1 -• 1 ou in UJ „ 001 1 0010 0010 1 0010 0010 0010 20 to c >o M CC — 2: 1 1 0100 0100 1 0100 0100 0100 «_i Q: 2: tr 000 0000 0000 c 0000 c 0000 t-CE 13 (/) 2: CC UJ * >o «c e * vO * .o * + s: z • z a -i cr CO 1^ >o in -* n 1 iij UJ UJ Ul 1 a: _i q: M a I iLr^ -69- fVJ o ^ o in o >o o «o o <0 O c o in o n o o o o o o o o o o o o o o o o o o o o ♦ o •♦ o ♦ o ♦ o ♦ o ♦ ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o o o o o o o o o in o in o r^ o (VI o oo o m o r- o in o (J. . r^ • o . c^ . in . a) • o . • o • o • o • o • o • o • o CT> ■* o 1^ •-• —> — • — • (V o .-i c — c ^ o (V o (V' c ^ o o o o o c o o o o o o o o o o o o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ « o ■» o ♦ o ♦ o o o o o o o o o o o o o o o o o ^ o fVI o W o ^ o ^ o o o t^ o in o • o • o . o . o • o • o n • • o « o o o o o a o o o o o OS o r^ o o o -^ o ^ o oo o o o ♦ o ♦ o ■» o ♦ ♦ o ♦ o ♦ o ♦ o ♦ o o o o o o o oo o o r- o (VJ O o o (VJ o o o (VJ ^ ^ —1 o o OsD O rt o ^ O 1^ O 1^ •~- oc ir OD if m ■* ■* o- m (t ri 'C (T . <* . 00 . IT, • in • m • (VJ • •* • m • in in OJ (T (^ oc CD r- in ■O O h- o 00 O <0 O o o ^ o (VI o (VI o O O o O o o o o —1 o -H o ~ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o + o * (VJ * o •♦ o ♦ o ♦ o ♦ o ^ O O O 0> O -1 o in o oo O CO O -C O (VI — ^ ro r- o >* -H s* V* >0 CO f>- OD -4 (VI (y r^ > (VI . in • f^ • o . (VI . -* • OD . in -H (^ o o o o ^ V* O IT -^ — o I o ♦ o ♦ ♦ O +0 or- o r^ o) O (VI O r-« O * o (^ O 1^ O (VI O (VJ o . o » O • O . O . o • o • l/l (VJ O (B o r) o n O .* o .* O ^ o ^ o (n z - UJ >- "1 •-> .-• —t --) —■ ^H —1 -70- „ stl^ fVJ .-* -I 0. ^ ,-, CD in n n r) r) in oc — n r5 1 1 c- 1 ♦ ♦ 00 ♦ ♦ « in (vj — -* 1 ^ in 1 a in 1 >o in •♦ ^ ^ « r) * ♦ ro in ♦ n n ♦ « t^ r- 1^ r^ ^ in ^ r>- r) (^ >C -< - s: a u ♦ 00. ♦ >D <0 • « 1^ ^ • ■» >o >o • ♦ ^.C . ♦ St . « in XI . ♦ o> h- . - a li. 1- — . r^ r- -* • ^ r^ • r ^ r . * St I' • -* -* — ' • V* «»< r\j . in r- >» . 2: UJ »- in • >* CD • 00 m • 1- rvj "fyj h- . tu 00 • ^ (VJ • •X) • (\J _ -. m -* ^ (VJ r) (VI in --• — OJ ^ -* in 00 in in - ♦ ^ ♦ r- in ♦ in in ♦ St r^ 00 (^ cr «*) r- in CD a (VJ (VJ (VJ (VJ r^ «c z a uj ♦ r- <£- • ■» ^ . ♦ CC -H . « (VI rt • •♦ in * • « (VI in . ■» — . « (S in . — q: iij »- fM >0 • O) CVl • (D n -H . ^ . .* r- St- . ,£. (VI . -. ^ . ,0 o> -* . xiij »- tvi . ^ CD • f^ in .(VJ a r- • n r^ .(VJ .in B ..0 a >- ^ (V s» ^f\i -* (VJ a in — < cr in (VJ u >- <\J -*• ^ .-■ in in m CD -* in -c (VI St in in in >o in St — ^ (V c c ■» c ♦ C c « c c ♦ c c ♦ c ♦ ♦ c ♦ n (VI — in 1 vo in 1 o 1 in »^ 1 (Vi 00 1 00 • n (E 1 r- rt — s: ♦ 05 (VJ ■» in (VI ♦ (VJ ♦ cr St •♦ (VJ (n ♦ r>- - in -^ --•r- ro 0^ (VI -» 00 - • orr, » ♦ h- (VI . (vicr . ♦ Oh-* 0(VI . (VJ «^ • 0- r- . ~ tr uj *- -* «c • in • (VI <^ ^ . ^ -£ . m (VJ in • r^ h- >* • ro >o . in r^ -* . r-l z iij 1- • -* •-• • CD c> .cr (^ • f- (o .in (U .(n CD .(n ^ 1 (O (VJ 1 (VJ (VI 1 r- (VJ 1 >-. — 3: 1 s* f^ ■» .C St ■» f^^- ■♦ t^ rt ♦ (^ ^ ■» (VJ ^0 1 -t OJ 1 r- — — s: q: <»• 00 ^ .-« (Vi in o (VI in St CD rvj (VI 00 r- -* ^ X tt liJ !<£>». 1 -* >c . 1 (O ^ • 1 (VI {^ . ^0 . (O c . ♦ —. a . - q: u •- r- (7> .(H 00 «£> . ^ soin .t^ i>- . <\j cr h- .in 00 . ^ CO . in CD ct) . ^ 3: UJ »- t- .in -* . (VJ o (VI (T in (VJ ro iij CJ> —1 in c (VJ •» OlD ■♦ sD ro ■» r- 1 in 1 (VI 1 St ~ - X a St ~* 1^ (^ r- (VJ -t en (VJ IT (^ o n • ♦ incvj . ♦ in (»> (VJ • 00. ♦ (VJ^ . — a ii. 1- <£ ir> . IT -H (V • (^ ^ . ^ 1^ sj • (^ ro (^. . X tf • ir 0.0 ff ..-. (VI .>* in •<» CO .vO 00 »Ch (J- . in a 1- '1 (O (^ CO r^ 1^ fs- in — cr -* cr f^ r- iij n ,-1 n ro CD sD in in -* -* -t ro in -* _■ 00 c s* ^ c rt r. (VI in rr, - vD CD -^O ^ ^ rt -< ^00 (VI (VJ ^ (VJ — (» (VJ ♦ ♦ ♦ ♦ 1 00 1 1 ♦ -H — n 1 n (VJ 1 (7- (\j 1 in (V) 1 (O in 1 ro •♦ in -* ♦ ^ - 3: 1 CD 1 00 r-- 1 —• + cr ^ ♦ -• CJ> ♦ (^ IT ♦ fvi ty ♦ fVJ 3: a -H r5 CD in a (VJ (VJ z: ct uj ♦ (^ CO . ♦ in . « st r^ . * tn • ♦ (\j . -•(^ . sO St . ♦ CD (X) . - CC IJ ►- >o . in fo CD . ro in CD .r^ %c . -H CD ^ . ro in <» . c> in (VJ . 3; UJ •- (VJ . (VJ ^ • r- (\J • (VI -• . (VJ in • (VJ CD • r^ . >* — »- in ^ ^H ^ (VI (^ tn in I/) < — in CD r^ -♦ -H 000 ^00 000 Dl/1 CD _l r- (V, c ♦ ♦ ♦ ZQ 2 ~ -H ^- -t 1 en in 1 ro vO 1 in ■^ ♦ fs- r- ♦ (VI CD ♦ St cr « St cr •♦ ro oz . cr - 3: 1 ^ vO 1 in 1 a) vO ♦ 00 in ♦ in — ♦ cr St ♦ (VJ ♦ cr r^ uo (J to „ - 3: cc cr (VJ vO ro CO t^ ^ in r- in * ^ ro • l/)UJ I/) UJ a UJ i- cr * . (VI ro . in St .in cr . in ro (D . -« 00 cr . St St cr (VI m • (T .in oc 3: Ul »- in • (O —t « -* -» . «o r^ . in in . -* ro . ro if) (^ (5 q: t- cr St cr sO •c st (VJ cr (VJ UJO UJ UJ . • »-(t) CD 1- r~ vC ^^ fO oc -* —I _■ 3— ♦ 1 rs- (V ~t IT in z . r^ 1 -^ ro .-1 c ru ro St in Ct- St •-■00 (VI 1 — -* -* ^ ^00 (VJ ro in yO M tt — 3 1 cr 1 V) 3: tr UJ ♦ CO in . > cr r^ . « cr -c . « CO > ■» 00 . * (VJ in CD ro . 11 H M in a Ul 1- (VI • ^ cr • cr . -o • (VI -» . in cr . — n -C oD . cr IliJ 3: a UJ t- « rvj » • CO . m • -t * r-l (VJ ~-3 UJ UJ cr Ul t- IT. —^ —1 (O CD ^ .-> cr >- 3: _J Ul t- (VJ (VJ OJ (VJ •c UJI- ►HO (/) . • 3: 1-2 en n >o •-tr o in ■* 00 ID UJ Ul UJ .-• _IO £r M q: I O -* O IP o i o o o o o ♦ O ♦ C3 ♦ O ♦ O +0 o .-• o o o o o o >r o • ,-« • fVI • a o ■-< o t^ o -o o >0 o o o o o o o o ♦ o ♦ o ♦ o o o o o ♦ --< o o o o o o o o o o o o o o o o o o o ♦ o ■» o ♦ o ♦ o ♦ o ♦ o ♦ o * in ■» o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o ♦ o o t^ o - o 0^ o r^ o o o •C O in o o- • tn . r- • o . o • o • o • o • o • o • o ^o ^o ^o oo ~ o ■* o ■» o* o ■* -* ♦O ♦O ♦O +0 .-< oo oo r-o -♦o — •♦O O'O CDO oo o o o o Co -HO .-1 c o o O •♦ O •» C3 ■♦ ♦ O +0 ♦ o (CO no CD O P) o O O .-I o (VI • — I • o ♦ O) O -H O 0 O ^ (T IT a. 4 n n -* a n ff r xf * a o> . ^ . o • 00 • f^ • n '- o in 1^ oo in CD cr .o n o >o -• in St h- 00 r- 00 O OJ >- fVJ • >c • I-- • (VJ • (\J • s* • -* • -« • 'I 'in (VI (VI no -* o o I ■♦ o o o -* in «o • n o ♦ o ♦ o ♦ O ♦ o ♦ ♦ o ♦ o ♦ O o r- O 1^ O (VI o in o .-• O) vO r^ o so in (^ o >o in o (J> .c 0> o o> (VJ I- — (T o o ^ oeo a <£> cr (V (r (Vi (^ h- a n o • o • <^ • (T • cr • a> • (r • (/) (VJ (T CO o n on tr •♦ (> •* o> r^ (T ^ <^ n • • • • • • O V ^ V, r-l •■ O I • I •• o -o •--. •a »c -Z. ^HO r-in r-«-H ^HO IV I \ IX I "v ■ V •v IV. I \ o o o o tM o • r-H » O o - •■ o • o • o o » • o o • ► o o • o • • o O O o o o o o o o o o o o o o o • O » O » O • o o » o • • o o - o - - o o • • o 2 ^o o> o >0 O o o o o 00 o (> o C> O o V (VJ v O >v O X o \ (Vl X (\J X (VI X \ -> \ •-< \ -. \ -^ X -• V -< X ^ X '-^ Cf rvi vO tvj o tvj Cf tVJ (y (\j (J> (VI (D (VJ (O tVJ M IM (Vl (VI CVI ru (VI .-•o ino if)o ino oo ox (vjx (vjx (vjx ox XO XO XO XO X(VI inm — I r) —im .-im .-i rvj O O O O O o ox OX OX X OJ X (VJ X (VI -. (VI r-1 (VJ -^M cr oo o o oo oo oo CO oo oo o»»-«-»»»»-» 2^ .o »o .o .o »o »o -o -o O CO CO oo CO oo CO CO oo l-l(VJ .».»»..>>•• »-0 »0 "O "O •o •o •o ••.* »o tj»Z oo oo oo oo oo oo o(Vj oo DC XXXXXXXX az xn xfo x(o xi»i x(^ xn X(^ x-* I- O^ Ort Or-1 O^ O^ O •-• Ort o(Vj 1/5 LJ c t-OOOCOOOO ai •••...•• oo OOOOOOOO -)>-" o(^eor^^Olf)^*(n c no r) o o o o o o -^ O —t O » (VJ vO O ro o o o o o o o • in - o^ •a -o ".o •o •■00 •o o>^ >co r^o h-o OO a)o o- ^ o>-»-< (VJX ox ox ox ox (VJX MX (\JX XCT- xo xo xo xo xo xo xco rocs ^r-« o^ ^^ ^'-^ o*«— I Or-* CO"-* (D^H rvj tvj tvi tvj fVJ UJ OO OO OO OO <00 OO OO OO Qro>>.»>»>* _(\J .-H •o •o -o •o •o •o •o •© <: OO CO CO OO OO OO OO OO o«* • • • • • • . » •-"O O •O •O •O •O •O •CD •O ".O ►— ZZ OO OO OO OO OO OO OO OO o xxxxxxxx :d X-* X-* X-* X-* X-* X-* xcD xoo O O •-> O-i O-^ O'-* O— • O— OO OO o oc a. •-OOOOOCOC CCI ••...... OO OO OO oooo in m APPENDIX C. SDCBIT LISTING PROGRAM SIXBIT COMMON /A/ 1 BRK(4) FLUX ( 120) FTERMW(40 ) H2HOLD(20) IGAINA( 20»20) IL0SS(4 ) rLOSSD( 20,20) 8 K0(51 ) 9 KAB503 (120) X L3(51 ) 1 N ( 5 1 ) NUSRUS(400) QRK{20) RR1(21 ) TAB02 ( 87) TfRMW( 100 ) YNEW(20) Y3SUN( 20) COMMON /B/ 2 3 A( 51 ) .CRK(4) ,FLXINF( 12 .G( 20) ,H3HOLD(20 , IGAINB( 20 ,IL05SA(20 , INFORM( 30 »KABS(20,1 ,! ) 9000 RETURN END ) ,20) ,20) ) 20) 00) 20) 4, IF , IT A, JZ , NS , TM , TS 2,R3, SECT LONG, Z( JZl 100 LAG5 CHXX END PEC AX2 EC RRl, X K5,K ) ,JZ RK(4) (20) RATE(20) G(20) ATOM( 10,20 GAINC(20,2 LOSSB(20,2 PULL (20 ) ABSN2( 120) 1(51 ) L2(21) FTERM(20) ROHLD(71,2 2(51 ) R3(21) EMP( 20) HOLD(20,20 P R I N T ( 2 , 2 (20) , HCOUNT, , IFLAG6, , I UN IT , , JZMAIN, , P I , TMIN , , ZBOTUM, RR2,RR3,RX ) 0) IFLAGl IFLAG7 JCNTRL LAT T TMONTH ZSTEP ,B(51) ,FF(20) ,FTERM( 40) ,HBAR (20) , IGAIN ( 4) , IGAIND (20,20) ,ILOSSC(20,20) ,K( 1000) ,KABS02 (120) ,L2{51 ) ,LL3(21) ,NTERM( 20) ,00LD( 20) ,R3( 51 ) ,SUMIN( 10,20) ,TERM( 100) ,YLAST( 20,20) ,YSTART (20,20) ,ZPRINT (20) IFLAG2, IFLAG8, JFLAG2, LONG , TDAY , TOLD , ZTOP 1 ,RX2,RX3' 5MAX ,KABS02 , KABS03 , K ABSN: 1=1,6) 00 = 7, JZEND) EP + 1 STEP GO TO 200 T CHECK OUT IN SUBROUTINE HEIGHTS. TH! HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE S WHE HE HE HE HE IGHTS I6HTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHTS IGHT- i SUBROUTINE COMMON /A/ 1 BRK(^) 2 FLUX(120) 3 FTERMW(40) 4 H2HOLD(20) IGAINA(20»20) IL0SS(4) ILOSSD(20,20) K0(51 ) KABS03( 120) L3(51> 1 N(51 ) 2 NUSRUS(400) 3 QRK(20) 4 RRl (21) 5 TAB02(87) 6 TERMW(IOO) 7 YNEW(20) 8 Y3SUN(2t^> COMMON / 2 3 5 6 7 LGEBRA A( 51 ) ,CRK (4) .FLXINF( 120) ,G( 20) ,H3HOLD{20 > IGAINB(20.20 ) ILOSSA(20»2C) INFORM( 30) »KABS( 20>120) »KRK( 80) .LLl ( 21 ) .NASRAS{ 1000) ,PHOLD( 20, 20) ♦ Rl ( 51 ) ,RR2(21 ) »TAB03 ( 87 ) ,Y( 20) ,YOLD( 20) »Y4SUN( 20) ALPHA , EQUIL IFLAG4 IFLAG3 ■ IFLAG9, IJZX H IFLAG5 ITCHXX JZEND NSPEC TMAX2 TSEC ,ARK(4 ) ,F ( 20) ,FRATE (20) ,GG(20) , I ATOM(10.2 , IGAINC(20» , IL0SSB(2C, , JPULL(20 ) ,KABSN2( 120 ,L1(51 ) ,LL2(21) .NFTERM( 20 ) ,PR0HLD( 71 . .R2( 51 ) ,RR3( 21 ) , TEMP ( 20 ) ,YHCLD(20,2 ,YPR INT ( 20, ,Z (20 ) , HCOUMT, IFLAG6, lUNIT , J Z MA IN, PI , TMIN , ZBOTUM, ,8(51 ) ,FF( 20) ,FTERM(40) ,HBAR ( 20 ) ) , I GAIN (4) D) ,IGAIND(20 D) ,IL0SSC(20 ,K( 1000) ,KABS02 ( 12 ,L2(51 ) ,LL3(21 ) ,NTERM( 20) 0) ,QOLD(20) ,R3( 51 ) ,SUMIN( 10, ,TERM( 100) ) ,YLAST(2C, 0) ,YSTART(20 ,ZPRINT ( 20 IFLAGl, IFLA&2 20) ,20) ) IFLAG7, JCNTRL, LAT T : TMONTH, ZSTEP . IFLAG81 JFLAG2i LONG , TDAY : TOLD , ZTOP JFLAG6, JOSHUA, NFREAC, NREAC , THOUR , TMAXl , TPRINT, TQUIT , INTEGER HC0UNT,PH0LD,R1 ,R2,R3,RRl,RR2,RR3»RXl,RX2,RX3, 1 TMONTH, TDAY, THOUR , TMIN , TSEC, TX REAL K0,N,KABS,K»KRK,LAT,L0NG,K5,K5MAX ,KABS0 2,KABS0 3,KABSN2 READ (60 ,999) ( T EMP ( JZ 1 ), JZ 1 = 1 , JZEND ) DO 10 Il=l,NREAC IXX = (I 1-1 )*JZEND DO 10 JZ2=1, JZEND IJZ = IXX + JZ2 K(IJZ) = K0( I 1 )*(TEMP( JZ2 )/A{ I 1 ) )**N( I 1 )*EXPF( B( I 1 )/TEMP( JZ2 ) ) PROHLD{ I 1, JZ2 ) = K ( I JZ) 10 CONTINUE DO 20 Ml = l, NSPEC NTERM(Ml) = 20 CONTINUE DO 30 IZERO = 1,1000 NASRAS( IZERO) = 30 CONTINUE DO 260 I2=1,NREAC LXl = LI ( 12 ) LX2 = L2( 12) LX3 = L3( 12) RXl = Rl( 12) RX2 = R2 ( 12 ) RX3 = R3( 12) IF (LXl .LT, 3) GO TO 130 IF (RXl .NE. LXl ) GO TO 110 LXl-= 1 ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA -83- 210 220 230 240 250 260 RXl = 1 GO TO 130 IF (RX2 .N LXl = 1 RX2 = 1 GO TO 130 IF (RX3 .N LXl = 1 RX3 = 1 IF (LX2 .L IF (RXl . N LX2 = 1 RXl = 1 GO TO 160 IF (RX2 .N LX2 = 1 RX2 = 1 GO TO 160 IF (RX3 .N LX2 = 1 RX3 = 1 IF (LX3 .L IF (RXl .N LX3 = 1 RXl = 1 GO TO 200 IF (RX2 .N LX3 = 1 RX2 = 1 GO TO 200 IF (RX3 .N LX3 = 1 RX3 = 1 IF (LXl .L NTERM(LX1 ) NASRA5( ( LX IF (LX2 .L NTERM(LX2) NASRAS( (LX IF (LX3 .L NTERM(LX3) NASRAS( ( LX IF (RXl .L NTERM(RX1 ) NASRAS( (RX IF (RX2 .L NTERM(RX2) NASRASI (RX IF (RX3 .L NTERM(RX3) NASRAS( (RX CONTINUE DO 310 I IZ NUSRUS( I IZ E. LXl) GO TO 120 E. LXl) GO TO 130 T. 3) GO TO 160 E. LX2) GO TO 140 E. LX2) GO TO 150 E. LX2 > GO TO 160 T. 3) GO TO 200 E. LX3) GO TO 170 LX3) GO TO 180 LX3) GO TO 200 GO TO 210 ERM(LX1 )+l NREAC+NTERM(LX1 ) GO TO 220 ERV,(LX2 )+l NREAC+NTERM( LX2 ) GO TO 230 ERM(LX3)+1 NREAC + NTERM( LX3 ) GO TO 240 ERM(RX1 ) + l NREAC + NTERM(RX1 ) GO TO 250 ERM(RX2)+i NREAC+NTERK(RX2 ) GO TO 260 ERM(RX3)+1 NREAC + NTERM(RX3 ) I2+NREAC ER0=1 ,400 ERO) = ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA -84- CONTINUE DO 320 M2=1.NSPEC NFTERM(M2) =0 CONT INUE DO 350 111=1, NFREAC 340 350 999 NFTERM(LL1 ( I NUSRUS( ( LLl ( IF (RRl (III) NFTERMIRRl { I NUSRU,S( (RRl ( IF (RR2( III) NFTERM(RR2 ( I NUSRUS( (RR2 ( IF (RR3( III) NFTERM(RR3 ( I NUSRUS( ( RR3 ( CONTINUE FORMAT (7(E8.2,2X RETURN END LT LT LT = NFTERM( LLl ( I II ) ) 3) *NFREAC+NFTERM( LL 3) GO TO 330 = NFTERM(RR1( III)) 3) *NFREAC+NFTERM{ RR 3) GO TO 340 = NFTERM(RR2 (III) ) 3)*NFREAC+NFTERM( RR 3) GO TO 350 = NFTERM(RR3 (III)) 3 ) *NFRLAC+NFTERM( RR ^- 1 Kill I 1 I 1+NFREAC Il+NFREAC I 1+NFREAC ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBRA ALGEBR- -85- SUBROUTINE START COMMON /A/ 1 BRK(4) FLUX( 120) FTERMW( 40 ) H2H0LD( 20 ) IGAINA(20,20) I LOSS (4 ) ILOSSD(20,20) K0(51 ) KABS03( 120) L3(51 ) 1 N(51) NUSRUS( 400) QRK( 20) RR1(21) TAB02 ( 87) TERMW( 100 ) YNEW(20) Y3SUN(20) COMMON /B/ , A(51 ) .CRK(4) »FLXINF( 120) .G( 20) .H3HOLD{20 ) » IGAINB(20»20) . ILOSSA(20,20) » INFORM(30) »KABS(20»120 ) .KRK ( 80) .LL1(21 ) .NASRAS( 1000) »PHOLD( 20»20) .Rl( 51) .RR2( 21 ) »TAB03(87 ) .Y(20) »YOLD(20) .Y4SUN(20) ALPHA , EOUIL • IFLAG3, IFLAG4. IFLAG9* IJZX H JFLAG6. JOSHUA, JZ NFREAC» NREAC » NS THOUR , TMAXl » TM TPRINTt TOUIT » TS INTEGER HC0UNT,PH0LD»R1 ,R2,R3, 1 TMONTH»TDAY, THOUR, TMIN,TSEC»T REAL KO,N,KABS,K , KRK, LAT, LONG, DATA ( IFLAG2=0) , ( IFLAG3 = 0) , ( IF 1 ( (ARK( I ), 1 = 1, 4)=0. 5, 0.292893, 2 ( (BRK( I) , 1 = 1 ,4)=2. 0,1.0, 1.0,2 3 ( (CRK( I ), 1 = 1 ,4)=0. 5, 0.292893 , 4 ( (HBAR( JZ) ,JZ=1 ,20) = 20(1.0) 5 ( ( JPULL( JZ) , JZ=1,20) = 20(0)) 6 ( ( (YHOLD(M,JZ) ,M=1,2) ,JZ = 1,20 7 ( ( ( YPRINT(M,JZ) ,M = 1 ,2) »JZ=1 ,2 8 ( ( (PHOLD( M,JZ) ,M=1,2) ♦JZ=1,20 9 ( ( (YHOLD(M,JZ) ,M=3»20) »JZ=1 ,2 1 ( { (SUMIN( ITCH,JZ) ,ITCH = 1 ,10) , 2 ( ( (PHOLD(M, JZ) ,M=3,20) ,JZ=1,2 DATA ( (FLXINF(LAMDA ) ,LAMDA=1,6 LAG5 CHXX END PEC AX2 EC ,ARK(4 ) ,F ( 20 ) ,FRATE (20) ,GG(20) , I ATOM ( 10,20) ,IGAINC(20,20) , ILOSSB(20,20) , JPULL (20) ,KABSN2( 120) ,L1(51 ) ,LL2(21) ,NFTERM(20) ,PROHLD(71,20) ,R2(51 ) ,RR3(21 ) ,TEMP( 20) ,YHOLD(20,20) ,YPR INT(20,20) ,Z (20) , HCOUNT, IFLAG6, lUNIT , JZMAIN, P ! , TMIN , ZBOTUM, ,B(51) ,FF( 20) ,FTERM(40) ,HBAR (20) , IGAIN( 4) ,IGAIND(20,20) ,ILOSSC(20,20) ,K(1000) ,KABS02 ( 120) ,L2( 51 ) ,LL3(21 ) ,NTERM{ 20) ,OOLD(20) ,R3(51 ) ,SUMIN ( 10,20) ,TERM( 100) ,YLAST( 20,20) ,YSTART (20,20) ,ZPRINT (20) FLAGl, IFLAG2, IFLAG8, JFLAG2, LONG , TDAY , IFLAG7 JCNTRL LAT T TMONTH, TOLD ZSTEP , ZTOP RR1,RR2,RR3,RX1,RX2»RX3, X K5.K5MAX,KABS02»KABS0 3,KABSN2 LAG6 = 0) , ( IFLAG7 = 0) ,(PI=3.1415 9265^ 1. 707107, 0.166666667), .0) , 1.707107,0.5) , ) , 3.79E08, 2.65E09, 5.88E09, 1.20E10, 7.30E10, 2.54E111 3.58E13. 5.68E13, 6.12E14, 9.38E14, 2.12E15, 2.23E15 4.25E15, 4.45E15 5.20E15, 5.20E15 5.40E15, 5.40E15, 5.40E15, 5 DATA ( (FLX INF(LAMDA) ♦LAMDA = 61 1 5.40E15, 5.40E15, 5.40E15, 5 3.80E09, 8. 2.16E10, 1. 1.03E12, 1. 7.10E13, o. 1.09E15, 1. 2.40E15, 2. 4.73E15, 4. 5.20E15, 5. 40( 1 .0) ) , = 40( 1.0) ) , 40( 1 ) ) , = 360 (0.0) ) , 1,20) = 200( 0) ) = 360 (0 ) ) 3*52E09, 1 .03E10, 3.74E12, ] .45E14, 1 .49E15, 3.06E15, 5.00E15, 5.30E15, 7.71E09, 3.94E11 , 8.88E12, 2.21E14, 1.83E15, 3.53E15, 5.20E15, 5.40E15, 5.08E09^ 2.72E101 1.6 4E131 3.12E14, 2.05E15, 3.93E15, 5.20E15, 5.40E151 ,30E15, 5.25E15. START START START START START START START START START START START START START START START START START START START START START START START START START START START START START START .START START START START START START START START START START START START START START START START START START START START START START START START -86- 2 5.20E15, 5.15E15 5.10E15, 5.05E15 » 5.00E15, 4.95E15, 4.90E15, START 3 4.85E15. 4.80E15 4.72E15, 4.64E15. 4.56E15, 4.48E15, 4.44E15, START 4 4.39E15» 4.34E15 4.30E15, 4.25E15. 4.21E15, 4.17E15, 4.13E15, START 5 4.09E15. 4.05E15 4.00E15. 3.95E15. 3.90E15, 3.85E15, 3.80E15. START 6 3.77E15. 3.74E15 3.71E15. 3.68E15 » 3.65E15, 3.62E15, 3.58E15, START 7 3.54E15, 3.50E15 3.46E15, 3.42E15. 3.39E15, 3.36E15 , 3.33E15. START 8 3.30E15, 3.26E15 3.22E15. 3.18E15. 3.14E15. 3.10E15, 3.06E15. START 9 3.01E15. 2.97E15 2.93E15, 2.88E15 ) START DATA ( (KABSN2 (LAMDA ) .LAMDA = 1,60) = START 1 5.0E-18, l.lE-17 1.2E-17, 1.4E-17, 1 .5E-17, 1.6E-17, 1.8E-17, START 2 1.9E-17, l.OE-18 1.0E-18» l.OE-ie, 2.0E-21 , 2.0E-21, 2.0E-21. START 3 2.0E-21, 45(0. OE- -00) ) START DATA ( (KABSN2(LAMDA) »LAMDA= 61.120) = START 1 60(0. OE- 00) ) START DATA ( (KABS03 ( LAMDA) ,LAMDA= 1,60) = START 1 1.5E-20. 3.0E-20 5.0E-20, 9.0E-20. 1 .7E-19. 3.2E-19, 6.0E-19. START 2 l.OE-18, 2.0E-18 3.6E-18, 7.0E-18, l.OE-17, l.OE-17. l.OE-17, START 3 3.6E-18, 1.6E-18 8.3E-19, 6.0E-19, 4.5E-19, 3.2E-19, 5.0E-19, START 4 l.OE-18, 2.3E-18 5.0E-18, 9.0E-18, 1 .OE-17, 5.0E-18, 2.3E-18, START 5 1.3E-18» 7.0E-19 2.0E-19, e.OE-20. 2.0E-20, l.OE-21, l.OE-21, START 6 l.OE-21, O.OE-00 O.OE-OO, O.OE-00 , o.oE-on , o.oE-on . O.OE-OO, START 7 O.OE-00, O.OE-00 1 .OE-24, 1.3E-22, 2.0E-22, 3.3E-22, 5. OE-22 , START 8 8.0E-22, l.lE-21 1.3E-21, 1.5E-21, 1 .RE-21 , 2.1E-21, 2.5E-21. START 9 2.9E-21, 3.4E-21 4.nE-21, 4.7E-21) START DATA ( (KABS03 ( LAMDA ) ♦LAMDA = 61,120) = START 1 4.8E-21* 4.0E-21 3.3E-21, 2.SE-21, 2.4E-21 , 2.0E-21, 1.7E-21, START 2 1.5E-21, 1.2E-21 l.lE-21. 7.8E-22, 5 .OE-22, 3.2E-22. 2.0E-22, START 3 46(0. OE- 00) ) START DATA ( (KABS02(LAMDA) ,LAMDA = 1,120) = START 1 6.0E-19, 6.0E-18 1.4E-17, 2.0E-17, 2.3E-17, 2.3L-17, 2.1E-17, START 2 1.5E-17» 9.0E-18 3.8E-18, l.lE-18, 3.4E-19, l.OE-20. 7.2E-18, START 3 1.3E-17, 8.8E-18 3.3E-18» 5.3E-19, 1.8E-20, 2.8E-23 . 1.6E-23, START 4 l.OE-23, 6.6E-24 4.7E-24. 1.5E-24, 5.0E-25, 94( 0.0 ) ) START DATA ( (TAB02( JZTAE 5) , JZTAB=1 ,87) = START 1 1.61E09, 1.70E09 1.79E09, 1.88E09, 1.98E09, 2.n8E09, 2.20E09, START 2 2.32E09, 2.46E09 2.6nE09, 2.74E09 , 2.89En9, 3.04E09, 3.22E09, START 3 3.41E09, 3.60E09 3.81E09, 4.02E09 , 4.26E09, 4.52E09, 4.78E09, START 4 5.15E09. 5.52E09 6.00E09, 6.63E09, 7.26E09, 8.15E09, 9.04E09, START 5 I.OIEIO. 1.16E10 1.31E10, 1.53E10, 1.75E10, 2.04E10, 2.50E10, START 6 2.96E10, 3.80E10 4.64E10, 6.00E10, e.oiEic, 1.03E11 , 1.38E11, START 7 1.74E11. 2.26E11 3.12E11. 3.98E11, 5.61E11, 7.24E11, 1.08E12, START 8 1.54E12, 2.00E12 3.06E12, 4.14E12, 6.04E12, 9.50E12, 1.30E13, START 9 1.90E13, 2.76E13 3.98E13. 5.70E13, 8.06E13, 1.02E14, 1.54E14, START 1 2.10E14, 2.78E14 3.68E14, 4.78E14, 6.16E14, 7.86E14, 9.98E14, START 2 1.27E15» 1.60E15 2.02E15, 2.58E15 , 3.30E15, 4.26E15 , 5.52E15, START 3 7.22E15, 9.52E15 1.27E16, 1.69E16, 2.28E16, 3.06E16, 3.80E16, START 4 5.64E16. 7.66E16 7.66E16) START DATA ( (TAB03{ JZTAE n , JZTAB=1 ,87 ) = START 1 l.OOEOO, 1.30E00 1.90EOO, 2.40E00, 3.10E0n, 4.00E00, 6.0OE00, START 2 8.00EOO, 1.05E01 1.40E01, 1.95E01, 2.60E01, 3.50E01 . ^.50E01. START 3 6.20E01. 9.00E01 1.30E02» 1.80E02, 2.40E02, 3.40E02. A.80E02 , START 4 7.00E02. 1.10E03 1.70E03 . 2 . 2 E 3 , 3.20E03, 5.00E03. 7.20E03, START 5 1.00E04, 1.50E04 2.00E04, 3.00E04, 4.00E04, 5.80E04, 8.00E04, START 6 1.01E05. 1.60E05 2.10E05. 3.00E05, 3.90E05, 5.00E05, 7.00E05, START -87- 7 9.20E05, 1.20E06. 1.70E06, 2.00E06* 2.60E06. 3.50E06* ^.50E06. START 8 5.70E06» 7.50E06, 9.00E06. 1.10E07» 1.30E07, 1.7nE07, 2.10E07, START 9 2.60E07, 3.50E07, ^.8nE07, 7.00E07» 1.00E08* 1.4nE08, 1.90E08. START 1 2.90E08, 4.00E08, 5.60E08, 8.50E08» 1.30E09* 2.00E09, 2.90E09* START 2 4.10E09, 6.00E09, I.OOEIO, 1.50E10. 2.20E10, 3.60E10» 5.60E10, START 3 8.50E10, 1.40E11, 2.50E11, 6.00E11. 1.40E12. 3.00E12* 4.00E12. START 4 5.00E12* 6.00E12* 6.00E12) START IFLAG4 = START IFLAG5 = START DO 10 M1 = 3.NSPEC '■ START READ (60,999) M2 . ( I ATOM { I TCH , M2 ) , I TCH = 1 , 1 ) START 10 CONTINUE START READ (60,994) T , TPR I NT ,TM AX 1 , TMAX2 , TQU I T , I FLAG8 , START 1 IFLAG1,EQUIL, JOSHUA START READ (60,993) LAT , LONG , TMONTH , TDAY , THOUR . TM I N , TSEC START DO 20 M3=2,NSPEC START READ {IUNIT,998) M4 START IF (M4 .EQ. 3) IFLAG6=1 START IF (M4 .EG. 4) IFLAG7=1 START IF (M4 .EG. 99) GO TO 30 ' START YNEW{M4) = 2.0 START READ (IUNIT,997) ( YHOLD ( M4 , JZ2 ) , JZ 2 = 1 » JZEND ) START 20 CONTINUE START PRINT 996 START IFLAG9 = 1 START GO TO 9000 START 30 IF (lUNIT .EG. 13) REWIND 13 START IF (IFLAG6 .EG. 1 .AND. IFLAG7 .EG. 1) GO TO 40 START IF (Z(l) .LE. 200.0 .AND. Z(JZEND) .GE. 30.0) GO TO 40 START PRINT 995 START IFLAG9 = 1 START GO TO 9000 START 40 DO 60 JZ3=1»JZEND START DO 50 JZTAB=1 ,87 START ZTAB = 202-2*JZTAB START IF (ZTAB .GT. Z(JZ3)) GO TO 50 START Y3SUN(JZ3) = (TAB02( JZTAB-1 )-TAB02 ( JZTAB) )*(Z( JZ3)-ZTAB)/2.0 START 1 + TAB02(JZTAB) START Y4SUN(JZ3) = (TAB03 (JZTAB-1 )-TAB03 (JZTAB) )»(Z(JZ3 )-ZTAB)/2.0 S'.ART 1 + TAB03(JZTAB) START GO TO 55 START 50 CONTINUE START PRINT 995 START IFLAG9 = 1 START GO TO 9000 START 55 IF (IFLAG6 .EG. 0) YHOLD(3,JZ3) =■ Y3SUN(JZ3) START IF (IFLAG7 .EG. 0) YH0LD(4,JZ3) = Y4SUN(JZ3) START 60 CONTINUE START 80 DO 85 M6=5.NSPEC START IF (YNEW(M6) .GT. 1.0) GO TO 85 START DO 83 JZ6=1,JZEND START YHOLD (M6,JZ6) = 5 . 0E-6*YHOLD ( 3 , JZ6 ) START 83 CONTINUE START 85 CONTINUE START DO 110 ITCH=1,10 START CHECK = 0.0 START DO 100 JZ5=1.JZEND START SUMIN( ITCH»JZ5) = 0.0 START DO 90 M5=3»NSPEC START SUMIN( ITCH»JZ5) = SUM I N ( I TCH , JZ 5 ) + YHOLD ( M 5 . J Z5 ) * I ATOM ( I TCH ,M5 ) START 90 CONTINUE START CHECK = CHECK + SUM I N ( I TCH , JZ 5 ) START 100 CONTINUE START IF (CHECK .GT. 1.0) GO TO 110 START ITCHXX = ITCH-1 START 60 TO 120 START 110 CONTINUE START 120 TOLD = T START DO 150 M = 1 ♦NSPEC START DO 140 JZ = l.JZEND START YPRINT{M,JZ) = YHOLD(M,JZ) START YLAST(M,JZ) = YHOLD(K,JZ) START 140 CONTINUE START 150 CONTINUE START 999 FORMAT( I2.2X,10( U .IX) ) START 998 FORMAT ( 12 ) START 997 FORMAT (7(E8.2,2X)) START 996 FORMAT (IH ,*N0 END OF START PROFILES FOUND.*) START 995 FORMAT (IH ♦»THE ALTITUDES ARE OUTSIDE THE BOUNDS OF 30-200 KM SO START lYOULL HAVE TO PUT IN YOUR OWN PROFILES FOR 02 AND 03.*) START 994 FORMAT ( 5 ( E8 . 2 . 2X ) , 2 ( 1 1 . 2 X ) . E8 . 2 » 2 X . 1 1 ) START 993 FORMAT ( 2 ( E9 . 2 . 2X ) , 5 ( I 2 . 2 X ) ) START 9000 RETURN START END START SUBROUTINE COMMON /A/ 1 BRK(4) FLUX(120) FTERMW(40 H2HOLD(20 IGAINA(20 IL0SS(4 ) ILOSSD(20 K0(51) KABS03( 12 L3(51 ) 20) 20) 1 N(51) NUSRUS(40 QRK( 20) RRK 21 ) TAB02(87) TERMW( 100 YNEW(20) Y3SUN( 20) COMMON /B/ 0) 2 3 4 5 6 7 INTEGER HC 1 TMONTH.TD REAL K0,N, IF (IFLAG3 RAD = PI/1 DEC = 23.4 IFLAG3 = 1 PRINT 998 JCNTRL = THOLD = TP TPRINT = TZEN = 15. ALPHA = (1 1 + CO CALL HANDL JCNTRL = 1 TPRINT = T GO TO 100 10 IF (JOSHUA TZEN = 15. ALPHA = (1 1 + CO 100 IF (ALPHA IF (ALPHA DO 200 JZl Y3SUN( JZl ) Y4SUN( JZl) ZENANG A( 51 ) .CRK ( 4) .FLXINF( 12 »G( 20) .H3HOLD(20 » IGAINB( 20 . ILOSSA(20 , INFORM( 30 .KABS( 20.1 ,KRK(80) .LL1(21) .NASRAS( 10 .PHOLD( 20. .Rl { 51 ) .RR2(21) .TAB03{ 87) ) .Y( 20) ,YOLD(20) .Y4SUN( 20) ALPHA , EQUIL IFLAG3 IFLAG9 JFLAG6 NFREAC, NREAC THOUR , TMAXl TPRINT, TOUIT OUNT.PHOLD.Rl ,R AY, THOUR. TMIN.T KABS.K.KRK.LAT. .NE. 0) GO TO 80.0 45*SINF( ( 2.0*PI RINT 0) ) ,20) ,20) ) 20) 00) 20) H IFLAG I JZX JOSHU 4, .ARK (4 ) .F (20) ,FRATE(20) ,GG(20) ,1 ATOM( 10,20) ,IGAINC(20.20) , ILOSSB(20,20) , JPULL (20) ,KABSN2( 120) ,L1 (51 ) .LL2(21) .NFTERM( 20) ,PROHLD( 71 ,20 ) ,R2(51 ) ,RR3(21 ) ,TEMP( 20) .YHOLD(20,20) .YPRINT(20.20) .Z ( 20) , HCOUNT, IFLAG6, lUNIT , JZMAIN, P I TWIN , ZBOTUM, IFLAG5 , , ITCHXX, A, JZEND , , NSPEC , , TMAX2 , , TSEC , 2,P3,RR1,RR2,RR3,RX1 SECTX LONG,K5,K5MAX,KABS0 2.KA 10 IFLAGl IFLAG7 JCNTRL LAT T TMONTH ZSTEP ,B(51 ) ,FF(20 ) ,FTERM(40) ,HBAR(20) . IGAIN(4) , IGAIND(20 ,ILOSSC(20 ,K( 1000) ,KABS02 ( 12 ,L2( 51 ) ,LL3(21 ) ,NTERM( 20 ) .QOLD(20) ,R3( 51 ) ,SUMIN( 10. .TERM( 100) ,YLAST(20, .YSTART ( 20 ,ZPRINT (20 IFLAG2 IFLAG8 JFLAG2 LONG TDAY TOLD ZTOP ,20) ,20) 20) ,20) ) RX2,RX3' iS03,KAESN2 30.4*TMONTH+TDAY-80.0 ) )/365.0) *ABSF(FLOATF(TSEC+60»(TMIN+6 0*THOUR))-4 3200.0)/3600.0 ,0/RAD)*ACOSF(SINF( L AT*R AD ) *S I NF ( DEC ^R AD) SF(LAT*RAD)*COSF(DEC*RAD)*COSF{TZEN*RAD)) E HOLD .EG. 1) GO TO 100 *ABSF(FLOATF(TSEC+6 0* ( TM I N+6 0*THOUR ) ) -4 3200.0)736 00.0 .0/RAD)*ACOSF(SINF (LAT*RAD)*SINF(DEC*RAD) SF( LAT*RAD)»COSF(DEC*RAD ) *COSF ( TZEN*RAD ) ) .GT. 90.0+SGRTF{Z( 1) ) ) GO TO 9000 .GT. 90.0) GO TO 300 = 1 , JZEND = YH0LD(3,JZ1) = YH0LD(4,JZ1) ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG ZENANG -90- 200 CONTINUE ZENANG 300 IF (Y3SUN(2) .GT. Y3SUN(1)) GO TO 310 ZENANG PRINT 997 ZENANG IFLAG9 = 1 ZENANG GO TO 9000 ZENANG 310 IF (Y3SUN(1) .GT. 0.0) GO TO 320 ZENANG H2H0LD(1) = 1.0 ZENANG GO TO 330 ZENANG 320 H2H0LD(1) = ( Z ( 2 ) -Z ( I) ) /LOGF { Y3 SUN ( 1 ) / Y3SUN ( 2 ) ) ZENANG 330 IF (YASUN(2) .GT. Y^SUN(l)) GO TO 3'^0 ZENANG PRINT 996 ZENANG IFLAG9 = 1 ZENANG GO TO 9000 ZENANG 340 IF (Y4SUN(1) .GT. 0.0) GO TO 350 ZENANG H3H0LD(1) = 1.0 ZENANG GO TO 400 ZENANG 350 H3H0LD(1) = ( Z ( 2 ) -Z ( 1 ) ) /LOGF ( Y4SUN ( 1 ) / Y4SUN ( 2 ) ) ZENANG 400 ADD3 = Y3SUN( 1 )*H2H0LD( 1) ZENANG ADD4 = Y4SUN( 1 )*H3H0LD( 1 ) ZENANG DO 410 JZ2=2.JZEND ZENANG ADD3 = ADD3 + ( Y3SUN ( JZ2- 1) + Y 3 SUN ( JZ2 ) ) ♦ ( Z ( JZ2 - 1 ) -Z ( JZ2 ) ) / 2 . ZENANG ADD4 = ADD4 + ( Y4SUN ( JZ2- 1 ) +Y4 SUN ( JZ2 ) ) * ( Z ( JZ2 -1 ) -Z ( JZ2 ) ) / 2 . ZENANG H2HOLD(JZ2) = ADD3/Y3SUN( JZ2 ) ZENANG H3HOLD(JZ2) = ADD4/Y4SUN( JZ2) ZENANG 410 CONTINUE ZENANG 997 FORMAT (1H0.*THE TOP SCALE HEIGHT FOR 02 IS NEGATIVE*) ZENANG 996 FORMAT (1H0.*THE TOP SCALE HEIGHT FOR 03 IS NEGATIVE*) ZENANG 998 F0RMAT(1H1»* STARTING VALUES *,/» ZENANG 1 *(REACTION TERM BLOCK CONTAINS RATE COEFFICIENTS*) ZENANG 9000 RETURN ZENANG END ZENANG- -91- SUBROUTINE PRODUC COMMON /A/ A(51) 1 BRK(4) FLUX ( 120) FTERMW( 40 ) H2H0LD( 20 ) IGAINAI 20»20 ) IL0SS(4) ILOSSD(20.20) K0(51 ) KABS03( 120) L3(51 ) 1 N(51 ) 2 NUSRUS(400) 3 QRK{ 20) 4 RR1(21) 5 TAB02(87) 6 TERMW(IOO) 7 YNEW(20) 8 Y3SUN(20) COMMON /B/ 2 3 ♦CRK(4) ,FLXINF( 120) ♦G( 20) .H3H0LD( 20 ) * IGAINB(20»20 ♦ILOSSA(20»20 . INFORM( 30 ) .KABS( 20»120 ) »KRK ( 80 ) .LL1(21 ) »NASRAS( 1000) .PHOLD(20,20) »R1 ( 51 ) .RR2(21 ) »TAB03(87) ,Y( 20) ♦YOLD( 20) .Y4SUN(20) ALPHA , EQUIL « f I FLAG3 ^ IFLAG9, JFLAG61 NFREAC, THOUR . TPRINT IFLAG4 IJZX ■ JOSHUA NREAC ■ TMAXl TQUIT IFLAG5 ITCHXX JZEND NSPEC TMAX2 TSEC INTEGER HC0UNT,PH0LD,R1 .R2.R3»RR1, 1 T MON TH,T DAY, THOUR, TM IN ♦T SEC. TX REAL KO,N,KABS»K,KRK , L AT , LONG , K 5 , K CALL Q9EXUN IF (NFREAC .EG. 0) GO TO 9000 ALPHAR = ALPHA*PI/180.0 REARTH = 6371,0 IF (ALPHA .GT. 90.0) GO TO 100 THIS CALCULATION DONE USING EON. 53 VOL. 12, NO. 8, 1964. PI2 = (PI/2.0)**2 XX = (REARTH+Z( JZMAIN) ) /H2H0LD( JZM AXLE = ( 1.0-PI2»(0. 115 + 1. 0/LOGF(PI F02 = EXPF ( ( ALPHAR»*2/2.0 )/ ( 1.0-AL XX = (REARTH+Z( JZMAIN) )/H3H0LD( JZM AXLE = ( 1.0-PI2*(0. 115+1. 0/LOGF(PI F03 = EXPF ( (ALPHAR GO TO 300 100 IF (ALPHA .GT. 90. + SQRTF ( Z ( JZM A I N V = (6370,0+Z( JZMAIN) )*SI NF ( ALPHAR IF (V .GE. Z( 1 ) ) GO TO 130 IF (V .LE. Z(JZEND)) GO TO 140 DO 110 JZ1 = 1 , JZEND IF (V .GT. Z( JZl ) ) GO TO 120 110 CONTINUE GO TO 140 120 H02 = (H2H0LD( JZ1)-H2H0LD( JZl-1 ) )* RK(4) (20) RATE (20 ) G(20) ATOM( 10, GAINC(20 LOSSB(20 PULL{20 ) ABSN2( 12 1(51) L2(21) FTERM(20 R0HLD(71 2(51 ) R3(21) EMP( 20) HOLD(20, PRINT(20 (20) » HCOUNT , IFLAG6 , lUNIT , JZMAIN , PI , TMIN , ZBOTUM RR2 ,RR3 , 10) i20) ^20) ,B(51) ,FF(20) ,FTERM(40) ,HBAR (20 ) , IGAIN{ 4) , IGAIND(2n ,ILOSSC(20 ,K ( 1000 ) ,KABS02 ( 12 ,L2(51 ) ,LL3(21 ) ,NTERM( 20) ,QOLD( 20 ) ,R3( 51 ) ,SUMIN ( 10, ,TERM( 100) ,YLAST( 20, ,YSTART (20 ,ZPRINT (20 IFLAGl, IFLAG2. IFLAG7, IFLAG8' JCNTRL, JFLAG2. '0) i20) 20) ,20.) ) LAT T 5 TMONTH, ZSTEP , LONG TDAY TOLD ZTOP RXl ,RX2»RX3' ^ ,KAB502,KABS03,KABSN2 N SWIDER, PLANET. SPACE SCI., AIN) *XX/2.0) ) ) /PI2*»2 PHAR**2»(0.115+AXLE*ALPHAR**2 AiN) *XX/2.0 ) ) ) /PI2**2 2/2.0)/(1.0-ALPHAR»*2*(O.115+AXLE*ALPHAR*»2 ) ) ) GO TO 1000 )-6370.0 (V-Z(JZ1-1 ) ) /(Z( JZl )-Z( JZl-1 ) PRODUC f PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC , PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC ) ) ) PRODUC PRODUC PRODUC ) ) ) PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC ) PRODUC -92- MV-Z ( JZl-1 ) ) / (Z( JZl )-Z( JZl-1 [ERF+1 ) ) 1 +H2H€)LD(JZ1-1 ) H03 = (H3H0LD( JZ1)-H2H0LD( JZl-1 ) 1 +H3H0LD( JZl-1 ) GO TO 150 130 H02 = H2H0LD( JZl ) H03 = H3H0LD( JZl ) GO TO 150 140 H02 = H2H0LD( JZEND) H03 = H3H0LD( JZEND) 150 X02 = (ALPHA-90.0)/SORTF(HO2) X03 = (ALPHA-90.0) /SQRTF(H03) !F (X02 .GT. 2.5) GO TO 180 ERF02 = X02 FXX = 1.0 DO 160 IERF=1,50 FXX = FXX»IERF TXX = (-l)**IERF»X02**(2*IERF+l )/(FXX* (2^ ERF02 = ERF02 + TXX IF (ABS(TXX) .LT. ABS ( ERF02* 1 . OE-4 ) ) GO TO 170 160 CONTINUE 170 ERF02 = ERFO2*2.0/SQRTF(PI) GO TO 190 180 ERF02 = 1.00 190 IF (X03 .GT. 2.5) GO TO 220 ERF03 = X03 FXX = 1.0 DO 200 IERF2=1,50 FXX = FXX*IERF2 TXX = (-1)**IERF2*X03**(2*IERF2+1)/(FXX»(2*IERF2+1)) ERF03 = ERF03 + TXX IF (AB5(TXX) .LT. ABS ( ERF03* 1 . OE-A ) ) GO TO 210 200 CONTINUE 210 ERF03 = ERFO3*2.0/SORTF(PI ) GO TO 230 220 ERF03 = 1.0 230 F02 = {101.4/SQRTF(H02) )*(1.0 + ERF02) F03 = (101.4/SQRTF(H03) )*(1.0+ERF03) 300 DO 310 LAMDA1=1,120 FLUX(LAMDAl) = FLX I NF ( LAMDA 1 ) *EXPF ( 1 -KABS02(LAMDA1 ) ♦Y3SUN ( JZMA I N ) *H2H0LD ( JZMA I N ) * 1 . 0E5*FO2 1 -KABS03(LAMDA1 ) *Y4SUiN ( JZMA I N ) *H3H0LD ( JZMA I N ) * 1 . 0E5*FO3 ) 310 CONT I NUE DO 330 II=1,NFREAC FRATE( II) = 0.0 DO 320 LAMDA2=1.120 l-RATE(II ) = FRATE( II ) 320 CONTINUE 330 CONTINUE GO TO 9000 1000 DO 1010 II = l.NFREAC FRAtE( II) = 0.0 1010 CONTINUE 90C0 CALL R9EXUN RETURN END KABS( I I ,LAMDA2 ) »FLUX ( L AMDA2 ) PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC PRODUC- -93- SUBROUTINE PHYSICS COMMON /A/ A(51) 1 BRK(4) 2 FLUX(120) 3 FTERMW(40) 4 H2HOLD(20) 5 IGAINA(20.20) 6 IL0SS(4) 7 ILOSSD(20.20) 8 KOI 51 ) 9 KABS03( 120) X L3(51 ) 1 N(51 ) NUSRUS(400) ORK( 20) RRl(21) TAB02(87) TERMWI 100) YNEW(20) Y3SUN( 20) COMMON /B/ .CRK(4) .FLXINFC 120) »G(20) .H3HOLD(20) , IGAINB(20.20 » ILOSSA(20.20 ,INFORM( 30 ) .KABS(20.120 ) .KRK(80) »LL1(21) .NASRAS( 1000) »PHOLD(20»20) »R1(51) »RR2(21 ) ♦TAB03( 87 ) . Y ( 2 ) »Y0LD(20) .Y4SUN( 20) ALPHA , EQUIL , H IFLAG5 ITCHXX ,ARK(4) »F (20) »FRATE(2 ,GG(20 ) , I ATOM! 1 , IGAINC( ♦ILOSSB( ,JPULL (2 ,KABSN2( ♦ LI (51 ) ,LL2(21 ) »NFTERM( ,PROHLD( ♦ R2(51 ) ,RR3(21) ,TEMP( 20 ,YH0LD(2 ,YPRINT( ,Z (20) . HCOU IFLA lUNI JZMA P I TMIN ZBOT 0) 20) 20) ) 20) 2 IFLAG3, IFLAG4, 3 IFLAG9, IJZX , 4 JFLAG6» JOSHUA, JZEND 5 NFREAC, NREAC , NSPEC 6 THOUR , TMAXl , TMAX2 7 TPRINT, TQUIT » TSEC INTEGER HC0UNT,PH0LD,R1 »R2.R3,RR1,RR2.RR 1 TMONTH,TDAY»THOUR,TMIN.TSEC.TX REAL ICO,N,KABS,K ,KRK , LAT . LONG » K5 » K5MAX ,K DO 800 M1=3»NSPEC I I IX = (Ml-3 )*NREAC IIFX = (Ml-3 )*NFREAC NNSTOP = NTERM(Ml) + NFTERM(Ml) NTX = NTERM(Ml) DO 700 JZ1 = 1, JZEND DO 100 IL00K1=1,4 IL0SS( ILOOKl ) = IGAIN( ILOOKl ) = ICO CONTINUE DO 600 IL00K2=1»4 ILMAX = IGMAX = TLMAX = TGMAX = 0.0 DO 500 NN=1»NNST0P IF (NN .GT. NTX) GO TO 110 NX = NASRAS( I I IX+NN ) GO TO 150 110 NX = NUSRUS( I IFX+NN-NTX ) 150 ILX = IL00K2-1 IF (NN .GT. NTX .AND. NX .GT. NFREAC) GO IF (NN .LE. NTX .AND. NX .GT. NREAC) GO DO 200 IL00K3=1 , ILX IF (NN .GT. NTX) GO TO 180 IF (NX .EG. ILOSS( IL00K3) ) GO TO 500 GO TO 200 ) 0,2( 20,; NT , G6, T . IFLAGl IFLAG7 JCNTRL LAT T TMONTH ZSTEP .B(51) .FF(20) ,FTERM(40) ,HBAR(20) , IGAIN (4) ,IGAIND(20 ,ILOSSC(20 »K(1000) ,KABS02( 12 ,L2(51 ) ,LL3(21 ) ,NTERM( 20) ,QOLD(20) ,R3(51 ) ,SUMIN( 10, ,TERM( 100) ,YLAST ( 20, ,YSTART (20 .ZPRINT (20 IFLAG2 IFLAG8 JFLAG2 LONG TDAY TOLD ZTOP .20) ■ 20) 20) ,20) UM 3»RX1,RX2»RX3 ABS02»KABS03.KABSN2 TO 300 TO 300 PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS -9k- 180 200 250 310 320 350 400 420 500 610 620 630 700 800 IF (NX+NREAC .EQ CONTINUE IF (NN .GT. NTX) IF ( PROHLD(NX ♦JZ GO TO 500 IF (PROHLD(NX+NR TLMAX = PROHLD(N ILMAX = NX+NREAC GO TO 500 TLMAX = PROHLD(N ILMAX = NX GO TO 500 IF (NN .GT, NTX) NX = NX-NREAC GO TO 320 NX = NX-N'^REAC 00 400 IL00K4=1 . IF (NN .GT. NTX) IF (NX .EQ. IGAI GO TO 400 IF (NX+NREAC .EU CONTINUE IF (NN .GT. NTX) IF {PROHLD(NX»JZ GO TO 500 IF (PROHLD(NX+NR TGMAX = PROHLD(N IGMAX = NX+NREAC GO TO 500 TGMAX = PROHLD(N IGMAX = NX CONTINUE IF (IL00K2 .EQ. IF (TLMAX .GE. IF (TGMAX .GE. GO TO 600 ILOSS(l) = ILMAX IGAIN(l) = IGMAX BIGL = TLMAX BIGG = TGMAX CONTINUE IF (ILOSS(l) .EQ IL0SSA(M1 , JZl ) = IL0SSB(M1, JZl) = IL0SSC(M1,JZ1 ) = IL0SSD(M1, JZl) = GO TO 620 IL0SSA(M1 , JZl )=I IF ( IGAIN( 1 ) .EQ IGAINA(M1 . JZ 1 ) = IGAINB(M1 , JZl ) = IGAINC(M1, JZl ) = IGAIND(Mi , JZ 1 ) = GO TO 700 IGAINA(M1 . JZl )= I CONT INUE CONTINUE RETURN END IL0S5( I L00K3 ) ) GO TO 500 GO TO 260 1) .GT. TLMAX) GO TO 270 EAC»JZ1) .LE. TLMAX) GO TO 500 X+NREAC.JZl ) GO TO 310 ILX GO TO 350 N( IL00K4) ) GO TO 500 . IGAIN( I L00K4 ) ) GO TO 500 GO TO 410 1 ) .GT. TGMAX ) GO TO 420 EAC.JZl) .LE. TGMAX) GO TO 500 X+NREAC.JZl ) 1) GO TO 550 .1*BIGL) ILOSS( IL00K2 ) = ILMAX .1*BIGG) IGAIN( IL00K2 ) = IGMAX . 0) GO TO 610 IL0S5( 1 )+100 ILOSS( 2 ) +100 IL0SS( 3)+100 IL0SS( 4)+100 L0SSB{M1, JZl )=IL0SSC{M1 ,JZ1 . 0) GO TO 630 IGAIN{ 1 )+lCO IGAIN (2 )+100 IGAIN( 3)+100 IGAIN(4 )+100 L0SSD(M1 ,JZ1 GAINB(M1,JZ1)=IGAINC(M1,JZ1)=IGAIND(M1,JZ1) PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSICS PHYSIC- -95- SUBROUTINE HANDLE COMMON /A/ 1 BRK(4) 2 FLUX(120) ?! FTERMW(40) 4 H2HOLD(20) 5 I6AINA( 20»20 ) 6 IL0SS(4) 7 ILOSSD(20,20) 8 K0(51 ) 9 KABS03( 120) X L3( 51 ) 1 N( 51 ) 2 NUSRU5(400) 3 QRK( 20) A RRl(21) 5 TAB02(87) 6 TERMW(IOO) 7 YNEW(20) 8 Y3SUN(20) COMMON /B/ A( 51 ) CRK(4) FLXINFC 120) G( 20) H3HOLD(20 ) IGAINB(20»20) ILOSSA(20,2C) INFORM( 30 ) »KABS( 20»120 ) ,KRK{80) .LL1{21 ) »NASRAS( 1000 ) »PHOLD( 20»20 ) .Rl ( 51 ) .RR2(21 ) ♦TAB03( 87) ,Y(20) .YOLD(20 ) »Y4SUN( 20) ALPHA , EQUIL , H IFLAG4 I JZX JOSHUA NREAC TMAXl TQUIT 2 IFLAG3» IFLAG4, IFL 3 IFLAG9, IJZX , ITC 4 JFLAG6» JOSHUA, JZE 5 NFREAC, NREAC , NSP 6 THOUR » TMAXl , TMA 7 TPRINT, TQUIT , TSE INTEGER HC0UNT,PH0L0»R1 »R2.R3,R ] TMON TH.TDAY, THOUR, TM IN ,TSEC,TX REAL K0,N,KABS,K,KRK,LAT,LONG,K TX=TPRINT+TSEC+TMIN*60+THOUR»36 TMONTH = TX/25920G0 TX = TX-TMONTH*2592000 TDAY = TX/B6400 TX = TX - TDAY*86400 THOUR = TX/3600 TX = TX - THOUR*3600 TMIN = TX/60 TX = TX - TMIN*60 TSEC = TX ITSEC = TSEC + 100 ITMIN = TMIN + 100 ITHOUR = THOUR + 100 ITDAY = TDAY + 100 PRINT 999, JCNTRL, ITHOUR, ITMIN, ITIME = KLOCK( 1 ) TIME=ITIME/1000. PRINT 966, TIME PRINT 998, T PRINT 997, ALPHA PRINT 994 NTREAC = NREAC + NFREAC IF (NTREAC .LT, 8) GO TO 10 PRINT 993 ARK(4) F (20) FRATE(20) GG(20) I ATOM( 10,2 IGAINC{20, ILOSSB(20, JPULL (20) KABSN2( 120 Ll(51 ) LL2(21) NFTLRM( 20) PR0HLD(71, R2(51) RR3(21) TEMP( 20) YHOLD(20,2 YPRINT(20, Z{20) , HCOUNT, IFLAG6, lUNIT , JZMAIN, PI TMIN , ZBOTUM, ,B(51 ) ,FF( 20) ,FTERM(40) ,HBAR (20) ) ,IGAIN(4) 0) ,IGAIND(20 0) ,ILOSSC(20 ,!<( 1000) ,KABS02 ( 12 ,L2( 51 ) ,LL3(21 ) ,NTERM( 20) D) ,Q0LD(20) ,R3( 51 ) ,SUMIN( 10, ,TERM(100) ) ,YLAST{20, 0) ,YSTART(20 ,ZPRINT (20 IFLAGl, IFLAG2 20) 20) 20) ,201 ) IFLAG7 JCNTRL. LAT T TMONTH. Z5TEP . IFLAGS^ JFLAG2. LONG , TDAY , TOLD . ZTOP R1,RR2,RR3,RX1,RX2,RX3 5, K 5 MAX ,KABS02,KABS03,KABSN2 00+TDAY*86400+TMONTH*2 5 9 2000 ITSEC, 1 TDAY, TMONTH HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE -96- IF (NTREAC PRINT 992 IF (NTREAC PRINT 991 IF (NTREAC PRINT 965 IF (NTREAC PRINT 964 IF (NTREAC PRINT 963 IF (NTREAC PRINT 962 IF (NTREAC PRINT 961 10 DO 20 JZ = PRINT 982» IF (NTREAC PRINT 989» IF (NTREAC PRINT 988 , IF (NTREAC PRINT 987» IF (NTREAC PRINT 989. IF (NTREAC PRINT 988, IF (NTREAC PRINT 987, IF (NTREAC PRINT 989, IF (NTREAC PRINT 988, 20 CONTINUE PRINT 986 IF (NSPEC PRINT 985 IF (NSPEC PRINT 984 IF (NSPEC PRINT 983 30 DO 40 JZ2= PRINT 990, IF (NSPEC PRINT 989, IF (NSPEC PRINT 988, IF (NSPEC PRINT 987, 40 CONTINUE IF (JCNTRL PRINT 975 IF (NSPEC PRINT 974 IF (NSPEC .LT. 16) GO TO 10 .LT. 24) GO TO 10 .LT. 32 ) GO TO 10 .LT. 40) GO TO 10 .LT. 48) GO TO 10 .LT. 56) GO TO 10 .LT. 64) GO TO 10 1 , JZEND Z(JZ), HBAR(JZ), { PROHLD( II ,JZ ) , I 1 = 1 »7 ) .LT. 8) GO TO 20 (PROHLD( 12, JZ) , 12 = 8, 15 ) .LT. 16) GO TO 20 ( PROHLD( I 3,JZ ) , 13=16,23 ) .LT. 24) GO TO 20 (PROHLD( 14, JZ) ,14=24,31 ) .LT. 32) GO TO 20 (PROHLD( 12, JZ) ,12 = 32, 39) .LT. 40 ) GO TO 20 (PROHLD( I3,JZ) ,13 = 40,47) .LT. 48 ) GO TO 20 (PROHLD( 14, JZ), 14=48, 55) .LT. 56) GO TO 20 (PROHLD( 12, JZ) , 12 = 56,63 ) .LT. 64) GO TO 20 (PROHLD( 13, JZ) ,13 = 64,72) .LT. 10) GO TO 30 .LT. 18) GO TO 30 .LT. 26) GO TO 30 1 , JZEND ZPRINT(JZ2), (YPRINT(I5,JZ2),I5=2,9) .LT. 10) GO TO 40 ( YPRINT ( 16, JZ2 ), 16=10,17 ) .LT. 18) GO TO 40 (YPRINT ( 17, JZ2 ), 17=18,25 ) .LT. 26) GO TO 40 (YPRINT ( 18, JZ2 ) , 18 = 26,33) .EG. 0) GO TO 9000 .LT. 10) GO TO 50 .LT. 18) GO TO 50 HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE -97- 26) GO TO 50 10) GO TO 70 18) GO TO 70 PRINT 973 IF (NSPEC .LT PRINT 972 50 DO 60 JZ3=1.JZEND PRINT 971 »ZPR1NT ( JZ3) « ( ILOS 1 IL0SSC(M1 ,JZ3 ) »IL0SS0(M1, IF (NSPEC .LT. 10) GO TO 60 PRINT 970, ( ILOSSA (M2 ,JZ3 ) , 1 ILO.SSD(M2,JZ3) ♦M2=10.17 ) IF (NSPEC .LT. 18) GO TO 60 PRINT 969, ( I LOSSA (M3,JZ3 ) , 1 IL0SSD(M3,JZ3) ,M3=18,25 ) IF (NSPEC .LT. 26) GO TO 60 PRINT 968, ( IL0SSA(M4,JZ3 ) , 1 lL0SSD(M/i, JZ3) ,M4 = 26,33 ) 60 CONTINUE PRINT 967 IF (NSPEC .LT PRINT 974 IF (NSPEC .LT PRINT 973 IF (NSPEC .LT. 26) 60 TO 70 PRINT 972 70 DO 80 JZ4=1,JZEND PRINT 971 ,ZPRINT ( JZ4) , ( IGAI 1 IGAINC(M5, JZ4) ,IGAIN0(M5, IF (NSPEC .LT. 10) GO TO 80 PRINT 970, ( IGAINA( M6, JZ4 ) , I 1 IGAIND(M6, JZ4) ,M6=10,17 ) IF (NSPEC .LT. 18) GO TO 80 PRINT 969, ( IGAINA(M7,JZ4 ) , I 1 IGAIND(M7,JZ4) ,M7=18,25 ) IF (NSPEC .LT. 26) GO TO 80 PRINT 968, ( IGAINA( MB, JZ4) , I 1 IGAIND(M8,JZ4) ,M8=26,33 ) 80 CONTINUE IF (JZEND .GT WRITE (13,399 GO TO 110 90 IF (JZEND .GT WRITE (13,898 GO TO 110 100 WRITE ( 13,397 110 CONTINUE RFWIND 13 897 FORMAT( I2,/,2(7(F8.2,2X) ,/) 898 FORMAT ( I? , /,7 ( E8.2 ,2X ) ,/ , 7( 899 FORMAT ( 12 ,/,7 (E8.2,2X ) ) 961 FORMAT ( IH , * 1M(66) TERM(67) TERM 2(71)*) 962 F0RMAT(1H , ♦ 1(58) TERM(59) T£RM( 263)») SA(M1,JZ3),IL0SSB(M1,JZ3) JZ3 > ,M1=2 ,9) LOSSB(M2,JZ3) ,IL0SSC(M2,JZ3) LOSSB(M3,JZ3) ,ILOSSC (M3,JZ31 ILOSSBt M4,JZ3 ) , I LOSSC ( M4 , JZ 3 1 7) GO TO 90 (M, ( YHOLD(M, NA(M5,JZ4),IGAINB(M5,JZ4) JZ4) ,M5=2 ,9 ) GAINB(M6,JZ4),IGAINC(M6,JZ4) GAINB(M7,JZ4),IGAINC{M7,JZ4) GAINB(M8,JZ4) ,IGAINC(M8,JZ4) JZ) , JZ=1 ,7 ) ,M=3,N5PEC ) ,99 14) GO TO 10 (M, ( YHOLD(M, JZ),JZ=1,14),M=3,NSPEC),99 (M, ( YHOLD( M,JZ) ,JZ=1,20) ,M = 3, NSPEC) ,99 ,&( E8.2 ,2X ) ) E8.2,2X ) ) TERM(64) TtRM(65) (68) TERM(69) TERM(70) TERM(56) TERM(57) 60) TERM(61) T£RM(62) HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE HANDLE TERHANDLE TERMHANDLE HANDLE TERMHANDLE TERM(HANDLE HANDLE 963 FORMAT ( IH ♦ ♦ TERM( 48 ) TERM(49) TEHANDLE 1RM(50) TERM(51) TERM(52) TERM( 53 ) TERM( 54) TERHANDLE 2M(55 )*) HANDLE 964 FORMAT ( IH » * TERM (40) TERM ( 41 ) TERHANDLE 1M(42) TERM(43) TERM( 44) TERM(45) TERM ( 46 ) TERMHANDLE 2 ( 4 7 ) ♦ ) HANDLE 965 FORMAT (IH . ♦ TERM ( 32 ) TERM( 33) TERMHANDLE 1(34) TERM(35) TERM( 36) TERM( 37 ) TERM( 38 ) TERMIHANDLE 239)*) HANDLE 966 F0RMAT(1X,*C0MPUTER TIME= »,F9.3,* SECONDS*) HANDLE 967 FORMATdHl .*MAJOR PRODUCTION REACTIONS * »/ ,* HEIGHT NO. HANDLE 12 NO. 3 NO. 4 NO. 5 NO. 6 HANDLE 2 NO. 7 NO. 8 NO. 9 - > HANDLE 968 FORMATdH ,8X,8( 3X, 12, */*,2( 12,*,* ,12)) HANDLE 969 F0RMAT(1H .7X»8(3X,I2> «/-»• ,2 ( 12,*,* ,12)) HANDLE 970 FORMATflH »6X»8(3X.I2, */*,2( 12,*,* ,12) ) HANDLE 971 FORMAT (1H0,F6.1»2X,8(I2,*/».2(I2»* ,*) , 12, 3X) ) HANDLE 972 FORMAT (IH ♦* NO. 26 NO. 27 NO. 28 HANDLE 1 NO. 29 NO. 3r NO. 31 NO . 32 NO . 33*)HANDLE 973 FORMAT ( IH ,* NO . 18 NO. 19 NO. 20 HANDLE 1 NO. 21 NO. 22 NO. 23 NO. 24 NO. 25*) HANDLE 974 FORMAT( IH ,* NO. 10 NO. 11 NO. 12 HANDLE ]N0. 13 NO. 14 NO. 15 NO. 16 NO. 17*) HANDLE 975 FORMAT( 1H1,*MAJ0R DF5TRUCTI0N REACTIONS -*,/,* HEIGHT NO. HANDLE 12 NO. 3 NO. 4 NO. 5 NO. 6 HANDLE 2 NO. 7 NO. 8 NO. 9 i^ ) HANDLE 982 FORMAT ( lH0,f^6.1 ,2X,E1 0.4 ,3X ,7 ( ElO .3,3X) ) HANDLE 9 83 FORMAT ( IH ,* Y ( 26 ) Y(27) Y ( 2 8 ) YHANDLE 1(29) Y(30) Y{ 31 ) Y( 32) Y ( 33 )*) HANDLE 984 FORMATdH ,* Y ( 1 8 ) Y( 19 ) Y(20 ) Y( HANDLE 121) Y(22) Y(23 ) Y (24 ) Y( 25 )*) HANDLE 985 FORMATdH , * Y (1 ) Ydl ) Y( 12 ) Y ( l HANDLE 13 ) Y( 14) Y ( 15 ) Y ( 1 6 ) Yd7)*) HANDLE 986 FORMATdHl,//,* SPECIE S DENSITIES- *,/,♦ HEIGHT Y{2) YHANDLE 1(3) Y (4 ) Y(5 ) Y (6 ) Y(7 ) Y(HANDLE 18) Y(9)*) HANDLE 9 87 FORMAT ( IfH ,12X,8(E10.: ,3X ) ) HANDLE 988 FORMATdH ,11X,S(£10.; ,3X ) ) HANDLE 989 FORMAT( IH ,10X,8(E10.: ,3X ) ) HANDLE 990 FORMAT ( IHO , F6 . 1 , 2 X , F 1 1 . 8 , 3 X , 7 ( E 10 . 3 , 3 X ) ) HANDLE 991 FORMATdH , * TERM( 24 ) TERM(25) TEHANDLE 1RM(26) TLRM{27) TERM{ 28 ) TERM( 29 ) TERM( 30) TERHANDLE 2M( 31 )*) HANDLE 992 F0RMAT(1H , » TERM ( 16 > TERM ( 17 ) TERHANDLE 1M(18) TERM(19) TERM( 20 ) TERM ( 71 ) TERM( 22 ) TERMHANDLE 2(23)*) HANDLE 993 F0RMAT(1H ,» TERM( 3 ) TERM( 9) TERM (10) TERMHANDLE 1(11) TERM (12) TERM( 13) TERM( 14) TERM( 15 )*) HANDLE 994 FORMAT ( IHO ,*REACTI0N TERMS * ,/,* HEIGHT STEP SIZE TERMHANDLE 1(1) TERM( 2 ) TERM( 3 ) TLRM(4 ) TERM( 5 ) TERM(HANDLE 26 ) TERM ( 7)* ) HANDLE 997 F0RMAT{1H ,»ZENITH ANGLE = ♦,F6.2, * DEGREES*) HANDLE 998 F0RMAT(1H ,*RUNNING TIME = *,E10.3 ,* SECONDS*) HANDLE 999 FORMAT ( 11 ,*LOCAL TIME (HOURS, MINUTES, SECONDS) = *, 12,*/*, 12, */* , I2HANDLE 1 ,* ON DAY NUMBER *, I 2 , *, MONTH NLiMBER ♦, I 3 ) HANDLE 9000 RETURN HANDLE END HANDLE -99- SUBROUTINE RUNKUT COMMON /A/ A(51) 1 BRK(4) 2 FLUX(1?0) 3 FTERMW(40) 4 H2HOLD(20) 5 IGAINA(20>20 6 I LOSS (4) 7 ILOSSD(20»20 8 KOISI ) 9 KABS03(120) X L3(51 ) 1 N(51) 2 NUSRUS(40 0) 3 QRK( 2^) 4 RRl (21) 5 TAB02(87) 6 TERMW(lOO) 7 YNEW(20) 8 Y3SUN{20) COMMON /B/ 2 ALPHA 1FLAG3 IFLAG9 JFLAG6 NFREAC THOUR TPRINT CRK(4) FLXINFC 12 G( 20) H3H0LD( 20 1 G A I N B ( 2 ILOSSA(20 INFORM( 30 KABS(20»1 KRK(80) LL1(21) NASRAS( in PHOLD( 20» Rl (51) RR2 (21 ) TAB03( 87) Y(20) YOLD{ 20) Y4SUN( 20) EQUIL IFLAG IJZX JO SHU NREAC TMAXl TQUIT 20) 20) H IFLAG5 I TCHXX JZEND NSPEC TMAX2 TSEC ,*RK(4) ,F(20) .FRATE (20 ,GG(20 ) ,1 ATOMdO ,IGAINC(2 »IL0SSB(2 ,JPULL (20 .KABSN2( 1 ,L1 (51 ) ,LL2(21 ) ,NFTERM( 2 ,PROHLD( 7 ,R2(51 ) ,RR3(21 ) ,TEMP( 20) ,YHOLD (20 ,YPRINT(2 ,Z (20) , HCOUN IFLAG I UN IT JZMAI PI TMIN ZBOTU 10) .20) .20) •n) 20) . IFLAGl . IFLAG7 ■ JCNTRL LAT . T TMONTH . Z5TEP ♦R(51) ,FF(20) ,FTERM( 40 ) ,HBAR (20) , IGAIN ( 4) , IGAIND(20 ,IL0SSC(20 »!<( 1000) ,KABS02 ( 12 ,L2(51 ) ,LL3(21 ) ,NTERM( 20) ,Q0LD(2n) ,R3(51 ) ,SUMIN( 10 » .TERM( 100 ) ,YLAST( 20, ,YSTART(20 ,ZPRINT (20 IFLAG2 IFLAG8 JFLAG2 LONG TDAY TOLD ZTOP '20) '20) 20) ,2o: ) 2,R3,RR1,RR2»RR3,RX1,RX2,RX3, ■SECTX . LONG, K 5, K 5 MAX ,KABS0 2,KABS03,KABSN2 ,0E-6) GO TO 20 ,0E-6) H = TPRINT-T+TOLD INTEGER HC0UNT,PH0LD,R1 ,R: 1 TMONTH, TDAY, THOUR, TMIN, T; REAL K0,N,KABS,K,KRK,LAT CALL Q9EXUN H = HBAR ( JZMAIN ) *3.0 10 TCHECK = T-TOLD+H IF (TCHECK .LT. TPRINT-1 IF (TCHECK .GT. TPRINT+1.( IFLAG2 = 1 20 DO 30 Ml=2, NSPEC YOLD(Ml) = Y{M1) QOLD(Mi) = QRK(Ml) 30 CONTINUE 40 DO 90 JRK=1 ,4 DO 50 11 = 1 , NREAC TERM( Il + NREAC) = K ( I 1* JZE ND- I JZX ) * Y ( L 1 ( I 1 ) ) *Y ( L2 ( 1 1 ) ) * Y ( L3 ( I 1 ) ) TERM(Il) = T£RM( I 1+NREAC ) 50 CONTINUE DO 52 II 1 = 1, NFREAC FTERM( III) = FRATt ( I FTERM( i I 1 + NFREAC ) = FTERMI 52 CONTINUE JRKX = ( JRK-1 )*NSPEC DO 70 M2=3, NSPEC MJRK = JRKX+M2 IIIX = (M2-3)*NREAC DY = 0.0 ( LLl (III)) 1 ) RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT NREAC) GO TO 55 -NN ) ) ^NN) ) FTERM(NUSRUS( I IFX + NF ) ) DY NNSTOP = NTERM(M2) DO 60 NN=1»NNST0P IF (NASRAS( I I IX+NN) .LE DY = DY+TERM(NASRAS( III; GO TO 60 55 DY = DY - TERM(NASRAS( I 60 CONTINUE IIFX = (M2-3 )*-NFREAC NFSTOP = NFTERM(M2) DO 65 NF=1»NFST0P IF (NUSRUS( I IFX + NF) .LE. NFREAC) GO TO 63 DY = DY + FTERM(NUSRUS{ I IFX + NF) ) GO TO 65 63 DY = DY - 65 CONTINUE KRK(MJRK ) 70 CONTINUE DO 80 M3=3»NSPEC MJRK = JRKX+M3 XXRK = ARK( JRK)» (KRK (MJRK )-BRK ( JRK )*QRK(M3 ) ) Y(M3) = Y(M3) + H»XXRK QRK(M3) = QRK(M3 )+3.0*XXRK-CRK ( JRK )*KRK(MJRK) 80 CONTINUE 90 CONTINUE K5MAX = 0.01 DO 100 M4=3,NSPEC IF (Y{M4) .LT. 0.0) Y{M4) = 0.0 IF (ABS(KRK(MA+NSPEC) ) .LT. l.OE-10 .OR. 1 ABS( 1.0-KRK(M4)/KRK(NSPEC+M4) ) .LT. l.OE-2 .OR. 2 ABS(KRK( M4 )+KRK( NSPEC+M4 ) +KRK ( 2*NSP EC + M4 ) +1<',R K (3*NSPEC+M4) )*H 3 .LT. Y0LD(M4)*1.0E-4) GO TO 100 K5 = ABS( (KRK(NSPEC+M4)-KRK(2*NSPEC+M4) ) / ( KRK ( M4 ) -KRK ( NSPEC+M4 ) IF (K5 .GT. K5MAX) K5MAX = K5 100 CONTINUE IF (K5MAX .LT. 0.10) GO TO 120 H = H*0.08/K5MAX IFLAG2 = DO 110 M5=2.NSPEC Y(M5) = Y0LD(M5) QRK(M5) = Q0LD{M5) 110 CONTINUE GO TO 40 120 T = T+H H = H*0.08/K5MAX HCOUNT = HCOUNT+1 IF ( IFLAG2 .EQ. 0) GO TO 10 CALL R9EXUN RETURN END RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT RUNKUT- SUBROUTINE WEIRD COMMON /A/ 1 BRK(4) 2 FLUX(120) 3 FTERMW(40) 4 H2HOLD(20) A( 51 ) .CRK(4) .FLXINFI 12( .G( 20) •H3H0LD( 20 : 5 IGAINA( 20.20) .IGAINB(20! 6 IL0SS(4) 7 IL0SSD( 20»20) 8 K0( 51 ) 9 KABSO3(120) X L3(51 ) N(51 ) NUSRUSC 400) QRK( 20) RRl{21 ) TAB02(87) TERMW( 100 ) YNEW(20) 8 Y3SUN(20) COMMON /B/ 2 3 4 5 6 7 . ILOSSA(20 ♦ INFORM( 30 ♦KABS( 20,1 ,KRK ( 80) .LL1(21) ,NASRAS( 10 .PHOLD( 20. ,R1( 51 ) .RR2(21) .TAB03(87) .Y(20) .YOLD(20) .Y45UN( 20) EQUIL [FLAG^ IJZX ALPHA IFLAG3 IFLAG9 JFLAG6. JOSHU/ NREAC TMAXl TQUIT H NFREAC THOUR TPRINT INTEGER HCOUNT.PHOLD.Rl .R2. 1 TMONTH.TDAY. THOUR .TMIN.TSE REAL KO.N.KABS.K .KRK.LAT.LO AKRACY = l.OE-3 HMIN = l.OE-6 YNEW(l) = 1.0 CALL Q9EXUN H = HBAR ( JZMAIN)*3.0 IF (H .LT. TPRINT- (T-TOLD H = TPRINT-( T-TOLD) IFLAG2 = 1 10 DO 100 11=1. NREAC TERM(Il) = K( I 1*JZEND-IJZ TERM( Il+NREAC) = TERM(Il) 100 CONTINUE DO 200 111 = 1. NFREAC FTERM(IIl) = FRATE ( I II )*Y ( LLl ( I I 1 ) ) FTERM( I Il + NFREAC ) = FTERM(IIl) 200 CONTINUE DO 500 M1=3.NSPEC F(Ml ) = 0.0 G(M1 ) = 0.0 NNSTOP = NTERMCMl) II IX = (M1-3)*NREAC DO 300 NN=1. NNSTOP NASX = NASRAS( I I IX+NN) IF (NASX .GT. NREAC) GO TO 250 .ARK (4 ) ,F (20) .FRATE(20) ,GG(20) ,1 ATOM( 10. . IGAINC(20 .1 LOSSB(20 . JPULL (20 ) .KABSN2( 12 .LI (51 ) ,LL2(21 ) ,NFTERM( 20 .PROHLD( 71 .R2(51 ) .RR3(21 ) ,TEMP( 20) ,YHOLD(20. ,YPRINT(20 .Z (20) . HCOUNT IFLAG6 lUNIT J Z MA IN PI TMIN ZBOTUM .B(51 ) .FF(20) ,FTERM(40) .HBAR ( 20) ) ,IGAIN(4) D) .IGAIND(20 D) ,ILO5SC(20 .K (1000) ,KABS02 ( 12 .L2 ( 51 ) .LL3( 21 ) ,NTERM( 20) 3) .OOLD(20) ,R3(51 ) .SUMIN ( 10, .TERM( 100) ) .YLAST(20. D) .YSTART(20 .ZPRINT (20 IFLAGl. IFLAG2 ■ 20) ■ 20) 20) 20) .20) ) IFLAG5 I TCHXX JZEND N5PEC TMAX2 TSEC R3,RR1,RR2.RR3 C.TX NG.K5,K5MAX ,KA IFLAG7, JCNTRL, LAT T 1 TMONTH, ZSTEP , IFLAG8. JFLAG2. LONG . TDAY • TOLD ■ ZTOP .RXl ,RX2,RX3^ iS02.KAB503.KABSN2 ) ) GO TO 10 X)*Y(L1(I1))*Y(L2(I1))*Y(L3(I1)) WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD G(M1) = G(M1) + TERM(NASX) WEIRD GO TO 300 WEIRD 250 F(M1) = F(M1) + TERM(NASX) WEIRD 300 CONTINUE . WEIRD IIFX = {M1-3)»NFREAC WEIRD NFSTOP = NFTERM(Ml) WEIRD DO 400 NF=1,NFST0P WEIRD NUSX = NUSRUSt I IFX + NF) WEIRD IF (NUSX .GT. NFREAC) GO TO 350 WEIRD G(M1) = G(M1) + FTFRM(NUSX) WEIRD GO TO 400 WEIRD 350 F(M1) = F(M1) + FTERM(NUSX) WEIRD 400 CONTINUE WEIRD IF (Y(M1) .LT. l.OE-100) GO TO 450 WEIRD G(M1) = G(M1)/Y(M1) WEIRD GO TO 500 WEIRD 450 G(M1) = 0.0 WEIRD 500 CONTINUE WEIRD 600 DO 700 M2=3.NSPEC WEIRD IF (G(M2) .LT. l.OE-lOn .AND. F(M2) .LT. l.OE-100) GO TO 675 WEIRD IF {G(M2)*H .LT. l.OE-8) GO TO 675 WEIRD IF (G(M2>»Y(M2 )*1.0E8 .LT. F(M2)) GO TO 650 WEIRD IF (F(M2)*1.0E8 .LT. G ( M2 ) ♦ Y ( M2 ) ) GO TO 625 WEIRD YNEW(M2) = ( Y (M2 )-F ( M2 ) /G ( M2 ) )*EXPF (-G (M2 >*H) + F(M2)/G(M2) WEIRD GO TO 700 WEIRD 625 YNEW(M2) = Y ( M2 ) *EXPF ( -G ( M2 ) *H ) WEIRD GO TO 700 WEIRD 650 YNEW(M2) = Y ( M2 ) + F(M2)*H WEIRD GO TO 700 WEIRD 675 YNEW(M2) = Y ( M2 ) WEIRD 700 CONTINUE WEIRD DO 800 I2=1.NREAC WEIRD TERMW(I2)=K(I2*JZEND-IJZX)*YNEW(L1(I2))*YNEW(L2(I2))*YNEW(L3(I2)) WEIRD TERMW( I2 + NREAC) = TERMW(I2) WEIRD 800 CONTINUE WEIRD DO 900 112 = 1. NFREAC WEIRD FTERMW(II2) = FRATE ( I I2)*YNEW(LL1{ I 12) ) WEIRD FTERMW( I I2 + NFREAC ) = FTERMW(II2) WEIRD 900 CONTINUE WEIRD DO 1200 M3=3.NSPEC WEIRD FF(M3) = 0.0 WEIRD GG(M3) = 0.0 WEIRD NNSTOP = NTERM(M3) WEIRD IIIX = (M3-3)*NREAC WEIRD DO 1000 NN2=1. NNSTOP WEIRD NASX = NASRASC I I IX+NN2 ) WEIRD IF (NASX .GT. NREAC) GO TO 950 WEIRD GG(M3) = GG(M3) + TERMW(NASX) WEIRD GO TO 1000 WEIRD 950 FF(M3) = FF(M3) + TERMW(NASX) WEIRD 1000 CONTINUE WEIRD IIFX = (M3-3 )*NFREAC WEIRD NFSTOP = NFTERM(M3) WEIRD DO 1100 NF2=1. NFSTOP WEIRD -103- NUSX = NUSRUS( IIFX + NF2) WEIRD IF (NUSX .GT. NFREAC) GO TO 1050 WEIRD GG(M3) = GG(M3) + FTERMW(NUSX) WEIRD GO TO 1100 WEIRD 1050 FF(M3) = FF{M3) + FTERMW(NUSX) WEIRD 1100 CONTINUE WEIRD IF (YNEW(M3) .LT. l.OE-100) GO TO 1150 WEIRD GG(M3) = GG(M3 ) /YNEW(M3 ) WEIRD GO TO 1200 WEIRD 1150 GG(M3) = 0.0 WEIRD 1200 CONTINUE -, WEIRD DELMAX = 0.0 WEIRD DO 1400 M4=3,NSPEC WEIRD IF (Y(M4) .LT. l.OE-100) GO TO 1400 WEIRD IF (ABS( Y( M4)*G(M4)-F(M4) )*EXPF(-G (M4)*H)*1.0E2 .LT. F ( M4 ) ) WEIRD 1 GO TO 1400 WEIRD IF (GG(M4)*H .GT. 300.0) GO TO 1400 WEIRD IF (GG(M4) .LT. l.OE-100 .AND. FF(M4) .LT. l.OE-100) GO TO 1275 WEIRD IF (GG(M4)*H .LT. l.OE-8) GO TO 1275 WEIRD IF (GG(M4) *YNEW(M4 )*1.0E8 .LT. FF(M4)) GO TO 1250 WEIRD IF ( FF(M4) *1.0E8 .LT. GG ( M4 ) * YNEW ( M4 ) ) GO TO 1225 WEIRD YBACK = (YNEW(M4)-FF(M4) /GG(M4 ) )*EXPF(GG(M4)*H)+FF (^4)700(^^4) WEIRD GO TO 1300 WEIRD 1225 YBACK = YN EW ( M4 ) *EXPF ( GG { M4 ) *H ) , WEIRD GO TO 1300 WEIRD 1250 YBACK = YNEW(M4) - FF(M4)*H ", WEIRD GO TO 1300 WEIRD 1275 YBACK = YNEW(M4) WEIRD 1300 DEL = ABS( YBACK/Y(M4) - 1.0) WEIRD IF (DEL .GT. DELMAX) DELMAX = DEL WEIRD 1400 CONTINUE WEIRD IF (DELMAX .LE. AKRACY) GO TO 1500 WEIRD IF (H .LT. 1.01*HMIN) GO TO 1500 WEIRD H = H/2,0 WEIRD IF (H .LT. HMIN) H = HMIN WEIRD IFLAG2 = WEIRD GO TO 600 WEIRD 1500 DO 1700 M6=3»NSPEC ' WEIRD FBAR = ( F( M6 )+FF (M6 ) ) /2.0 WEIRD GBAR = (G(M6 )+GG(M6) ) /2.0 WEIRD IF (GBAR .LT. l.OE-100 .AND. FBAR .LT. l.OE-100) GO TO 1575 WEIRD IF (GBAR*H .LT. l.OE-8) GO TO 1575 WEIRD IF ( &BAR»Y (M6 )*1 .OEB .LT. FBAR) GO TO 1550 WEIRD IF (FBAR*1.0E8 .LT. GBAR*Y(M6)) GO TO 1525 WEIRD Y(M6) = ( Y(M6 )-FBAR/GBAR) *EXPF (-GBAR*H ) + FBAR/GBAR WEIRD GO TO 1700 WEIRD 1525 Y(M6) = Y( M6)*EXPF(-GBAR»H) WEIRD GO TO 1700 WEIRD 1550 Y(M6) = Y(M6) + FBAR*H WEIRD GO TO 1700 WEIRD 1575 Y(M6) = Y(M6) WEIRD 1700 CONTINUE WEIRD T = T+H WEIRD HCOUNT = HCOUNT + 1 WEIRD IF (IFLAG2 .EQ. 1) GO TO 9000 IF (DELMAX .GT. l,0E-8) GO TO 1800 H = H»1.0E2 GO TO 1900 1800 H = H*1.5 IF (H .LT. HMIN) H 1900 IF (H .LT. TPRINT H = TPRINT IFLAG2 = 1 GO TO 10 9000 CALL R9EXUN RETURN END SCOPE LOAD HMIN (T-TOLD) (T-TOLD) WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD WEIRD -105- SUBROUTINE MICROPLT COMMON /A/ 1 BRK(4) FLUX( 120) FTERMW(40 ) H2HOLD(20 ) IGAINA(20.20) IL0SS<4) ILOSSD(20,20) 8 K0(51 ) 9 KABS03( 120) X L3(51 ) 1 N(51) 2 NUSRUS(ifOO) 3 QRK(20) h RRl (21 ) 5 TAB02(87) 6 TERMW( 100 ) 7 YNEW (20 ) 8 Y3SUN(20) COMMON 2 3 i¥ 5 6 7 A(51 ) ■ CRK (4) .FLXINF( 12 .G( 20) iH3HOLD(20 I IGAINB(20 . ILOSSA(20 , INFORM( 30 .KAB5(20, 1 .KRK(80) iLLl (21 ) .NASRAS( 10 .PHOLD{20. ■ Rl (51 ) >RR2(21 ) .TAB03( 87 ) .Y(20) .YOLD(20) .Y4SUN( 20) EQUIL IFLAG I JZX ♦ARK(4) ,F (20) ,FRATE (20) ,GG(20) , IATOM( 10,2 ♦ IGAINC(20, ,1 L0SSB(2C, ,JPULL (20) ,KABSN2( 120 »L1{51 ) ,LL2(21 ) »NFTERM(20) ,PR0HLD(71, ,R2(51 ) ,RR3(21 ) ,TEMP{ 20) »YHOLD(20»2 »YPRINT( 20, ,Z (20) H , HCOUNT, IFLAG5, IFLAG6, ITCHXX, lUNIT , JZEND , JZMAIN, NSPEC , PI , TMAX2 , TMIN , TSEC , ZBOTUM, 3,RR1,RR2,RR3,R ,TX C,IX,I Y 18 , MICROPLT (2) , IXHOLD(20,20) ) ,(NCYCLE=0) ALPHA , IFLAG3, IFLAG9, JFLAG6, JOSHU NFREAC. NREAC THOUR , TMAXl TPRINT, TOUIT INTEGER HCOUNT, PHOLD, Rl ,R 1 TMONTH, TDAY, THOUR , TMIN, T COMMON/DD/IN, IOR,IT,IS,IC COMMON /DD ID/ ID (4) , IFL.LU DATA(ID=32HG W ADAMS X-39 DIMENSION TEXTK 10) ,TEXTV DATA(TEXTV=16H HEIGHT(KM) IFL = IFL+1 MINY = 2 0*( INT(Z(JZEND) )/20 MAXY=2 0*(INT(Z(l)-l)/20+l NTICS=(MAXY-MINY)/10+1 IF(NCYCLE) GO TO 20 NCYCLE=1 $ MAXPOWR=0 C LOOP 100 COMPUTES RANGE OF YHOLD. DO 100 ISPEC=3, NSPEC J=-l 10 J = J+1 $ IF( YHOLD( ISPEC, JZEND)/ lO.-i 100 IF ( J.GT.MAXPOWR ) MAXPOWR=J C LOOP 800 PLOTS UP TO 5 CURVES/FRAME 20 NOPLTS =(NSPEC+2)/5 DO 800 IPL0T=1 , NOPLTS IN=IS=1 CALL DDBOX( 70, 1020,70,1020) C LOOP 200 DRAWS VERTICAL TICS AND HEIGHTS. ,8(51 ) ,FF(20) ,FTERM(40) ,HBAR(20) ) ,IGAIN(4) 0) ,IGAIND(20 0) ,ILOSSC(20 ,i< ( 1000) ,KA3S02 ( 12 ,L2(51) ,LL3(21 ) ,NTERM( 20 ) 0) ,QOLD(20) ,R3( 51 ) , SUM IN (10, ,TERM( 100) ) ,YLAST(20, 0) ,YSTART(20 ,ZPRINT(20 IFLAGl, IFLAG2 20) 20) 20) ,20) ) IFLAG7, JCNTRL, LAT T TMONTH, ZSTEP , IFLAG3^ JFLAG2. LONG , TDAY , TOLD . ZTOP Xl,RX2,RX3i U.GT.l) GO TO 10 MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT -106- 25 200 DO 200 J=1,NTICS IY=70+9 50*(J-1)/(NTICS-1) IF(MOD( J»2) .EQ.O) GO TO 25 LABEL=MINY+10*( J-1 ) $ ENCODE ( 8 »920 » TEXT ) LABEL FORMAT( I4,4X ) IX=16 $ CALL DDTAB $ CALL DDT ABNA8 ( 1 » T EXT . 1 ) IF( J.EQ. l.OR.J.EQ.NTICS) GO TO 200 IX=70 $ ICC=1R- $ CALL DDTAB $ CALL DDTABCC IX=1020 $ ICC=1R- S CALL DDTAB $ CALL DDTABCC CONTINUE VERTICAL LABELS AND TIC MARKS CALL DDXY $ IY=1020 $ CALL DDTAB CALL DDTAE C DRAW VERTICAL LINES AND LABEL GRAPH. IN=0 $ CALL DDSEGVEC LIMIT=MAXP0WR-1 DO 300 M = l , LIMIT IX=70+950*M/MAXPOWR $ IY=70 $ 3 00 CALL DDXY IN=I0R=1 $ IS=2 $ IX=0 $ IY=460 CALL DDTABNA8 (2 »TEXTV.l ) S IOR=0 LABEL=10 S ENCODE(8»920.TEXT ) LABEL DO 400 M = l , LIMIT IX=36+950*M/MAXPOWR $ IY=45 $ IS=1 CALL DDTABNA8( l.TEXT tl ) ENC0DE(8.920,TEXT1 ) M IX=IX+20 + M/10»8 $ IY=55 $ IS=0 S CALL DDTAB 400 CALL DDTABNA8( l.TEXTl ,1 ) ENCODE(48.940»TEXT1 ) ALPHA »THOUR .TMIN ,TSEC »TMONTH,TDAY 940 FORMAT (*ALPHA = »F5. 1 »6X, 31 3»* M0NTH*I3* DAY*I3) IY=3 $ IX=150 S IS=2 $ CALL DDTAB CALL DDTA3NA8 (6 »TEXT1 . 1 ) IX=470 $ ICC=12B S CALL DDTAB $ CALL DDTABCC IX=518 $ ICC=12B $ CALL DDTAB S CALL DDTABCC C LOOP 500 SCALES LOG(YHOLD) TO PLOTTER UNITS. DO 500 J=l,400 V=YHOLD(J) $ IF(V.LT.l.) GO TO 500 IXHOLD( J )=ALOG10( V)*950./MAXPOWR + 70. 500 CONTINUE C LOOP 600(1) PLOTS EACH CURVE. IS=IN=0 LIMIT =MIN0(5.NSPEC+3-5»IP LOT) DO 600 1 = 1 .LIMIT IND=5*IPLOT+I-3 ENCODE(8.960»TEXT ) IND 960 FORMAT( I2.6X) IXTEMP = IXHOLD( IND, 1 ) $ I YTEMP= ( Z (1 ) -M I N Y ) *95 0/ ( MAX Y-M I N Y ) + 70 < C LOOP 600(J) PLOTS SPECIES SYMBOL AT EACH POINT AND DRAWS LINES C BETWEEN SYMBOLS. DO 600 J=1.JZEND IX=IXTEMP-8+IND/10*4 $ IY=IYTEMP CALL DDTAB $ CALL DDT ABNA8 ( 1 . TEXT , 1 ) IX=IXTEMP MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT 600 800 IF( J.EQ.JZFND) GO TO 600 CALL DDCONVEC IXTEMP = IXHOLD( IND» J + 1 ) IYTEMP=(Z ( J + 1 )-MINY )*950./(MAXY-MINY) + 70. XDIF=IXTEMP-IX $ YDIF=IYTEMP-IY D=SQRTF(XDIF*XDI F+YDIF*YDIF) IDX=8*XDIF/D $ IDY=8*YDIF/D IX=IX+IDX $ IY=IY+IDY $ CALL DDXY IX= IXTEMP-IDX $ IY= lYTEMP-IDY i CALL DDXY CONTINUE CALL DDFRAME RETURN END MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROPLT MICROP- -108- PENN STATE UNIVERSITY LIBRARIES AQDQD7EDED13S