DESCRIPTION OF THE COMPUTER PROGRAM Appendix F of Development of a Simulation Model for Estimating Ground Level Concentrations of Photochemical Pollutants Prepared by Systems Applications, Inc. Beverly Hills, California 90212 for the Air Pollution Control Office of the Environmental Protection Agency Durham, North Carolina 27701 ------- ~ESCRIPTION OF THE COMPU'l'ER PROGRAM Appendix F of. Development of a Simulation Model for Estimating Ground Level Concentrations of Photochemical Pollutants Steven D. Reynolds Report 71S2\I26 August 1971 Prepared by Systems Applications, Inc. Beverly Hills, California 90212 for the Air Pollution Control Office of the Environmental.protection Agency Durham, North Carolina 27701 under Contract CPA 70-148 ------- ACKNOWLEDGlotENTS ., . We wish to acknowledge the contributions of Mr. Clarence L. Nelson and Mr. Philip J. W. Roberts in the preparation of the Aircraft Emissions Program and the subroutines UPWIND and PLUME of: the Atmospheric Pollution Simulation Program respectively. ------- CON'i'ENTS '-' INTRODU-::'l'IOH . . . . . . . . . . . . . . . . r. THE ATMOSPHERIC POLLUTION SIMULATION P~OGf~ A". ~ne Structure of ths ~tmosphe~i~ Pollution Simulation Frog'ram ~ . . . 13. Progrfu~ Operation. . . . . . II. THE A:.;(CHAFr E:l!SSIONS PROGRAM . . . . . . . REFERENCES . . . . . . . . . . . . . . . p.age . . . F- 1 . F- 2 . . F- 3 . . . . . F- 8 . . . F-l8 . F-22 ------- INTRODUCTION The purpose of this Appendix is to describe in. a general way the various facets of the computer program that comprise the urban airshed model. The development of this program is an involved undertaking, as a wide variety of calculations are required both in the numerical solution of the governing equations and in the handling and manipulation of source .and meteorological data. . Emphasis was placed on the preparation of an efficient program and on insuring the applicability of the program in any urban airshed. Computational efficiency is vital in hro respects. First, it is important to minimize the number of calculations that must be carried' out in any large-scale computing effort, particularly when these calcula- tions are a part of iterative loops or are of a repetitive nature. Second, the numerical solution of three dimensional, coupled, time dependent partial differential equations is, apart from data input, handling, and processing considerations, a very time consuming calcula- tion, and efficiency is paramount. With regard to the general applica- bility of the program, an airshed is essentially defined by its terrain, meteorological, and source inputs. Since these parameters are all treated as input data to the program, no difficulty should be encounter- ed in applying the program in a variety of urban airsheds. At this time, the most cumbersome aspect of program usage is the large effort required in the preparation of the meteorological input data (see Appendix C). Inccrporation of an interpolation scheme into the program which automatically constructs wind and inversion maps from the raw monitoring station data would alleviate this difficulty. Such an addition to the existing program is required prior to the validation of the model for a large number of different meteorological coni'.itions. It is also worthwhile to remind the reader that, in general, a sizable effort is needed for the preparation of the source emissions inventory for any large urban area. Ho~.,ever" it should be noted that, while the spatial and temporal distributions of pollutant emissions are generally the same from day to day (with the exception of the weekend), the meteoro- logy is likely to undergo daily variations. Thus, the emissions inven- tory must be undertaken only once. If the model is to be used in testing control strategies, it will become necessary to alter the source input data. Such alterations should present few difficulties, however, from the standpoint of program operation. Two programs are currently being used in the Los Angeles Basin simulation effort. The main program, the Atmospheric Pollution Simula- tion Program, calculates pollutant concentrations as 'a function of space and time using source and meteorological inputs. We discuss both the structural and ele operational aspects of thi~ program in Section I of this Appendix. Aircraft emissions are treated as input to this program via magnetic tape. The computer program which generates this tape, The Aircraft Emissions Program, is discussed in Section II. ------- 1:. THE ATMOSPHERIC POLLUTION SIl.fiJLATION PROGRAM The Atmospheric Pollution Simulation Program accepts source and meteorological data as inputs to produce pollutant concentration/ time histories at discrete points distributed throughout the airshed. The code consists of a MAIN program and several subroutines which perform specialized tasks. In its present configuration, the program utilizes two magnetic tape units and temporary disc storage. Current- ly, the program size is approximately 350K bytes. Inputs to the program are contained on cards and two magnetic tap~s. The following information is supplied via punched cards: . the horizontal extent and topography of the airshed . the time span over which the simulation is to take place the spatial and temporal grid' ~ize " the time interval for concentration'printout the chemical species to be sinrolated the maximum number of iterations allowed in Step I . of the numerical solution of the governing equations ,freeway and non-frce""ay vehicle mileage distributions fixed source emissions the chemical reaction rate constants the initial concentration distribution . . . . . . . Tbe two magnetic tapes mentioned above contain meteorological and aircraft emissions data respectively. The meteorological data consists of wind speed and direction, mixing depth, and temperature for each colUmn of grid points. Each set of data is applicable for a certain time interval (say, between 0630 and 0730 PST), and this time interval is recorded on the tape immediately preceding each set of data. Before ~ simulation begins, the program reads the tape until it locates a set of data applicable at the starting time of the simulation. Similarly, the aircraft emissions tape contains sets of data preceded by a time interval. This tape is also automatically positioned before the simu- . tat ion begins. The aircraft emissions tape ,is generated by the Aircraft Emissions Program, discussed in Section II, and consists of pollutant emissions rates (ppm/minute) into each cell of the airshed. The primary output of the program is a printout of the hourly- ~veraged ground-level concentration distribution for each chemical species simulated. These concentration maps, in stan~ardized output format, ~an be easily, compared with observed pollutant levels. It should be noted, how- ever, that the number of concentration values that are computed is gen- erally quite large. Thus, we do not attempt to display the entire con- centration field at any time. It is felt that generally the ground-level concentrations will be of most interest since measurements of pollutant levels 'aloft are scarce. With this in mind, we have included an output option which allows the user to monitor the ~rogress of the simul~tion with an appropriate level of detail. The program prints the instantaneous F-2 ------- "--- g~ound~level concentrations and user-selected vertical concentration profiles at regular time intervals. The user chooses ground-level coordinates of points of interest, and the program prints all concen- trations between the ground and inversion base at those points. In simulations performed thus far, we have sampled the instantaneous concentrations every twenty minutes and have chosen to exanline the vertical concentration profiles only above monitoring stations. We have found this display of output to be quite satisfactory for this phase of model development. In future phases of progr~m development, we plan to investigate alternative means of presenting the results of a simulation. Sucha p~esentation will prob~ly include the use of graphical techniques to disp~ay ground-level concentration conto~ maps and vertical concen- tration profiles. Moreover, the present b;eatment of data input to the program is by no means final, and, in fact, only represents a con- venient means for satisfying current needs. When the model is com- p~eted, treatment of input will be re-evaluated. At that time we will be. able to make recommendations as to how the ,input should be handled. A foreseeable change in the input structure will occur, for example, if'we develop an 'interpolation scheme to construct the wind field directly from the monitoring station dat~. The remainder of this discussion will be segmented into two parts. In Section A we describe the overall structure of the program and dis- cuss the role played by each subroutine subprogram. Section B is devoted to a genc~alized, step-by-step description of how the program op~rates during a typical simulation run. A'. The Structure of the Atrnospheric Pollution Simulation Program In this section we describe the various subprograms the main program, giving details of how they operate and appropriate portions of the main text and appendices for cussions of the underlying concepts involved. In Figure cate schematically the overall structure of the program, Y-l:we list each subprogram and summarize its function. that comprise referring to detailed dis- F-l we indi- and in Table MAIN MAIN is responsible for two aspects of a simulation: (1) (2) the control of model initialization the n~merical integration of the governing air shed equations Model initialization consists of defining the physical characteristics of theairshed (horizontal ~xtent and topography), reading source emissions data, positioning data tapes, est~lishing the grid on which Figure F-l and Table F-l follo~' . F-3 ------- "IJ I ~ SOURCE I , GJ FI~URE F-l. STRD~URE 9F THE ATMPSPHERIC P9LWTIPN snruLATION PROGRAM BCCONC AVCOOC ( ) START GJ j i V ICCONC 1 . 1 \ GUCOF (METEOR \ / . - -r-1 SIMDIG PLUME ------- TABLE F-l. ROUTINE MAIN AVCONC BCCONC DIFCOF ICCONC METEOR PLUME PRINTl SIMDIG SOURCE UPWIND SUMMARY OF SUBPROGRAMS AND THEIR FUNCTION FUNCTION initialization of the model and the numerical solution of the governing airshed equations calculates and prints hourly-averaged ground level pollutant concentration maps supplies boundary concentrations at points of horizontal inflow into the model calculates the turbulent eddy diffusivity at each point in a given column of grid points sets up the initial concentration distribution reads lower wind field and inversion base elevatioIE and alters the vertical concentration distribution due to temporal changes in the inversion base elevations calculates the spatial distribution of power plant et:\issions prints instantaneous ground-level pollutant concentration maps as well as user-selected vertical concentration profiles solves system of linear equations ~mose coefficient matrix is tridiagonal calculates automotive and fixed source pollutant emissions (exc~pt power plants) calculates the upper wind field F-5 ------- . the nUInerical integration is to take place, and entering the initial concentration distribution. However, not all of these functions take place in MAIN. The sequencing involved in model initialization is described in more detail in section B. When model initialization is complete, the nUInerical integration procedure may be initiated. Before an integration time step can be taken, source and meteoro- logical variables. must be updated. Subroutines SOURCE, METEOR, 'and PLUME are called to supply vehicular and fixed source emissions, meteorological data, and power plant emissions respectively. Also, aircraft emissions are obtained from the temporary disc area, and subroutineDIFCOF is called to supply values of the eddy diffusivity. The numerical integration now proceeds with the sequential solution of Steps I, II, and III as described in Appendix.D. .. , When the simulation is complete, the program punches the final concentration distribution onto cards. ,This enables the user to re- start the integration at a later date. In early stages of model valida- tion, we have found it convenient to divide a 12-hoursimulation into three four-hour model runs, t~ereby allo~ring us to monitor the progress of the simulation. If, after a four-hour segment, the predicted concen- trations are found to be acceptable, we merely use the punched concen- trations as initial conditions to continue the simulation further in time. AVCONC Subroutine AVCONC is called from MAIN after every integration time step and is responsible for the calculation and printout of hourly- . averaged ground level concentration maps. Each time AVCONC is entered, an array is incremented by the current values of the ground level con- centrations. After one hour, this array is divided by the n~ber of time steps taken over that hour so as to produce an array of hourly- averaged pollutant levels. By manually superposing the actual monitor- ing station observations on these maps, one can easily evaluate the model's performance. BCCONC Subroutine BCCONC .supplies pollutant levels at points on the air shed boundary where horizontal inflow takes place. Generally, the actual con- centrations are not known, and hence. the values used are merely educated guesses. If levels were kno~m, then they could be input on cards or tape. ~n current simulation efforts involving CO, we have employed boundary con- centrations of 3 ppm at all points of horizontal inflow. DIFCOF In Step I of the numerical integration and during the cubic spline interpolation procedure, it is necessary to have values of the turbulent eddy diffusivity. Subroutine DIFCOF supplies values of the diffusivity F-6 ------- for.a vertical column of grid points. The coordinates of the column are transferred to DIFCOF, and diffusivities are calculated from knowledge of the meteorological conditions in the column (see Appendix C for a discussion of the algorithm employed). ICCONC The initial concentration distribution is established through sub- routine ICCONC. This subroutine assumes two forms, depending on the mode of initial concentration input. In the first mode, the initial ground-level concentrations are read from cards which were prepared from ground-level concentration maps that are based on actual observed pollutant levels. The values in the ground cells are extended verti- cally to the inversion to complete the definition of the initial con- centration field. The second mode of input a~ises when a simulation begins at a point in time at which a previous simulation terminated. As was mentioned in the discussion of the MAIN program, at the end of a simulation the entire concentration field is punched onto cards. Thus the second mode simply consists of a READ statement to read the punched cards. 14ETEOR The processing of the meteorological data is the domain of subroutine METEOR. During model initialization, METEOR is called by ~~IN, and the meteorological input tape is read until data applicable at the start of the simulation are encountered. As the simulation progresses, the meteor~ ological variables are updated as necessary by reading the tape. See Appendix C for a further discussion of the meteorological data. When inversion elevations change, it beco~es necessary to transfer the concen- tration field from the old grid to the new grid. A cubic spline proce- dure is employed which requires values of the diffusivity as supplied by DIFCOF. A system of linear equations must be solved, and subroutine SIMDIG performs this task. Finally, subroutines UPWIND and PLUME are called to supply the upper wind field and the spatial distribution of power plant emissions, respectively. PLUME Subroutine PLUHE computes the spatial and temporal distribution of p~1er plant emissions. Since the spatial distribution varies with changes in meteorological variables (every hour in this case), it. was convenient to call this subroutine from METEOR. A detailed discussion of the method used to distribute the power plant emissions may be found. in Section III of Appendix A. PRINTl The results of the numerical integration are printed in subroutine PRINTI. User-selected vertical concentration profiles and ground-level F-7 ------- concentration maps are printed at r,egular time intervals. 1'le have found that an interval of 20 minutes petween printouts gives a good indication of hourly trends and does not result in endless pages of unused computer output. Predicted vertical concentration profiles in the vicinity of monitoring stations are also printed out at the option of the user. SIMDIG Subroutine SI~IDIG solves a system of linear equations whose co- efficient matrix is tridiagonal. Thomas' algorithm is employed in the solution (see von Rosenberg (1969) for details). The three bands of the matrix and the y~own vector are transferred to the sQ~routine, and the solution is returned upon completion. , SOURCE Auto:notive and fixed source emissions (\'ith the exception of pO\oler plant emissions) are computed in subroutine SOURCE. All emissions are treated as surface fluxes. See Appendix A for a complete discussion of the source inventory. UPWIND Subroutine UPWIND computes the wind field aloft. ~~en the lower wind field is completely defined in HETEOR, control is passed to UPWIND. We refer the reader to Appendix C for a discussion of ,the algorithm employed. B. Program Operation In this section we discuss, in a general manner, the operation of the program during a typical si~ulation. Since neither the model nor the program have taken final form, a complete card-by-card description would be lengthy and of no particular value at this juncture. Figure F-2 should be referred to throughout this discussion, as a schematic repre- sentation of the order of the computational operations is shown therein. A simulation commences in ~~IN with the input of all information contained on punched cards, with the exception of the initial concentra- tion distribution, which is input later. The program then reads the aircraft emissions tape until it encounters data applicable at the time the simulation is to start. At this point the program calculates all the'parameters that will be used in the numerical integration of the governing airshed equations. For example, partial derivatives of the terrain elevation ah and ah ax ay F-8 ------- where h.(x,y) = terrain elevation at x,y are approximated by finite differences. After all necessary parameters are calculated, subroutine METEOR is called to position the meteorolog- ical data tape. As before, this is accomplished by reading the tape until data applicable at th~ start of the simulation is encountered. Control then passes to subroutine ICCONC where the initial conce~tra- tion field is established. Finally, subroutine SOURCE is ~ntered, and the automotive and fixed source emissions distributions are initialized for future use. . Model initialization is now complete, and PRINTl is called to print the initial concentration distribution. The HAIN program now enters the "computational loop," in which the governing partial differential equations are integrated numerically. Each pass through this loop results in advancing the simulaticn time by the value of one time step (1 to 5 minu~es). ;Before a forward step in time can be made, however, photochemical, met~orological, and source .variables must be updated. First, chemical reaction rate constants are calculated. (See Section VI-A in Appendix B for a discussion of the temporal variation of the photolysis rate constants.) Next, subroutine METEOR is called to update the meteorological variables. If the current- ly stored data are applicable at the time of the caJ.l, then control immediately returns to ~~IN. Recall that the meteorological variables are held constant for a period of one hour. If the data are no longer current, then the meteorological data tap~ is read, and new values of the meteorological variables are calculated. Upon return.from lI,ETEOR, subroutine SOURCE is entered, and new values of the automotive and fixed source emissions are computed. If the aircraft emissions rates stored in the temporary disc area are no longer applicable, then the next set of emissions is transferred from the tape to the disc area. This completes the definition of the photochemical,. meteorological, and source variables. We are now ready to proceed with the numerical integration of the govern- ing air shed equations. The numerical integration takes place in three sequential steps. In step I we move from column to column solving the p-direction convection- diffusion equation. The procedure followed in the solution of one column of points (or nodes) will be described. Power plant and aircraft emissions into the column are established first. The turbulent eddy diffusivity at each point in the column is then calculated in DIFCOF. Since Step I involves the iterative solution of implicit difference. equations, the computation primarily involves solving systems of linear equations of the form Ax == Y where A x y = tridiagonal matrix = solution vector = kno\.m vector F-9 ------- tie thus compute A and y, and then call SIMDIC: . which returns x. Each component of the solution vector represents th9' concentration' at one point in the colunm. ~nen the .iteration has converged at a given column, we repeat the calculation for the n~xt column. Step I is complete when these computations have been carried out for all columns. . Before Step II can be initiated, subroutine BCCONC must be called to provide pollutant concentration levels at points of horizontal inflow on east and west boundaries. upon return from BCCONC, Step II proceeds as described in Appendix D. When Step II is complete, BCCONCis again called,'Dut this time to supply con- centration levels for the north and sou~h boundaries. The compu- tation of Step III then follows. . At th~ conq~.usion of Step III, the numerical integration. procedure is compldted for a single t~e step. ..' Subroutine AVCONC is now entered, and the hourly-averaged ground- level .concentration array is incremented by the' ne\'l ground-level concentrations. If an hour has elapsed, then the array is divided by the number of time steps taken over that hour. AVCONC then pro- ceeds to print an average concentration map for each chemical species. Instantaneous ground-level concentration maps and vertical concentration profiles are printed at regular (say, every 20 minutes). The,program checks to see if be called to perform the output task. user-selected time intervals PlUNTl should We have no,'! made one .complete pass through the "coIi\putational loop." If the simulation is to be continued, then control passes to the beginning of the "loop," and the process described above is repeated. When the simulation is terminated, the entire concentra- . t~.on field .is punched onto cards~ This enables the user to re-start a simulation at some later time. F-lO ------- FIGURE F-2. Calculate Operating Parameters i f.1ETEOR FLOW DIAGRAM OF THE ATMOSPHERIC POLLUTION SIMULATION PRQGRN1 PRINTl Calculate Reaction .Rate Constants 0-- '-'''''~ . F-ll t40ve to Next Column. ., OC>I"'. . Input Aircraft Emissions (Tape) Output Air9raft . 'Emissions . (Disc) ------- .~ Set Up Linear Systems SnmIG No No F-12 Step II - Integrate in t-Dir. Step III - Integrate in T\-Dir. No No e ------- Stop '" Start AVCONC Increment Ave. Conc. Array Yes . Compute Average Cone. Output Ave. Cone. Start BCCONC 1 . procflss Inflow Cone. Return MAIN Start DIFCOF l- calCUlatle Vertical Diffusivity. 'Return F-13 ------- Return !-IAIN Reter." ) MAIN. . -~-- \.F-14 r;:-k. : Yes .~_. (0-- "j-----------~ I Cubic I : Spline I : Procedure: I Fol1m~s : L - -.- - _. - - - - - - Move to Next Colu:,\n Calculate Coefficient Matrix & Kno.,rn Vector SHIDIG Assign Cone. to New Grid. -~ ------- & 0-- Calculate Wind Components in Lm'ler t'lind Field Calculate. Vertical t1ind . Component PLUME F-15 start PLUME --- ...-- Advance to Next Power PUn t Calculate Amount & Spatial Dist. of Emissions No Return r-~ETEOR S,tart PRINTl ------- .C? Output vertical Profiles ~- Return HAIN Compute Temporal Distributions of Autos star;' SIMDI0 >J.I Conpute Auto Emissions 1 Compute Fixed Source Emissions Solve System of Linear Equations Return I4AIN ( R:urn) F-16 ------- ./ Sf~ Advance to Next Cell . Calculate Air Inputs to Cell L 'if . calCUlate] Air Outlet Velocities from Cell Return UETEOR ( F-17 no ------- II. THE AIRCRAFT EMISSIONS PROGRAM . The Atrnospheric Pollution Simulation Program, described in Section II reads aircraft emissions from a magnetic tape created by the Aircraft Emissions. Program. The a.ircraft emissions model .consists of two major p~rtsi ground operations and airborne operations. Emissions from these operations are treated as lumped volume sources, generated in the cell into which they are injected. We refer the reader to Section II of Appendix A for a discussion of the gove~ning model. .. The Aircraft subroutine AIRI. flight and ground ferred from MAIN. Emissions Program consists of a MAIN program and a , Subroutine AIRl calculates emissions for both the operation modes from aircraft and airport data trans- Input to the program:consists ofl . terrain elevations emissions rates tabulated accord~ngto aircraft class, opr.rating,mode,and chemical species the time spent by each class of aircraft in each operating mode individ~al airport data, including location of the airport on the grid and ground and flight operation characteristics meteorological data, such as inversion base elevations . . . The program uses this information to calculate aircraft emissions into those cells in which airplanes operate. The output of the program, i.e., the pollutant emissions rates'into each cell of the airshed, is recorded on the aircraft emissions tape for use by the Atmospheric Pollution Simulation Program. . Program operation starts in MAIN with the input of the first four items listed in the preceding paragraph, all of which are contained on p~nched cards. (See the schematic flow chart in Figure F-3.) We then begin the calculation of the aircraft emissions ~or each hour or, for that matter, for any other time interval of the simulated day. The emissions into each cell are stored in a large a.rray, which must be initialized to zero. The inversion base elevations are obtained from the meteorological data tape, and cell heights are calculated from (inversion base elevation - terrain elevation)/lO, where 10 is equal to the number of horizontal strata. We now calculate the emissions from each airport ,due to both ground and flight operations. Pertinent infor- mation about each mode of operation (e.g., take-off, taxi, landing, etc.) is transferred to AI~l, and the source array is incremented appropriately. After the emissions from each airport have been computed, the time of day and the source values are written on the magnetic tape. If fur tiler cal- culations are to be made, then the source array is again initialized, and the above sequence takes place for another time interval. F-IS ------- / It should be noted that this program can also function as a subroutine in the Atmospheric Pollution Simulation Program. We have operated the aircraft program separately to minimize the storage requirements for the Atmospheric Pollution Sinrolation Program. In future program developmentI efforts, we will investigate further the. question of integrating the Aircraft Emissions Program with the main simulation program. F-19 ------- FIGURE F-3. Start l-IAIN \!I. Input Airport Data Initialize Source Array to ro THE AIRCRAFT EMISSIONS PROGRM1 o F-20 (0 Calculate Cell Heights 1B Advan~e to Next Airport . l AJ:RL - Landing Mode o ------- (0 ! AIRl ~~~'p;Jt -. Time (Tape) J Output Emissions (Tape) ) eSTOP Ito Q ~v ~ V ~.NO ~~-- F-21 ( Start AIRl ) Yes calc~late J . Approach/ Take-off Emissions Ret urn MAIN G Calculate Ground Ope ra ti ons Emissions Return MAIN ------- . REFERENCES von Rosenberg, D. U.., "Methods for the Numerical Solution of Partial Differential Equations, II American .Elsevier Publishing' Company, NeW York (1969).' . . . J F-22 . . " ------- |