Homework on CMAQ

[ Home | Homework | ]


Homework Assignments


Due by
Details


23 Nov 2021

CMAQ installation, Part 1.  Lead by our TA, Tim Chui. 

Note that we will be using our laptops as terminals, for installation of CMAQ on the UBC/EOAS Optimum computer cluster.  CMAQ is too large to run on laptop computers.  We will provide accounts for you on Optimum. 

Don't forget that your CALPUFF homework on Slugs is due today.  Write your own code for this in your favorite computer language.

Bring your laptop. Make sure it is fully charged.

25 Nov 2021

CMAQ installation, Part 2. Lead by our TA, Tim Chui. 

Meanwhile, get started on your Smog-reactions HW, that is due next week.  Write your own code for this in your favorite computer language.

Bring your laptop.

30 Nov 2021


CMAQ installation, Part 2. Lead by our TA, Tim Chui.

Also, your homework on smog reactions is due today. (see below) 

Bring your laptop
.
30 Nov 2021
# ATSC 595D Atmospheric Dispersion Modeling -
# Homework exercise on smog reactions. Simplified from CMAQ.
# R. Stull - 5 Aug 2018, updated 26 Nov 2021.

# Adopted from Chapter 8 (section 8.2 Tropospheric Ozone) of free online textbook:
# "Atmospheric Chemistry" by István Lagzi, Róbert Mészáros, Györgyi Gelybó, Ádám Leelőssy
# "Copyright © 2013 Eötvös Loránd University, Budapest, Hungary.
#   This book is freely available for research and educational purposes."
# https://ttk.elte.hu/dstore/document/848/book.pdf
# (Accessed on 26 Nov 2021.)

# ========= Background Info ============

# Scenario: Simple example of ozone production from carbon monoxide and sunlight.
# Need to forecast concentrations of 7 chemicals that react with each other:
#   OH = hydroxyl radical
#   HO2 = hydroperoxly radical
#   NO = nitric oxide
#   NO2 = nitrogen dioxide
#   CO = carbon monoxide (a primary pollutant)
#   O1D = excited atomic oxygen (higher energy; more reactive)
#   O3 = ozone (as a secondary pollutant)

# Other 4 chemicals that are involved, but for which we don't need to forecast
#   O2 = molecular oxygen  (because ambient concentration is so large)
#   CO2 = carbon dioxide  (secondary product we don't care about here)
#   H2O = water vapour  (because ambient concentration is so large)
#   M = any non-reactive 3rd molecule, needed to remove energy from a reaction

# The 6 chemical reactions {as identified with CMAQ CB6 eq numbers}:
# {01'}  NO2 + h_nu     -> NO + O3
# {03}   O3 + NO        -> (O2) + NO2
# {09}   O3 + h_nu      -> (O2) + O1D
# {11}   O1D + H2O      -> 2 OH
# {25}   HO2 + NO       -> OH + NO2
# {123}  OH +CO +(O2)   -> HO2 + (CO2)
# where the items in parenthesis ( ) are items we can ignore in the rate eqs.
#(Note: {01'} is the combination of CB6 reactions {01} and {02}. )

# By neglecting the chemicals in ( ), this has the effect of reducing the
# order of the reaction. Also, we will assume h·nu is a known constant.
# The following are essentially 1st order reactions: {01, 09}.
# The following are essentially 2nd order reactions: {03, 11, 25, 123}.

# Let [ ] be the concentration of a chemical in units of molecules/cm^3.
# Rate equations predict the change in concentration [ ] of chemicals
# participating in a chemical reaction.  For example, consider:
#   1st order photolysis of A -> C.     Rate eq is:  d[C]/dt = j*[A]
#   2nd order reaction of A + B -> C.   Rate eq is:  d[C]/dt = k*[A]*[B]
# where
#   "photolysis rate constant" j has units of (s^-1)
#   "reaction rate constant"   k has units of ( (s * molecules/cm^3)^-1 )
# Although j varies with solar elevation, cloudiness, etc.,
# we will assume j = constant for simplicity in this exercise.

# For example, for reaction {25}, you have 1 production/loss term R
# and 4 rate eqs:
# {25.0}   Define R25 = k25*[HO2]*[NO]
# {25.1}   d[OH]/dt = R25
# {25.2}   d[NO2]/dt = R25
# {25.3}   d[HO2]/dt = -R25
# {25.4}   d[NO]/dt = -R25
# Note the minus sign for the last 2 eqs, which describe the loss of the
# reactants from the left side of reaction {25}.

# We can change each rate eq to a forward finite difference.  For example:
# For example, rate eq {25.1}    d[OH]/dt = R25  becomes
#
# [OH](new) = [OH](old) + delt * R25(old)
#
# where R25(old) = k25*[HO2]*[NO] (old) , and where delt is the timestep.
# Thus, each of the {25.x} eqs has a corresponding finite-difference form.


