.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/knote_ae_2015/plot_knoteae2015.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_knote_ae_2015_plot_knoteae2015.py: pykpp Tutorial ============== goals: Predict ozone at a monitor using HYSPLIT and pykpp authors: Barron Henderson, Robert Nedbor-Gross, Qian Shu # Overview 1. Create an idealized environment (Knote et al., 2015) 2. Create a simple trajectory model - Chemistry - CB05 - Emissions from NEI - Deposition from Knote et al., 2015 - Dilution following Knote et al., 2015 Assumes pykpp is installed python -m pip install http://github.com/barronh/pykpp/archive/main.zip Knote et al. 2015 ----------------- This example application is based on the Knote et al. (2015) simulation of a typical summer day with CB05. The example here uses the "Summer" PBL and temperature profile. At present, the example uses the a relatively "clean" initial environment (NOx(t=0) = 1ppb). Trajectory model ---------------- We'll use a KPP language to implement the model. Each piece will be saved as a string and then combined to produce an active model. - prepare an air pracel and model environment, - select gas-phase chemical reactions, - add emissions as chemical reactions, - add deposition as chemical reactions .. GENERATED FROM PYTHON SOURCE LINES 37-189 .. code-block:: default import time import numpy as np from pykpp.mech import Mech from io import StringIO import pandas as pd import matplotlib.pyplot as plt from pykpp.updaters import func_updater infile = StringIO(""" #INLINE PY_INIT TEMP = 298.15 P = 99600. t = TSTART = 0 * 3600. TEND = TSTART + 3600. * 24 DT = 60. MONITOR_DT = 3600 StartDate = 'datetime(2010, 7, 14)' Latitude_Degrees = 40. Longitude_Degrees = 0.00E+00 KDIL = 0. PBLH_OLD = 250. vd_CO = 0.03 vd_SO2 = 1.8 vd_NO = 0.016 vd_NO2 = 0.1 vd_NO3 = 0.1 vd_N2O5 = 4 vd_HONO = 2.2 vd_HNO3 = 4 vd_HNO4 = 4 vd_O3 = 0.4 vd_H2O2 = 0.5 vd_CH3OOH = 0.1 molps_to_molpcm2ps = Avogadro / 1200000**2 add_world_updater(func_updater(Update_M, incr = 360., verbose=False)) add_world_updater(func_updater(Update_THETA, incr = 360., verbose=False)) add_world_updater(interpolated_from_csv( 'summerenv.tsv', 'time', incr = 360., delimiter = '\\t', verbose=False )) add_world_updater(interpolated_from_csv( 'mean_emis.csv', 'time', incr = 360, verbose=False )) add_world_updater(func_updater( Monitor, incr = 7200., allowforce = False, verbose=False )) #ENDINLINE #INITVALUES CFACTOR = P * Avogadro / R / TEMP * centi **3 * nano {ppb-to-molecules/cm3} ALL_SPEC=1e-32*CFACTOR; TOTALNOx=1. M=1e9 O2=.21*M N2=.79*M H2O=0.01*M CH4=1750 CO=100 O3=30. NO = 0.1 * TOTALNOx NO2 = 0.9 * TOTALNOx SO2 = 1 N2O = 320 {B = 210.; what to do about B?} #MONITOR t/3600.; THETA; PBLH; TEMP; O3; #LOOKAT THETA; t/3600; PBLH; TEMP; O1D; OH; HO2; O3; NO; NO2; eNOx; TUV_J(6, THETA); emis_NO; emis_NO * molps_to_molpcm2ps / PBLH / 100 / CFACTOR * 3600; #include cb05cl.eqn EMISSION = ALD2 : emis_ALD2 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = ALDX : emis_ALDX * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = CH4 : emis_CH4 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = CL2 : emis_CL2 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = CO : emis_CO * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = ETH : emis_ETH * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = ETHA : emis_ETHA * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = ETOH : emis_ETOH * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = FORM : emis_FORM * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = HONO : emis_HONO * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = IOLE : emis_IOLE * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = ISOP : emis_ISOP * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = MEOH : emis_MEOH * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = NH3 : emis_NH3 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = NO + eNOx : emis_NO * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = NO2 + eNOx : emis_NO2 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = OLE : emis_OLE * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = PAR : emis_PAR * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = SESQ : emis_SESQ * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = SO2 : emis_SO2 * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = SULF : emis_SULF * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = TERP : emis_TERP * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = TOL : emis_TOL * molps_to_molpcm2ps / PBLH / 100 ; EMISSION = XYL : emis_XYL * molps_to_molpcm2ps / PBLH / 100 ; CO = DUMMY : vd_CO / PBLH / 100 ; SO2 = DUMMY : vd_SO2 / PBLH / 100 ; NO = DUMMY : vd_NO / PBLH / 100 ; NO2 = DUMMY : vd_NO2 / PBLH / 100 ; NO3 = DUMMY : vd_NO3 / PBLH / 100 ; N2O5 = DUMMY : vd_N2O5 / PBLH / 100 ; HONO = DUMMY : vd_HONO / PBLH / 100 ; HNO3 = DUMMY : vd_HNO3 / PBLH / 100 ; PNA = DUMMY : vd_HNO4 / PBLH / 100 ; O3 = DUMMY : vd_O3 / PBLH / 100 ; H2O2 = DUMMY : vd_H2O2 / PBLH / 100 ; MEPX = DUMMY : vd_CH3OOH / PBLH / 100 ; """) keywords = [ 'DUMMY', 'EMISSION', 'BNZHRXN', 'BNZNRXN', 'ISOPRXN', 'SESQRXN', 'SULRXN', 'TOLHRXN', 'TOLNRXN', 'TRPRXN', 'XYLHRXN', 'XYLNRXN' ] mech = Mech( infile, mechname='summercb05', incr=360, add_default_funcs=False, keywords=keywords ) # Create time start/stop matrix nhour = 24. start_end_ts = np.linspace( 0, nhour, int(nhour) * 12 + 1 ).repeat(2, 0)[1:-1].reshape(-1, 2) * 3600 # Capture initial values as "background" # concentrations mech.world['ybkg'] = mech.get_y(mech.parsed_world) # Create specific tolerance values for radical species # great care should be taken with O1D. # note that CMAQ EBI solver has no absolute tolerance and # uses 1 rtol for radicals and rxn counters. # atol = np.array([ 10 if (spc in ('O', 'O1D', 'HCO3', 'NTR') or spc[-3:] == 'RXN') else 1e-3 for spc in mech.allspcs ]) rtol = np.array([ 0.1 if (spc in ('O', 'O1D', 'HCO3', 'NTR') or spc[-3:] == 'RXN') else 1e-5 for spc in mech.allspcs ]) # Start timing the process runstart = time.time() # Update world for all interpolated values mech.world['t'] = start_end_ts.min() mech.update_y_from_world() mech.Update_World(forceupdate=True) .. rst-class:: sphx-glr-script-out .. code-block:: none Included /home/runner/work/pykpp/pykpp/pykpp/models/cb05cl.eqn /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring DUMMY warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring ISOPRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TRPRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring SULRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TOLNRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TOLHRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring XYLNRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring XYLHRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring BNZNRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring BNZHRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring SESQRXN warn('Ignoring %s' % spc) /home/runner/work/pykpp/pykpp/pykpp/parse.py:385: UserWarning: Ignoring EMISSION warn('Ignoring %s' % spc) Species: 65 Reactions: 223 {t:0,t/3600.:0,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:30} {t:0,t/3600.:0,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:30} .. GENERATED FROM PYTHON SOURCE LINES 190-196 Define a function to dilute with PBL ------------------------------------ - In a trajectory model, we need a PBL dilution - It should only dilute when PBL goes up - When PBL collapses, it should not concentrate - This will be called by the mechanism in between integrations .. GENERATED FROM PYTHON SOURCE LINES 196-222 .. code-block:: default def UpdatePBL(mech, world): """ Defining updater for PBL rise """ if 'y' not in world: return PBLH_OLD = world['PBLH_OLD'] ybkg = world['ybkg'] y = world['y'] PBLH = mech.world['PBLH'] if PBLH != PBLH_OLD: KDIL = np.minimum(1, PBLH_OLD / PBLH) mech.world['KDIL'] = KDIL world['PBLH_OLD'] = PBLH ydil = (KDIL * y) + (1 - KDIL) * ybkg nonzero = ybkg != 2.4195878568940516e-22 y[nonzero] = ydil[nonzero] else: mech.world['KDIL'] = 1 mech.update_world_from_y(y) # Add an updater that depends on the y vector mech.add_world_updater(func_updater(UpdatePBL, incr=360., verbose=False)) mech.world['PBLH_OLD'] = mech.world['PBLH'] .. GENERATED FROM PYTHON SOURCE LINES 223-226 Run Model --------- .. GENERATED FROM PYTHON SOURCE LINES 226-263 .. code-block:: default # Archive initial values mech.archive() # For each start/stop combination, update the world and integrate for t0, t1 in start_end_ts: # set current time to t0 t = t0 # Get y vector from world state y = mech.update_y_from_world() # integrate from t to t1 with initial state of y # ts, Y = mech.integrate( # t, t1, y0=y, solver='odeint', mxords=3, mxordn=3, atol=atol, # rtol=rtol, mxstep = 1000, hmax = 300, verbose=False # ) ts, Y = mech.integrate( t, t1, y0=y, solver='lsoda', atol=atol, rtol=rtol, nsteps=1000, max_step=300, max_order_ns=3, max_order_s=3, verbose=False ) # Update world to match new time mech.world['y'] = Y[-1] mech.update_world_from_y() mech.Update_World(forceupdate=True) # Update y vector for mixing of PBL # Set new time to last time of integration t = ts[-1] # Update world from y vector mech.update_world_from_y() mech.archive() # Optionally save to disk # mech.output() runend = time.time() print((runend - runstart), 'seconds') .. rst-class:: sphx-glr-script-out .. code-block:: none {t:3625,t/3600.:1,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:28} {t:7210,t/3600.:2,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:26} {t:7228,t/3600.:2,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:26} {t:10830,t/3600.:3,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:25} {t:14412,t/3600.:4,THETA:97,PBLH:2.5E+02,TEMP:2.9E+02,O3:23} {t:14434,t/3600.:4,THETA:97,PBLH:2.5E+02,TEMP:2.9E+02,O3:23} {t:18035,t/3600.:5,THETA:87,PBLH:3E+02,TEMP:2.9E+02,O3:23} {t:21612,t/3600.:6,THETA:76,PBLH:5.5E+02,TEMP:2.9E+02,O3:26} {t:21636,t/3600.:6,THETA:76,PBLH:5.5E+02,TEMP:2.9E+02,O3:26} {t:25238,t/3600.:7,THETA:65,PBLH:8.5E+02,TEMP:2.9E+02,O3:27} {t:28812,t/3600.:8,THETA:54,PBLH:1.1E+03,TEMP:3E+02,O3:28} {t:28839,t/3600.:8,THETA:54,PBLH:1.1E+03,TEMP:3E+02,O3:28} {t:32442,t/3600.:9,THETA:42,PBLH:1.3E+03,TEMP:3E+02,O3:28} {t:36032,t/3600.:10,THETA:32,PBLH:1.4E+03,TEMP:3E+02,O3:29} {t:36043,t/3600.:10,THETA:31,PBLH:1.4E+03,TEMP:3E+02,O3:29} {t:39644,t/3600.:11,THETA:22,PBLH:1.5E+03,TEMP:3E+02,O3:30} {t:43233,t/3600.:12,THETA:18,PBLH:1.5E+03,TEMP:3E+02,O3:33} {t:43246,t/3600.:12,THETA:18,PBLH:1.5E+03,TEMP:3E+02,O3:33} {t:46847,t/3600.:13,THETA:22,PBLH:1.5E+03,TEMP:3E+02,O3:35} {t:50435,t/3600.:14,THETA:31,PBLH:1.5E+03,TEMP:3E+02,O3:39} {t:50448,t/3600.:14,THETA:31,PBLH:1.5E+03,TEMP:3E+02,O3:39} {t:54052,t/3600.:15,THETA:42,PBLH:1.5E+03,TEMP:3E+02,O3:41} {t:57635,t/3600.:16,THETA:54,PBLH:1.5E+03,TEMP:3E+02,O3:44} {t:57655,t/3600.:16,THETA:54,PBLH:1.5E+03,TEMP:3E+02,O3:44} {t:61258,t/3600.:17,THETA:65,PBLH:1.5E+03,TEMP:3E+02,O3:45} {t:64836,t/3600.:18,THETA:76,PBLH:1.5E+03,TEMP:3E+02,O3:45} {t:64860,t/3600.:18,THETA:76,PBLH:1.5E+03,TEMP:3E+02,O3:45} {t:68465,t/3600.:19,THETA:87,PBLH:1.5E+03,TEMP:2.9E+02,O3:45} {t:72036,t/3600.:20,THETA:97,PBLH:1.5E+03,TEMP:2.9E+02,O3:45} {t:72066,t/3600.:20,THETA:97,PBLH:1.5E+03,TEMP:2.9E+02,O3:45} {t:75667,t/3600.:21,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:43} {t:79241,t/3600.:22,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:41} {t:79270,t/3600.:22,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:41} {t:82872,t/3600.:23,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:38} 15.545097827911377 seconds .. GENERATED FROM PYTHON SOURCE LINES 264-266 Make Plots ---------- .. GENERATED FROM PYTHON SOURCE LINES 266-289 .. code-block:: default data = mech.get_output() data['h_LST'] = data['t'] / 3600 data['NOx'] = data['NO'] + data['NO2'] data['jNO2*100'] = data['TUV_J(6,THETA)'] * 100 data['PBLH [km]'] = data['PBLH'] / 1e3 data['TEMP/25 [C]'] = (data['TEMP'] - 273.15) / 25. data['OH/10 [ppqv]'] = data['OH'] * 1e5 data['HO2 [pptv]'] = data['HO2'] * 1e3 fig, ax = plt.subplots() data.set_index('h_LST')[['emis_NO', 'jNO2*100', 'PBLH [km]', 'TEMP/25 [C]']].plot(ax=ax) fig.savefig('knote_physical.png') fig, ax = plt.subplots() data.set_index('h_LST')[['NOx', 'NO', 'NO2']].plot(ax=ax) ax.set(xlabel='hour [LST]', ylabel='ppb', ylim=(.1, 5), yscale='log') fig.savefig('chemical_nox.png') fig, ax = plt.subplots() data.set_index('h_LST')[['O3', 'OH/10 [ppqv]', 'HO2 [pptv]']].plot(ax=ax) ax.set(xlabel='hour (LST)') fig.savefig('chemical_ozone.png') .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_001.png :alt: plot knoteae2015 :srcset: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_002.png :alt: plot knoteae2015 :srcset: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_003.png :alt: plot knoteae2015 :srcset: /auto_examples/knote_ae_2015/images/sphx_glr_plot_knoteae2015_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 16.635 seconds) .. _sphx_glr_download_auto_examples_knote_ae_2015_plot_knoteae2015.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_knoteae2015.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_knoteae2015.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_