Siyuan Wang | Research - Ch.1(b) NOx at night

Ch 1(b). Nighttime NOx

OVERVIEW

- In the last chapter, we already set up the the model for the nighttime NOx chemistry. In this chapter, we'll keep working on it.
- We will look into the model infrastructure and start touching the source code from now on.
- We will also go through a few post-analysis tools, especially rate (budget) analysis, which is very very useful.

1. Hit CTRL + M

This brings you the main procedure window. Take a glance at it. Navigate yourself to a user-defined function called Sim1(). This is the major body of this model. See the schematic diagram.
There is a FOR loop, which loop thru each time step (at which you want the model to output the results). Within this loop, you can set a bunch of things, and then the solver is called and the ordinary differential equations are solved.
Before the FOR loop, the model firstly initialize, then set up temperature, pressure, RH, H2O, etc.
Then the model reads the initial concentrations and rate coefficients from the BIG TABLE. that's why if you provided these values there, you probably don't need to touch the source codes.
You may set a bunch of things inside the FOR loop, i.e. those time-dependent variables.

Now, in Sim1(), find this section. It's rough location is indicated by the big fat red arrow.

A couple of things:

TEMP is temperature (K), Pressure is the pressure (Pa), and RH is the relative humidity (%).
H2O is the water vapor concentration (molec/cm3) hard-coded for you. Note this this over liquid. If you work with ice, you need to touch this. Alternatively, if you have water vapor measurements or from other dataset, you can read it in. We'll talk about this later. RO2 is the total RO2 concentration (molec/cm3) - if you run mechanisms other than MCM, you probably don't need this.
KMT01, KMT02, ... are some of the complex rate coefficients from MCM, which are already loaded. If you don't run MCM, you probably don't need these either.
Jvalues are the photolysis frequencies (j-values, in s^-1). The parameterization from MCM is loaded for you, which depends on solar zenith angle (SZA). You may use other j-values parameterizations.

2. Scroll down & find this section...

just scroll a little lower. Again it's rough location is indicated by the arrow on the right.

You see this FOR loop. This loops thru each time step, i.e. the steps at which you asked the model to output the results. In the example we set up in the last chapter, there are 345 time steps. The model will call the solver for 345 times. In the ODE_ScratchBoard.ipf there is a function called ODE_Gen(kr, TimeScale, Conc, YDOT). This will be called many many times at each step until the system converges to a solution.
In the last chapter, we defined the initial concentrations and rate coefficients in the Big Table. These values are read into the model at the beginning of Sim1(). Can you find where?
Now, we can touch the source code to over-write the concentrations and / or rate coefficients. We can do this outside this FOR loop, so that these values will be set once. Alternatively, we can set these inside the FOR loop, and in this case these values will be set at each step. This is useful if you want some of your variables to be time-dependent. For example, you can let j-values to have some diurnal variation (will touch this in Chapter 4). Or you can read in measurements at each time step and constrain the model to these measurements at each step (Chapter 4 also).

3. Now copy-paste these two lines...

Concentration_Matrix[0][0] = (30) *Pressure*(6.0232E+8)/(8.314*TEMP) // NO Concentration_Matrix[0][4] = (40) *Pressure*(6.0232E+8)/(8.314*TEMP) // O3

Take a look at the flow: this is outside the FOR loop, but after reading the initial concentrations from the BIG TABLE. Therefore we over-write the initial concentrations of NO and O3 to 40 and 30 ppb. In this case the initial concentrations in the BIG TABLE don't affect the model.

I mostly setup the initial concentrations here within the main program to over-write whatever values has been read-in from the BIG TABLE. If this is the case, you don't need to set initial concentrations in the BIG TABLE.

4. Now scroll further down

...to the end of the MCM complex rate coefficient section, copy-paste these lines into the procedure window:

kr = 0 kr[0] = 1.4E-12*EXP(-1310/TEMP) // NO + O3 = NO2 kr[1] = 1.4E-13*EXP(-2470/TEMP) // NO2 + O3 = NO3 kr[2] = 1.8E-11*EXP(110/TEMP) // NO + NO3 = 2*NO2 kr[3] = KMT03 // NO2 + NO3 = N2O5 kr[3] = KMT04 // N2O5 = NO2 + NO3

Here we define the time-dependent rate coefficients (temperature-dependent), i.e. in case your temperature (or pressure, or RH) is time-dependent. Each element of the kr wave, kr[0] through kr[4], are loaded with new values at each step.
But have you noticed that i have a typo here? The last kr[3] should be kr[4]!!! In this case the model will not complain as the grammar is perfectly correct. The model runs perfectly fine but the results won't make sense. Gotta be careful! Sometimes bugs like this are very difficult to spot.
Now correct the typo for me, will you?

5. Rate (budget) analysis

Often times a species participates in many reactions. In this example, NO is involved in both R0 and R2. Which one is more important? To answer this question, we need to do some rate analysis.
There is a 2-dimensional wave called Rate_Matrix[m][n], in which the rate of reaction n (in molec/cm3/s) at time step m are stored.

Type in the command line:

display rate_matrix[][0] vs time_outputappendtograph rate_matrix[][2] vs time_output

As shown, in the first ~20 min, R0 is much faster than R2; after that R0 and R2 become equally important!
NO in this case is consumed only.

Now let's do the rate analysis for NO2! This is a bit more complicated, as this guy participated in a lot of reactions! It's removed in R1 and R3, and produced in R0, R2, and R4.
Let's plot all production pathways on one y-axis, and all removal pathways on another. Forgive my poor color-choice.
Note that you have to multiply R2 by a factor of 2, since two NO2 molecules are produced from this reaction (you don't have to do this if you're doing budget for NO or NO3)!
Noticed that R3 is pretty much equal to R4? What does this mean?

Imagine if you have hundreds or thousands of reactions! You can't possibly track all production and removal pathways by yourself unless you're a robot.
I have prepared this post-analysis tool to do that! Type in the command line:

Production_Removal("NO2")

and the program will look for all the removal and production pathways for you! You may customize this procedure for your own purposes (e.g., plot budget for other species, calculate branching ratios, ...)

6. This is it ;)

- I hope you are familiar with this model framework already. Go nuts!
- At this point the system is dark. We will have some lights in the next Chapter!

Google Sites

Report abuse