# ============== Your homework assignment ===========
# CAUTION: The notation both above and below use meta-code to specify
# physical processes and exercises.  The syntax might be incorrect
# for your computer language, so you will need to fix it.

#    For exercises A & B, write or type your answer on a piece
# of paper to submit.  (You don't need to use an eq. editor for this.)

# A(1). For each chemical reaction {} listed above, how many eqs apply,
#       and what is the resulting total number of equations?
#
#   For example,

#      Reaction {25} has 4 rate eqs + 1 R definition = 5 eqs.
#   Hint: don't write rate eqs to forecast any chemical in parenthesis (),
#   or for any chemical that we are considering to be constant (O2, H2O).
#
# A(2). For each one of the 7 chemicals that we want to forecast,

#   how many R production/loss terms apply?
#   For example,
#     NO has 3  R terms: 1 gain in {01} and 2 losses in {03} & {25}.


# B. Write all of the reaction eqs (similar to {25.0 to 25.4} shown above).
# But don't include reaction eqs for items in parenthesis ( ), or for
# any variable we assume is constant.


# C. In your favorite programming language, write into your computer
# program the finite difference form of the forecast eqs for each
# of the 7 chemicals of interest.
#    (1) CAUTION: During any one timestep, be sure to use all current
# concentrations to calculate ALL the R terms for all rate eqs.
#    (2) Note that some chemicals will have several R terms. Only
# after finding ALL the R terms will you use them in the
# forward difference with +sign for gains and -signs for losses.
# For example (using c_ to represent concentrations from now on):
#
# c_O3 = c_O3 + delt * [R02 - R03 - R09]
#
#   (3) CAUTION: truncate any negative concentrations to zero.


# D. Iterate your eqs forward in time on the computer, using the
# initial conditions IC and specifications given below.
#   Hint: use the IC values in all your finite difference eqs to step
# forward in time one timestep "delt" to get new concentration values
# of all 7 chemicals.
#    Then, once you have NEW concentration values for all 7 chemicals,
# (truncated to be non-negative)
# copy them into your array of OLD concentrations, use them to step
# forward another time step.  Repeat for each time step.
#    End at the elapsed time specified below.
# Hint: Don't save every timestep for plotting and other calculatoins. 
# Instead, save concentration values periodically (as specified)
# for all 7 chemicals, for the next 2 exercises.


# E(a). Using your output saved from D, convert the concentrations
# of the following 4 variables from c (molecules/cm^3) to mixing ratio
#  q (ppm) using the following approximate formula:
#
#         q(ppm) = c(molecules/cm^3) * (4.0e-14)
#
# E(b). Plot on the same one graph the following 4 concentration
# curves in units of ppm:  ppm_CO, ppm_NO, ppm_NO2, ppm_O3
# versus time in hours.  (Plot them only at intervals as given below
# in the specs. (Hint, since you saved values only at specified
# intervals, then these are the only values you need to plot.)


# F. Discuss the significance of your results from E.

# Please submit a printed copy of your results and your code.

# Note: CMAQ actually keeps track of over 200 chemical reactions,
#  resulting in over 1000 rate equations. We looked at only 6
#  of those reactions in this simplified exercise.


# ======== Specifications for your HW =========

# Rate constants for each reaction (slightly modified from cb6 rates for this exercise)
j01 = 1.0e-3     # (units 1/s)                      for {01'}  NO2 + h·nu   -> NO + O3
k03 = 1.73e-15   # (units (s · molecules/cm^3)^-1 ) for {03}   O3 + NO      -> (O2) + NO2
j09 = 1.0e-6     # (units 1/s)                      for {09}   O3 + h·nu    -> (O2) + O1D
k11 = 2.14e-10   # (units (s · molecules/cm^3)^-1 ) for {11}   O1D + H2O    -> 2 OH
k25 = 8.54e-12   # (units (s · molecules/cm^3)^-1 ) for {25}   HO2 + NO     -> OH + NO2
k123 = 2.28e-14  # (units (s · molecules/cm^3)^-1 ) for {123}  OH +CO +(O2) -> HO2 + (CO2)

# Time step info
delt_s = 0.0001     # time step (s) for your iterations
tend_h = 2.0        # duration of forecast (hours)
tsave_m = 2.0       # how often (minutes) to save results to plot

# Initial concentrations.
c_CO = 2.5e12     # molecules / cm^3
c_HO2 = 0.0       # molecules / cm^3
c_NO = 1.25e12    # molecules / cm^3
c_NO2 = 1.25e11   # molecules / cm^3
c_O1D = 0.0       # molecules / cm^3
c_O3 = 0.0        # molecules / cm^3
c_OH = 0.0        # molecules / cm^3

# Assume the following value is constant.
# (you don't need to forecast it, but you need to use it on the left side of the eqs)
c_H2O = 2.5e15    # molecules / cm^3  which is 0.01% of the air at sea level


# ====== End of assignment.  Good Luck. =========





Write this in your favorite computer language. 
Turn in printed copies of your output graphs and your code.

2 Dec 2021 Each student gives a 5-10 minute presentation to the class on one of the topics listed below, based on the info in the CMAQ overview website  https://www.epa.gov/cmaq/overview-science-processes-cmaq and on links it has to associated technical details. 
  1. Chemical: Photolysis chemistry
  2. Chemical: Multiphase chemistry
  3. Chemical: Airborne particle microphysics
  4. Chemical: Cloud-borne chemical processes
  5. Air-Surface Exchange: Emission processes (pick any 2)
  6. Air-Surface Exchange: Deposition processes
  7. Meteorological: Planetary boundary layer and surface layer processes
  8. Meteorological: Cloud and Radiation Processes
.
7 Dec 2021
# ATSC 595D Atmospheric Dispersion Modeling -
# Homework exercise on 1-D advection, demonstrating PPM method from CMAQ.
# PPM is Piecewise Parabolic Method.  For details, see:
# Chapter 7 of the CMAQ Science Document:
#      https://www.cmascenter.org/cmaq/science_documentation/pdf/ch07.pdf
# and the journal paper by Colella and Woodward, 1984, J. Computational Physics,
#      vol 54, pages174-201.
# But the best documentation is the actual fortran code for the
#      horizontal advection subroutine "hadv" from
#      https://github.com/USEPA/CMAQ/archive/5.2.1.zip and look in this path:
#      CCTM/src/hadv/yamo/hppm.F .  This version was updated on March 2018.
# R. Stull - 18 Aug 2018


# In this HW, a 1-D pollutant puff being advected in the x-direction by a constant
# wind u.  The two schemes you will compare are
#   (a) FTBS - Forward in time, Backward in space;
#   (b) RK3 - Runga-Kutte 3rd order, and
#   (c) PPM - Piecewise Parabolic Method (from CMAQ).


# Create the grid and initial conditions
imax = 1000             # number of grid points in x-direction
delx = 100.             # horizontal grid spacing (m)
delt = 10.              # time increment (s)
u = 5.                  # horizontal wind speed (m/s)

# 1) Calculate and display the Courant number.  You can display it in the
# graph in the question (3) if you would like.

# 2) (a) Create initial concentration distribution in the x-direction
conc <- rep(0.0, imax)   # initial concentration of background is zero
cmax = 10.0             # max initial concentration
conc[100:150] <- seq(0., cmax, len = 51)    # insert left side of triangle
conc[150:200] <- seq(cmax, 0., len = 51)    # insert right side of triangle
conc[20:40] <- seq(0., -0.5*cmax, len = 21)    # insert left side of triangle
conc[40:60] <- seq(-0.5*cmax, 0., len = 21)    # insert right side of triangle

#   (b) Plot the initial concentration distribution on a graph.

# 3) Also, on the same plot, show the ideal exact final solution, as given by
cideal <- rep(0.0, imax)   # initial concentration of ideal background is zero
cideal[800:850] <- seq(0., cmax, len = 51)    # insert left side of triangle
cideal[850:900] <- seq(cmax, 0., len = 51)    # insert right side of triangle
cideal[720:740] <- seq(0., -0.5*cmax, len = 21)    # insert left side of triangle
cideal[740:760] <- seq(-0.5*cmax, 0., len = 21)    # insert right side of triangle

# 4) Use "forward in time, backward in space" (FTBS) to advect concentration
# for the following number of time steps
nsteps = (imax - 300) / (u * delt / delx)
# and plot the resulting concentration on the same graph.

# 5) Repeat steps (2-4) to re-initialize, but plotting on a new graph, and using
# RK3 for the advection.  Use same number of time steps.

# 6) Repeat steps (2-4) to re-initialize, but plotting on a new graph, and using
# CMAQ PPM for the advection.  Use same number of time steps.
# Here is a code fragment that shows the PPM algorithm for our simplified
#
homework problem.  You can use or modify it to do this PPM exercise.

# 7) Discuss and compare the results of these three advection schemes.

Write this in your favorite computer language. 
Turn in printed copies of your output graphs and your code.

7 Dec 2021 Last Class day 

Presentation on how to visualize 3-D numerical simulations using the Vapor software.  

  • The period will start with Tim Chui leading us in the installation of this software. 
  • But most of the time will be a (Zoom) presentation by guest lecturer Dr. Nadya Moisseeva (Univ. of Hawaii), who has won awards for her visualization prowess using Vapor. 


Other Info





[ Home | Homework | ]

http://www.eos.ubc.ca/courses/atsc507/ADM/
Copyright © 2018, 2021 by Roland Stull.