#pragma rtGlobals=3 // Use modern global access method and strict wave access. //// ************************************************************************************** //// ************************************************************************************** // The following code provides the necessary macros/functions for the calcuation of N2O5 // uptake coefficients and ClNO2 yields using the N2O5 (NOAA) Uptake box Model (NUM) as // described in McDuffie et. al., 2018 (Heterogeneous N2O5 uptake during winter: Aircraft // measurements during the 2015 WINTER campaign and critical evaluation of current parameterizations). // Instructions on how to setup and run the NUM are provided below. // //Last Uptades: 1/18/2018 // // //************************************************************************************************ //**********************************MODEL INSTRUCTIONS******************************************** //************************************************************************************************ //Step 1. Run Macro "NUM_ChemMechanism_Setup" //Step 2. Copy the reaction mechansim from Table 2 in the Main Text (or below) to the Reaction Mechanism Table //**Note, if you use your own mechanism or change the order, the reactions rate constants defined in //"N2O5UptakeModel_RunMultiStep" and its daughter functions must be updated accordingly for accurate calculations //Chemical Mechansims // Reactant1 Reactant2 Product1 Product2 Product3 //0. NO2 O3 NO3 //1. NO2 NO3 N2O5 //2. N2O5 NO3 NO2 //3. N2O5 HNO3 HNO3 //4. N2O5 HNO3 ClNO2 //5. NO3 RONO2 //6. NO3 NO2 //7. NO3 NO NO2 NO2 //8. NO O3 NO2 //9. NO2 NO //10. NO3 NO2 //11. N2O5 NO2 NO3 //12. ClNO2 NO2 //13. O3 Dummy //14. N2O5 Dummy //15. NO2 Dummy //16. O3 Dummy //17. NO3 Dummy //18. HNO3 Dummy //19. ClNO2 Dummy //20. NO Dummy //21. O3BKG O3 O3BKG //Step #3 - Run "N2O5UptakeModel_Setup" Macro - will average waves to specified interval, and choice of calculating kNO3 and dtEmission waves (or set to constant vaues) // **Note - to calculate kNO3_wave, must have VOC data in the current data folder - check the "VOCNames", "AFactor", and "Bfactor" waves to make sure kinetic data is correct //Step #4 - Run - "N2O5UptakeModel_RunMultiStep" - this is the main model loop // Results are printed to the "Uptake_ResultsData_MultiStep" wave //************************************************************************************************ //Constant declarations CONSTANT kb = 1.381e-23 //Boltzman constant [J/K] CONSTANT c = 2.998e4 //speed of light [cm/microseconds] //************************************************************************************************ //Clear macro menu to avoid clutter Menu "Macros" End Menu Menu "N2O5 Simulation Uptake Macros" "NUM_ChemMechanism_Setup" //Step #1 "N2O5UptakeModel_Setup" //Step #2 run this macro to setup appropriate wave dependencies, global variables, and instrument waves "N2O5UptakeModel_RunMultiStep" //Step #3, calculates N2O5 uptake, phi, ktot for every point. t0 = specified time before sunset. Pre sunset time is divided End Menu //************************************************************************************************ //************************************************************************************************ // The following Functions/Macros: ************************************************************ // // 1. N2O5UptakeModel_Setup() // 2. NNUM_Calc_AvgWavesFn() // 3. NUM_Calc_kNO3waveFn() // 4. NUM_Calc_EmissionswaveFn() // 5. NUM_ConvertNitrateFn() //// 6. NUM_Calc_GammaPhiFn() // 7. NUM_Calc_DataErrorFn() // 8. NUM_Calc_SumNO3ReactivityWaves() // 9. NUM_Replace8ZeroFn() // 10.NUM_Calc_kfactorWaveFn() // 11.NUM_Calc_AdjustedPlumeAgeFn() // // are used to setup the waves required for the box model and calculate uptake coefficients // and yields from the model output // // Author:Erin McDuffie (2018) //************************************************************************************************ //************************************************************************************************ //************************************************************************************************ //************************************************************************************************ Macro N2O5UptakeModel_Setup(T, P, NO2, O3, NOy, N2O5, ClNO2, SA, pNitrate, wtime) //Runs the appropriate functions and macros to make the waves required for the N2O5 (NOAA) Uptake Model (NUM) Macro //run KinSim "Machanism_Setup" first and use the mechansim provided in Table 2 of the Main Text String T = "RAF_ATX" //Ambient Temperature (C) String P = "RAF_PSXC" //Ambient Pressure (mbar) String NO2 = "TDLIF_NO2_ppbv" //Measured NO2 (ppbv) String O3 = "ACD_O3_ppbv" //Measured O3 (ppbv) String NOy = "ARNOLD_NOy_ppbv" //Measured NOy (ppbv) String N2O5 = "Arnold_N2O5_ppbv" //Measured N2O5 (ppbv) String ClNO2 = "CIMS1_UW_CLNO2_ppbv" //Measured ClNO2 (ppbv) String SA = "AMS_WetSA" //Total We Aerosol Surface Area (um2/cm3) String pNitrate = "AMS_Nitrate_NO2_lt_1um_AMS" //Measured Particle Nitrate (ug/sm3) String wtime = "Start_UTC" //Time wave Prompt T, "Temperature Wave (C)" Prompt P, "Pressure Wave (mbar)" Prompt NO2, "NO2 measurement wave (ppbv)" Prompt O3, "O3 measurement wave (ppbv)" Prompt NOy, "NOy measurment wave (ppbv)" Prompt N2O5, "N2O5 measurement wave (ppbv)" Prompt ClNO2, "ClNO2 measurement wave (ppbv)" Prompt SA, "Wet aeroosl surface area wave (um2/cm3)" Prompt pNitrate, "Aerosol Nitrate wave (ug/sm3)" Prompt wtime, "Time Wave" Variable/G vNDensity Variable Interval NUM_Mechanism_Compile("ChemKinetic", 1, 1) //Compile Mechanism, Print the mechanism to the command window Make/o/N = (1,4) ResultsData //make wave to store results data Make/o/N = 4 Iterations //make wave to store the number of iterations //make dependent waves for all species in the list of "Names" make/o/n = 50 Model_O3_ppbv Model_NO2_ppbv, Model_NO3_ppbv, Model_N2O5_ppbv, Model_HNO3_ppbv, Model_ClNO2_ppbv, Model_RONO2_ppbv, Model_NO_ppbv, Model_O3bkg_ppbv, Model_Cl_ppbv model_o3_ppbv := (o3/vndensity)*1e9 model_no2_ppbv := (no2/vndensity)*1e9 model_no3_ppbv := (no3/vndensity)*1e9 model_n2o5_ppbv := (n2o5/vndensity)*1e9 model_hno3_ppbv := (hno3/vndensity)*1e9 model_clno2_ppbv := (clno2/vndensity)*1e9 model_rono2_ppbv := (rono2/vndensity)*1e9 model_no_ppbv := (no/vndensity)*1e9 model_O3bkg_ppbv := (o3bkg/vndensity)*1e9 Make/o/n = (numpnts(T), 3) Uptake_ResultsData_MultiStep = nan //Make Final Results Wave //Columns in this wave: // 1: k(N2O5) // 2: gamma(N2O5) // 3: Phi(ClNO2) //Error Codes reported: // -1: before sunset // -2: Data Error (At least one wave does not have a value at the given point AND/OR NOx > Total Nitrate // -3: model did not converge in the first loop (N2O5) // -4: model did not converge in second loop (ClNO2) NUM_Time2Sunset() //calculates wave dtSunset (time since sunset, hours) NUM_Time2Sunrise() //calculates wave dtSunrise (time since sunrise, hours) //calc T and P Duplicate/o $T Amb_T Amb_T += 273 //temperature wave in K Duplicate/o $P Amb_P //pressure measurement in mbar NUM_ConvertNitrateFn($pNitrate) Duplicate/o $NOy SolubleNitrate_ppbv SolubleNitrate_ppbv += $pNitrate Interval = NUM_Calc_AvgWavesFn($wtime, $NO2, $O3, $N2O5, $ClNO2, $SA, SolubleNitrate_ppbv, Amb_T, Amb_P) //makes averaged/expaned waves for user-defined interval (new waves have original names ending in _exp) String Starttime = "Start_"+num2str(interval)+"s" String Stoptime = "Stop_"+num2str(interval)+"s" String Pexp = "Amb_P_exp" String Texp = "Amb_T_exp" String NO2exp = Nameofwave($NO2)+"_exp" String O3exp = Nameofwave($O3)+"_exp" String SNexp = "SolubleNitrate_ppbv_exp" NUM_Calc_kNO3waveFn($wtime, $starttime, $stoptime, $Texp, $Pexp) //Calculates "kNO3_wave_exp" NUM_Calc_EmissionswaveFn($NO2exp, $O3exp, $SNexp, $wtime, $starttime, $stoptime, $Pexp, $Texp, kNO3_wave_exp, Interval) //Calculates "dtEmissions_wave" End Macro //************************************************************************************************ Function NUM_Calc_AvgWavesFn(wtime, NO2, O3, N2O5, ClNO2, SA, SolubleNitrate, T, P) //calculate the average time waves to be used in each procedure Wave wtime Wave NO2 Wave O3 Wave N2O5 Wave ClNO2 Wave SA Wave SolubleNitrate Wave T Wave P Variable Interval Prompt Interval,"Average Interval (e.g. 10 for 10s)" DoPrompt "Enter Average Interval", Interval NUM_MakeAvgTimeWaveFn(wtime, interval) String Starttime = "Start_"+num2str(interval)+"s" Wave wStartTime = $Starttime String Stoptime = "Stop_"+num2str(interval)+"s" Wave wStoptime = $Stoptime String newnameNO2 = nameofwave(NO2)+"_exp" String newnameO3 = nameofwave(O3)+"_exp" String newnameN2O5 = nameofwave(N2O5)+"_exp" String newnameClNO2 = nameofwave(ClNO2)+"_exp" String newnameSA = nameofwave(SA)+"_exp" String newnameSN = nameofwave(SolubleNitrate)+"_exp" String newnameT = nameofwave(T)+"_exp" String newnameP = nameofwave(P)+"_exp" Make/o/n=(numpnts(wtime)) $newnameNO2, $newnameO3, $newnameN2O5, $newnameClNO2, $newnameSA, $newnameSN, $newnameT, $newnameP NUM_DoAveragUsingStartStopFn(wtime, NO2, wStarttime, wStopTime) Wave NO2Avg = $(NameofWave(NO2)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, NO2avg, wtime, $newnameNO2) NUM_DoAveragUsingStartStopFn(wtime, O3, wStarttime, wStopTime) Wave O3Avg = $(NameofWave(O3)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, O3avg, wtime, $newnameO3) NUM_DoAveragUsingStartStopFn(wtime, N2O5, wStarttime, wStopTime) Wave N2O5Avg = $(NameofWave(N2O5)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, N2O5avg, wtime, $newnameN2O5) NUM_DoAveragUsingStartStopFn(wtime, ClNO2, wStarttime, wStopTime) Wave ClNO2Avg = $(NameofWave(ClNO2)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, ClNO2avg, wtime, $newnameClNO2) NUM_DoAveragUsingStartStopFn(wtime, Sa, wStarttime, wStopTime) Wave SAAvg = $(NameofWave(SA)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, SAavg, wtime, $newnameSA) NUM_DoAveragUsingStartStopFn(wtime, SolubleNitrate, wStarttime, wStopTime) Wave SNAvg = $(NameofWave(SolubleNitrate)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, SNavg, wtime, $newnameSN) NUM_DoAveragUsingStartStopFn(wtime, T, wStarttime, wStopTime) Wave TAvg = $(NameofWave(T)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, Tavg, wtime, $newnameT) NUM_DoAveragUsingStartStopFn(wtime, P, wStarttime, wStopTime) Wave PAvg = $(NameofWave(P)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, Pavg, wtime, $newnameP) killwaves/Z NO2avg, O3avg, N2O5avg, ClNO2avg, SAavg, Tavg, Pavg, SNavg Return interval End Function //************************************************************************************************ Function NUM_Calc_kNO3waveFn(wtime, wstarttime, wstoptime, T, P) //calculate a new kNO3 wave from VOC data or constant user-defined value Wave wtime Wave wstarttime Wave wstoptime Wave T Wave P String kNO3 Prompt kNO3,"Enter kNO3. **Use constant value (= value in s-1) or calculated wave (= -1)?**" DoPrompt "Enter kNO3 value or choice to calculate new wave", kNO3 String cmd Make/o/n = (numpnts(wtime)) kNO3_wave, kNO3_wave_exp If(str2num(kNO3) == -1) Print "NO3 Reactivity Values:" cmd = "NUM_Calc_SumNO3ReactivityWaves(\"Amb_T_exp\", \"Amb_P_exp\")" Execute cmd Wave kNO3_SumTimewave kNO3_wave = kNO3_SumTimewave NUM_DoAveragUsingStartStopFn(wtime, kno3_wave, wStarttime, wStopTime) Wave kNO3Avg = $(NameofWave(kNO3_wave)+"_avg") NUM_ExpandItFillFn(wStarttime, wStoptime, kNO3avg, wtime, kNO3_wave_exp) Else kNO3_wave_exp = str2num(kNO3) Endif End Function //************************************************************************************************ Function NUM_Calc_EmissionswaveFn(wNO2, wO3, wNOy, wtime, wstarttime, wstoptime, wP, wT, wkNO3, Interval) //calculate a new dtEmissions wave from chemical data or constant user-defined value Wave wNO2 Wave wO3 Wave wNOy //soluble nitrate Wave wtime Wave wStarttime Wave wStoptime Wave wP Wave wT Wave wkNO3 Variable Interval String Emissions Prompt Emissions, "Enter Air Age. **Use constant value (= value in hrs) or calculated wave (= -1)?**" DoPrompt "Enter Age value or choose to calculate new wave", Emissions String cmd Make/o/n=(numpnts(wtime)) dtEmissions_wave, averagek_fitfactor_exp Wave dtSunset If(str2num(Emissions) == -1) Print "Calculating..." NUM_Calc_kFactorWaveFn(wNO2, wO3, wNOy, dtsunset, wT, wP, wkNO3, 1e-7, 0, interval) //calculates the stoichiometric factor for low N2O5 loss rate Wave w_FitFactor Duplicate/o w_FitFactor, lowk_FitFactor Print "Completed 1 of 2 Calculations" Print "Calculating..." NUM_Calc_kFactorWaveFn(wNO2, wO3, wNOy, dtsunset, wT, wP, wkNO3, 1, 0, interval) //calculates the stoichiometric factor for high N2O5 loss rate Duplicate/o w_FitFactor, highk_FitFactor averagek_fitfactor Print "Completed 2 of 2 Calculations" averagek_fitfactor += lowk_fitfactor averagek_fitfactor /= 2 NUM_FilterWaveFn(averagek_fitfactor, 1e-10, 1e10) //filter wave to remove non-calculated points (so it does impact the average calculation) NUM_DoAveragUsingStartStopFn(wtime, averagek_fitfactor, wStarttime, wStopTime) Wave averagek_fitfactor_avg NUM_ExpandItFillFn(wStarttime, wStoptime, averagek_fitfactor_avg, wtime, averagek_fitfactor_exp) NUM_Calc_AdjustedPlumeAgeFn(wO3, wNOy, wNO2, wP, wT, averagek_fitfactor_exp) Wave Adjusted_age_exp dtEmissions_wave = Adjusted_age_exp Else dtEmissions_wave = str2num(Emissions) Endif End Function //************************************************************************************************ Function NUM_ConvertNitrateFn(wnitrate) //convert aerosol nitrate measurement from ug/sm3 to ppbv Wave wNitrate wNitrate /= 62e6 wNitrate *= 6.022e23 wNitrate /= (100)^3 //now in molec/cm3 wNitrate /= ((1013*1e-4)/(1.381e-23*273)) //divide by standard ndensity wNitrate *= 1e9 End Function //************************************************************************************************ Function NUM_Calc_GammaPhiFn(SA, T, ktot, ClNO2, error) //This fn uses the modeled rate coefficients to calculate gamma, phi, and the total nitrate produced - prints results to "ResultsData" wave for reference in the top level macro Variable SA //aerosol surface area (um^2/cm^3) Variable T //Temperature (K) Variable ktot //total N2O5 loss rate - from model Variable ClNO2 //observed ClNO2 (ppbv) Variable Error //was there an error during the fit to ClNO2? Wave RateCoefficients //from odel Wave ResultsData //made in "Model SetUp" Macro Variable mN2O5 = 1.79e-25 //molecular mass of N2O5 in kg Variable v_mean = sqrt(8*kb*T/(pi*mN2O5)) //mean molecular speed of N2O5 (m/s) Variable kClNO2 = RateCoefficients[4] // production rate of ClNO2, from model Variable gamN2O5 //N2O5 uptake coefficient Variable R // ratio of ClNO2/total Nitrate Variable Phi // ClNO2 Yield Variable kHNO3 Variable bkgNSN //calculate background nocturnal soluble nitrogen (HNO3 + NO3-) If(numtype(SA) == 0) //Only calculate uptake coefficients when SA data is available gamN2O5 = 4*ktot/(v_mean*100*SA*1e-8) //N2O5 uptake coefficient ResultsData[0][1] = gamN2O5 Else ResultsData[0][1] = -2 //report error code Endif If(Error == 0 || RateCoefficients[4] !=1e-9) //if model converged for ClNO2, then calculate phi and bkg nitrate Phi = kClNO2/ktot //ratio of the two product pathways R = 1/((2/Phi)-1) //R = ClNO2/bkgNSN bkgNSN = ClNO2/R ResultsData[0][2] = phi Else bkgNSN = -4 ResultsData[0][2] = -4 Endif ResultsData[0][0] = ktot //always report ktot (at this point in the main model code, an error would be thrown before now if model didn't converge on N2O5 ResultsData[0][3] = bkgNSN //always report NSN End Function //************************************************************************************************ Function NUM_Calc_DataErrorFn(O3, NO2, N2O5,ClNO2, SA, TotalNitrate, Age, kNO3) //Function returns an error ( = 1) if at least one of the data waves does not have a numeric value at the given point, if calculated plume age is less than zero, or if no kNO3 data Variable O3, NO2, N2O5, ClNO2, SA,TotalNitrate, Age, kNO3 Variable DataError = 0 If(numtype(O3) == 2 || O3 < 0.0) DataError = 1 Endif If (numtype(NO2) == 2 || NO2 < 0.0) DataError =1 Endif If (numtype(N2O5) == 2 || N2O5 < 0) DataError =1 Endif If (numtype(ClNO2) == 2 || ClNO2 < 0) DataError =1 Endif If(Age < 0 || NumType(age) == 2) DataError = 1 Endif If(numtype(kNO3) == 2) DataError = 1 Endif Return DataError //return data error for reference in main model End Function //************************************************************************************************ Macro NUM_Calc_SumNO3ReactivityWaves(T, P) //calculate kNO3 reactivity wave, provided the VOC waves are available and names correctly String T //temeprature in K String P //pressure in mbar Variable factor //for scaling ppt and ppb into mixing ratio Silent 1 SetDataFolder root: NewDataFolder/o kNO3Data SetDataFolder root:kNO3Data Make/o/t/n = 56 VOC_wnames_orig VOC_wnames_orig[0, 14] = {"Methane", "CO", "C3H8", "CH2O", "HCN", "i_Butane", "n_Butane", "Acetaldehyde", "i_Butene_1_Butene", "Methyl_Bromide", "i_Pentane", "n_Pentane", "Methanol", "Ethanol", "Acrolein"} VOC_wnames_orig[15, 29] = {"Isoprene", "DMS", "CFC11", "CFC113", "Propanal", "Acetone", "2_Methylpentane", "3_Methylpentane", "n_Hexane", "Dichloromethane", "CH3CN", "MACR", "MVK", "3_MethylFuran", "Butanal"} VOC_wnames_orig[30, 44] = {"MEK", "Chloroform", "Carbon_Tetrachloride", "EthylNitrate", "PropylNitrate", "MBO", "Benzene", "n_Heptane", "Bromodichloromethane","Dibromomethane", "tert_ButylNitrate", "2nN_ButylNitrate", "2n3_PentylNitrate", "Toluene", "Ethylbenzene_mp_Xylene"} VOC_wnames_orig[45, 55] = {"o_Xylene", "Tetrachloroethylene", "Dibromochloromethane", "Chlorobenzene", "Bromoform", "APinene", "Camphene", "BPinene", "d_Limonene_3_Carene", "124_Trimethylbenzene", "123_Trimethylbenzene"} Make/o/n = 56 AFactor_orig AFactor_orig[0, 14] = {1e-18, 4e-19, 7e-17, 5.6e-16, 0, 3.05e-12, 2.76e-12, 1.4e-12, 3.14e-13, 0, 2.99e-12, 8.7e-17, 9.4e-13, 2e-15, 1.72e-13} AFactor_orig[15, 29] ={ 3.15e-12, 1.9e-13, 0, 0, 6.5e-15, 3e-17, 1.8e-16, 2.2e-16, 1.1e-16, 0, 0, 3.4e-15, 6e-16, 0, 1.7e-12} AFactor_orig[30, 44] = {0, 0, 0, 0, 0, 4.6e-14, 3e-17, 1.5e-16, 0, 0, 0, 0, 0,7e-17, 4.5e-16} AFactor_orig[45, 55] = {4.1e-16, 0, 0, 0, 0, 1.19e-12, 6.6e-13, 2.51e-12, 1.065e-11, 1.8e-15, 1.9e-15} Make/o/n = 56 BFactor_orig BFactor_orig[0, 14] = {0, 0, 0, 0, 0, -3060, -3279, -1860, -938, 0, -2927, 0, -2650, 0, -1190} BFactor_orig[15, 29] = {-450, 520, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1500} BFactor_orig[30, 44] = {0, 0, 0, 0, 0, -400, 0, 0, 0, 0, 0, 0, 0, 0, 0} BFactor_orig[45, 55] = {0, 0, 0, 0, 0, 490, 0, 0, 0, 0, 0} Duplicate/o Bfactor_orig wBfactor Make/o/t/n = 10 wVOC_Classes_orig wVOC_Classes_orig = {"alkanes", "alkenes", "halocarbons", "aromatics", "biogenics", "alcohols", "ketones", "aldehydes", "alkylnitrates", "unclassified"} MoveWave VOC_wnames_orig, root: MoveWave Afactor_orig, root: MoveWave BFactor_orig, root: MoveWave wVOC_Classes_orig, root: SetDataFolder root: Duplicate/o VOC_wnames_orig wVOC_wnames Duplicate/o AFactor_orig wAfactor Duplicate/o Bfactor_orig wBfactor Duplicate/o wVOC_Classes_orig wVOC_Classes killwaves VOC_wnames_orig, AFactor_orig, Bfactor_orig, wVOC_Classes_orig Variable Names = numpnts(wVOC_wnames) Variable i = 0 Variable pnumber = 0 String wname, new_wname, wname_nonan, wname_copy Duplicate/o $T Amb_ND Amb_ND = (1e-4*$P)/(kb*($T)) //make NO3 reactivity waves (kNO3 * VOC) i =0 Do If(i == 0 || i ==1) factor = 1e-9 Else factor = 1e-12 Endif wname = wVOC_wnames[i] new_wname = "kNO3_"+wname wname_nonan = new_wname + "_" Duplicate/o $wname $new_wname NUM_Replace8ZeroFn($new_wname) //change any -888s to zero (these values are below detection limit and should not be interpolated or used in average) NUM_InterpNaNFn($new_wname, wname_nonan) //interpret between points (when NaN's - i.e. no data reported) to increase data coverage $new_wname = $wname_nonan $new_wname = $wname_nonan*factor*amb_ND*wAFactor[i]*exp(wBFactor[i]/($T)) //calculate kNO3 for every point NUM_Filter2WavesFn($new_wname, dtSunset, 0, 1e100) //filter for times before sunset NUM_Filter2WavesFn($new_wname, dtSunrise, -1e100, 0) //filter for times after sunrise killwaves $wname_nonan Wavestats/q $new_wname If(v_numNaNs >= numpnts($new_wname)) //if a wave has no data (all NaNs), change it's nan's to zero so it doesn't ruin the addition below NUM_ReplaceNaNZeroFn($new_wname, wname_nonan) $new_wname = $wname_nonan killwaves $wname_nonan Endif i+=1 While(i < numpnts(wVOC_wnames)) Make/o/n = (numpnts(start_utc)) kNO3_SumAlkanes, kNO3_SumAlkenes, kNO3_SumHalocarbons, kNO3_Sumaromatics, kNO3_SumBiogenics, kNO3_sumAlcohols, kNO3_sumKetones, kNO3_sumAldehydes, kNO3_SumAlkylnitrates, kNO3_SumUnclassified, kNO3_SumTimewave kNO3_Sumalkanes = kNO3_C3H8 + kNO3_i_Butane +kNO3_n_Butane + kNO3_i_Pentane + kNO3_n_Pentane kNO3_sumalkanes += (kNO3_2_Methylpentane + kNO3_3_Methylpentane + kNO3_n_Hexane+ kNO3_n_Heptane) wavestats/q kno3_sumalkanes print "Alkanes - min:", v_min, "Max:", v_max kNO3_sumalkenes = kNO3_i_Butene_1_Butene wavestats/q kNO3_sumalkenes print "Alkenes - min:", v_min, "Max:", v_max kNO3_sumhalocarbons = kNO3_Methyl_Bromide + kNO3_Bromoform + kNO3_Dibromomethane + kNO3_Dibromochloromethane + kNO3_Bromodichloromethane kNO3_sumhalocarbons += (kNO3_Dichloromethane + kNO3_Chloroform + kNO3_Carbon_Tetrachloride + kNO3_Tetrachloroethylene) kNO3_sumaromatics = kNO3_Benzene + kNO3_Toluene + kNO3_Ethylbenzene_mp_Xylene + kNO3_o_Xylene + kNO3_Chlorobenzene + kNO3_124_Trimethylbenzene + kNO3_123_Trimethylbenzene wavestats/q kNO3_sumaromatics print "Aromatics - min:", v_min, "Max:", v_max kNO3_sumBiogenics = kNO3_Isoprene + kNO3_MVK + kNO3_MACR + kNO3_APinene + kNO3_BPinene + kNO3_d_Limonene_3_Carene kNO3_sumbiogenics += (kNO3_Camphene + kNO3_MBO + kNO3_DMS) wavestats/q kNO3_sumbiogenics print "Biogenics - min:", v_min, "Max:", v_max kNO3_sumalcohols = kNO3_Methanol + kNO3_Ethanol wavestats/q kNO3_sumalcohols print "Alcohols - min:", v_min, "Max:", v_max kNO3_sumKetones = kNO3_Acetone + kNO3_MEK wavestats/q kNO3_sumketones print "Ketones - min:", v_min, "Max:", v_max kNO3_sumAldehydes = kNO3_Acetaldehyde + kNO3_Propanal + kNO3_Acrolein + kNO3_Butanal + kNO3_CH2O wavestats/q kNO3_sumaldehydes print "Aldehydes - min:", v_min, "Max:", v_max kNO3_sumAlkylnitrates = kNO3_EthylNitrate + kNO3_PropylNitrate + kNO3_tert_ButylNitrate + kNO3_2nN_ButylNitrate + kNO3_2n3_PentylNitrate wavestats/q kNO3_sumalkylnitrates print "Alkylnitrates - min:", v_min, "Max:", v_max kNO3_sumunclassified = kNO3_HCN + kNO3_CH3CN + kNO3_3_MethylFuran + kNO3_CFC11 + kNO3_CFC113 kNO3_SumTimewave = kNO3_sumalkanes + kNO3_sumalkenes + kNO3_sumalcohols + kNO3_sumketones+kNO3_sumhalocarbons+kNO3_sumaromatics+kNO3_sumaldehydes +kNO3_sumalkylnitrates + kNO3_sumBiogenics +kNO3_sumunclassified +kNO3_methane+ kNO3_co wavestats/q kNO3_sumtimewave print "Sum - min:", v_min, "Max:", v_max End Function //************************************************************************************************ Function NUM_Replace8ZeroFn(wname) //replaces -888's in a wave with zeros wave wname variable i = 0 For(i =0 ; i < numpnts(wname); i+=1) If(wname[i] == -888) wname[i] = 0 Endif Endfor End Function //************************************************************************************************ Function NUM_Calc_kfactorWaveFn(wNO2, wO3, wSN, wdt, wT, wP, wkNO3, ktot, pnumber, interval) //Calculates the kfactor by fitting the exponential decay of NO2 for every (ith) point (after NO2, O3 are fit at the input ktot value) //Uses kNO3_wave_exp for the values of kNO3 at point i //Uses ambeint temperature and pressure at each point Wave wNO2 //NO2 wave Wave wO3 //O3 wave Wave wSN //Solube nitrate wave Wave wdt //Time since sunset (hours) Wave wT //T (K) Wave wP //P (mbar) Wave wkNO3 //kNO3 (s-1) Variable ktot //total N2O5 loss rate (s-1) Variable pnumber //start point number Variable interval //calculation interval Variable/g vNDensity Variable i, k Variable dtpresunset = 0 //time is that elapsed since sunset Variable SumIter Variable n = 0 Variable counter= 0 Variable Error = 0 Variable Model_RConstant Variable Calc_RConstant Wave InitialConcentrations //Model wave Wave RateCoefficients //Model Wave Wave KK //Model Wave Wave CalTime //Model integration time Wave Iterations Wave dtSunrise Wave dtSunset Wave NO2 //from ODE integration Wave O3 //from ODE integration Variable K0, K1 //fit parameters Wave w_coef make/o/n = 3 mycoefs //coefficient wave for fit Make/o/n=(numpnts(wP)) w_FitFactor, wDataError wDataError = NUM_Calc_DataErrorFn(wO3,wNO2,0,0,0,wSN,0,wkNO3) k = 0 InitialConcentrations = 0 RateCoefficients = 0 i = pnumber wavestats/Q dtsunset If(v_min > 0) dtsunset *= -1 Endif Do If(dtsunset[i] > 0 && dtSunrise[i] < 0 && wDataError[i] == 0) vNDensity = (1e-6*100*wP[i])/(kb*wT[i]) RateCoefficients[0] = 1.4e-13*exp(-2470/wT[i]) //T dependent rate coefficients for NO2, O3, NO3, N2O5 system RateCoefficients[1] = NUM_IUPAC_NO3NO2(wP[i], wT[i]) //JPLFalloff(P, T) // //NO2 + NO3 --> N2O5 based on 2nd order kinetic falloff regime for given P and T RateCoefficients[2] = NUM_IUPAC_UniN2O5(wT[i], wP[i]) //RateCoefficients[1]/JPLEq(T)// //N2O5 --> NO2 + NO3 assuming equilibrium RateCoefficients[3] = ktot RateCoefficients[4] = 0 RateCoefficients[5] = wkNO3[i] RateCoefficients[6] = 4.0e-12*0.3e8 //IUPAC HO2+NO3 --> NO2+OH * observed wintertime concentration from Stone ACP 2014 RateCoefficients[7] = 1.8e-11*exp(110/wT[i]) //IUPAC NO + NO3 --> 2NO2 RateCoefficients[8] = 2.07e-12*exp(-1400/wT[i]) //IUPAC NO + O3 --> NO2 KK = RateCoefficients If(InitialConcentrations[0] <=0) InitialConcentrations[0] = 1e-9*wNO2[i]*vNDensity*exp(RateCoefficients[0]*wO3[i]*1e-9*vNDensity*wdt[i]*3600) Endif If(InitialConcentrations[1] <= 0) InitialConcentrations[1] = 1e-9*wO3[i]*vNdensity Endif InitialConcentrations[3] = 0 Do NUM_FitNO2_MultiStep(wdt[i], dtpresunset, wNO2[i], 0, 1) NUM_FitO3_MultiStep(dtpresunset, wdt[i], wO3[i], 0) SumIter = Iterations[0]+Iterations[1] n+=1 While(SumIter > 2 &&n < 50 ) If(n == 50) Error = 1 w_FitFactor[i] = -5 Endif //Determine fit factor with derived initial concentrations //NO2 + O3 --> NO3 //dNO2/dt = -k*NO2*O3 = -k'*NO2 //[NO2] = [NO2]0 *eexp(-k'*t) --this is what we are fitting when fitting an exponential to modeled NO2 //k' = k*O3 = fit parameter //so to calculate the k value, divide by O3 (final or initial?) NUM_Model_MultiStep(dtpresunset, wdt[i], 0) K0 = 0 //Hold y0 constant at 0 K1 = NO2[0] //Hold A constant at initial NO2 concentration mycoefs = {K0, K1, 1e-5} //Print K0 //Print K1 CurveFit/Q/G/H="110" exp kwCWave=mycoefs NO2 /X=CalTime/D ///D//[1,49] //where y =y0 + Ae-kt // print K2 // print model_O3_ppbv[0]-model_o3_ppbv[49] Model_RConstant = mycoefs[2]/(O3[0]) //K2 = k (inv tau) Calc_RConstant = RateCoefficients[0] // print calc_rconstant w_FitFactor[i] = Model_RConstant/ Calc_RConstant Else w_FitFactor[i] = -2 Endif counter += 1 //print i i += interval While(i< numpnts(wNO2)) End Macro //************************************************************************************************ Function NUM_Calc_AdjustedPlumeAgeFn(wo3, wNOy, wNO2, wP, wT, wkfactor) Wave wO3 //O3 wave Wave wNOy //Soluble Nitrate Wave wNO2 //NO2 Wave wP //P Wave wT //T Wave wkfactor //fit factor wave Variable age_s, age_hr, ndensity, i, NO2, O3, NOy, k NVAR kb Make/o/n= (numpnts(wO3)) Adjusted_Age_exp Do ndensity = ((1e-6*100*wP[i])/(kb*wT[i])) k = 1.4e-13*exp(-2470/wT[i]) //T dependent rate coefficients for NO2, O3, NO3, N2O5 system NO2 = wNO2[i] *(1e-9*nDensity) //change to concentration O3 = wO3[i] * (1e-9*nDensity) NOy = wNOy[i] * (1e-9*nDensity) age_s = ln(NO2/(NOy))/(-wkfactor[i]*k*O3) age_hr = age_s/3600 Adjusted_age_exp[i] = age_hr i+=1 While(i 0) dtsunrise *=-1 //set to 1 for day into night flights, -1 for night to day flights Endif Make/o/n = (numpnts(start_utc)) DataError DataError = NUM_Calc_DataErrorFn($O3, $NO2, $N2O5, $ClNO2, $SA,$NOy, dtEmissions_wave, kNO3_wave_exp) //makes a wave with 1's where there are missing measruements in the times series RateCoefficients[14, 22] = Dilution //Set Dilution Rate for all species KK = RateCoefficients InitialConcentrations = 0 //MAIN MODEL LOOP*** Do If(dtSunset[pnumber] > 0 && dtSunrise[pnumber] < 0) //check to make sure: 1) the sun is set, 2) all waves have real, non-negative values at the given point (filtering for measurement errors) If(DataError[pnumber] == 0 && dtEmissions_wave[pnumber] > 0) //Step #0.1 vNDensity = ((1e-6*100*Amb_P[pnumber])/(kb*Amb_T[pnumber])) //calculation for ambient number density for given T and P wave - molec/cm3 InitialConcentrations[9] = BkgO3*1e-9*vNDensity //Set the concnetration of background O3 (use-input) RateCoefficients[4] = 0 //Make sure ClNO2 production rate = 0 KK[4] = 0 RateCoefficients[5] = kNO3_wave_exp[pnumber] //Set the NO3 + VOC rate coefficient for this point //Step #0.2 Set the simulation time, pre-sunset time, and NOx initiailization form (NO or NO2), depending on the time since emission (i.e. age, dtEmission) If(dtEmissions_wave[pnumber] > dtSunset[pnumber]) dtSimTime = dtSunSet[pnumber] //aged plume --> sim start is sunset If(dtEmissions_wave[pnumber] > (dtSunset[pnumber]+PreSS_Start)) //if older than designated pre-sunset start time, start at designated hours before sunset dtPreSunset = PreSS_Start // hours - derived previously from daytime N2O5 Steady State analysis NOxForm = 1 //model NOx starting as NO2 Else dtPreSunSet = dtemissions_wave[pnumber] - dtsunset[pnumber] //if younger than designated start time, but older than sunset, presunset is equal to difference NOxForm = 0 //model NOx start as NO Endif Else dtSimTime = dtEmissions_wave[pnumber] //fresh plume --> sim start is air age (time of emission) NOxForm = 0 //model NOx as NO dtpresunset = 0 //only model plume after sunset (i.e. after emission) Endif //Steps # 1 & 2 PauseUpdate Error = NUM_Fit_NO2_O3_N2O5_MultiStep(dtPreSunSet, dtSimTime, $NO2[pnumber], $O3[pnumber], $N2O5[pnumber], Amb_T[pnumber], Amb_P[pnumber], bkgO3, NOxForm) //iterate (1) initial NO2 and O3 and (2) N2O5 loss rate until model resproduces observations ResumeUpdate If(Error == 0) //If model converged in the previous step, continue to ClNO2 iteration ktot = RateCoefficients[3] //initially assign all loss of N2O5 to HNO3 channel for the start of the following ClNO2 iteration wInit_NO2[pnumber] = InitialConcentrations[0] //save the derived initial concentrations wInit_O3[pnumber] = InitialConcentrations[1] wInit_NO[pnumber] = InitialConcentrations[7] //Step #3 PauseUpdate Error = NUM_Fit_ClNO2_MultiStep(dtpresunset, dtSimTime, $ClNO2[pnumber], ktot, bkgo3) //iterate production rate of ClNO2 until model reproduces observations ResumeUpdate //Step #4 NUM_Calc_GammaPhiFn(Amb_SA[pnumber], Amb_T[pnumber], ktot, $ClNO2[pnumber], error) //uses model results to calculate the N2O5 uptake coefficient, CLNO2 yield, and ppbv of total nitrate produced (nocturnal soluble nitrate) If(Error == 0) //If model converged, save data ModelOutput_RONO2[pnumber] = Model_RoNO2_ppbv[49] //track final values of RoNO2, HNO3, and ClNO2 to track where the nitrogen is going each point ModelOutput_ClNO2[pnumber] = Model_ClNO2_ppbv[49] ModelOutput_N2O5[pnumber] = Model_N2O5_ppbv[49] ModelOutput_HNO3[pnumber] = Model_HNO3_ppbv[49] NO3_ppbv_Model[pnumber] = model_no3_ppbv[49] w_dtsimtime[pnumber] = dtsimtime +dtpresunset Uptake_ResultsData_Multistep[pnumber][2] = ResultsData[0][2] //Phi Else print "**** MODEL DID NOT CONVERGE (ClNO2) ****" //if model did not converge for ClNO2, report errors ModelOutput_RONO2[pnumber] = -4 //track final values of RoNO2, HNO3, and ClNO2 to track where the nitrogen is going each point ModelOutput_ClNO2[pnumber] = -4 ModelOutput_N2O5[pnumber] = -4 ModelOutput_HNO3[pnumber] = -4 NO3_ppbv_Model[pnumber] = -4 Uptake_ResultsData_Multistep[pnumber][2] = -4 //Phi Endif Uptake_ResultsData_Multistep[pnumber][0] = ResultsData[0][0] // regardless of whether ClNO2 converged, save ktot results Uptake_ResultsData_Multistep[pnumber][1] = ResultsData[0][1] // save gamma N2O5 NSN[pnumber] = ResultsData[0][3] // save NSN w_dtsimtime[pnumber] = dtsimtime +dtpresunset // log the simulation time RateCoefficients[3] = ktot //reset for next iteration KK[3] = RateCoefficients[3] //reset for next iteration RateCoefficients[4] = 0 //reset for next iteration KK[4] = RateCoefficients[4] //reset for next iteration Else print "**** MODEL DID NOT CONVERGE (N2O5) ****" //If the model did not converge for N2O5, report errors Uptake_ResultsData_MultiStep[pnumber][] = -3 //report error in results wave if the model does not converge wInit_NO2[pnumber] = -3 wInit_O3[pnumber] = -3 wInit_NO[pnumber] = -3 ModelOutput_RONO2[pnumber] = -3 //track final values of RoNO2, HNO3, and ClNO2 to track where the nitrogen is going each point ModelOutput_ClNO2[pnumber] = -3 ModelOutput_N2O5[pnumber] = -3 ModelOutput_HNO3[pnumber] = -3 NO3_ppbv_Model[pnumber] = -3 NSN[pnumber] = -3 Endif Else //print "Data Error" wInit_NO2[pnumber] = -2 wInit_O3[pnumber] = -2 wInit_NO[pnumber] = -2 ModelOutput_RONO2[pnumber] = -2 //track final values of RoNO2, HNO3, and ClNO2 to track where the nitrogen is going each point ModelOutput_ClNO2[pnumber] = -2 ModelOutput_N2O5[pnumber] = -2 ModelOutput_HNO3[pnumber] = -2 NO3_ppbv_Model[pnumber] = -2 NSN[pnumber] = -2 Uptake_ResultsData_multistep[pnumber][] = -2 Endif Else // print "before sunset" wInit_NO2[pnumber] = -1 wInit_O3[pnumber] = -1 wInit_NO[pnumber] = -1 ModelOutput_RONO2[pnumber] = -1 //track final values of RoNO2, HNO3, and ClNO2 to track where the nitrogen is going each point ModelOutput_ClNO2[pnumber] = -1 ModelOutput_N2O5[pnumber] = -1 ModelOutput_HNO3[pnumber] = -1 NO3_ppbv_Model[pnumber] = -1 NSN[pnumber] = -1 Uptake_ResultsData_multistep[pnumber][] = -1 Endif counter += interval //user-defined interval pnumber += interval print pnumber-1 While(pnumber < numpnts($O3) && counter < numpnts($O3)) //loop while both the counter and point current point number are < the total points in the timeseries End Macro //************************************************************************************************ //********************************END OF MAIN BOX MODEL MACRO************************************* //************************************************************************************************ //************************************************************************************************ Function NUM_Fit_NO2_O3_N2O5_MultiStep(dtpresunset, dt, NO2, O3, N2O5, T, P, BkgO3, NOxForm) //Function iteratively solves chemical ODEs to fit measured NO2, O3, and N2O5 Variable dtpresunset Variable dt //Time since simulation start (hours) - timestep in hours Variable NO2 //Measured NO2 (ppbv) Variable O3 //Measured O3 (ppbv) Variable N2O5 //Measured N2O5 (ppbv) Variable T //Ambient Temperature (K) Variable P //Ambient Pressure (mbar) Variable bkgO3 //ppbv Bkg O3 Variable NOxForm //Adjust initial NO (0) or NO2 (1) concentration Variable SumIter, n, Error, Stop String cmd Wave RateCoefficients, InitialConcentrations, CalTime, Iterations, ErrorWave, KK NVAR vndensity variable marker = 0 //Initial estimates for concentrations and rate constants RateCoefficients[0] = 1.4e-13*exp(-2470/T) //T dependent rate coefficients for NO2, O3, NO3, N2O5 system RateCoefficients[1] = NUM_IUPAC_NO3NO2(P, T) // NO2 + NO3 --> N2O5 based on 2nd order kinetic falloff regime for given P and T RateCoefficients[2] = NUM_IUPAC_UniN2O5(T, P) //RateCoefficients[1]/JPLEq(T) // //N2O5 --> NO2 + NO3 assuming equilibrium RateCoefficients[6] = 4e-12*0.3e8 //Stone, ACP 14 - NO3 + HO2 --> NO2 +OH + O2 (using IUPAC rate coefficient and observed concentration from Stone) RateCoefficients[7] = 1.8e-11*exp(110/T) //IUPAC NO + NO3 --> 2NO2 RateCoefficients[8] = 2.07e-12*exp(-1400/T) //IUPAC NO + O3 --> NO2 RateCoefficients[3] = 1e-4 //Initial guess for ktot RateCoefficients[4] = 0 //make sure no clno2 is being formed KK = RateCoefficients If(NOxForm == 0) InitialConcentrations[0] = 0 InitialConcentrations[7] = 1e-9*NO2*vNDensity*exp(RateCoefficients[0]*O3*1e-9*vNDensity*dt*3600) //Estimated initial NO - molec/cm3 Else InitialConcentrations[7] = 0 InitialConcentrations[0] = 1e-9*NO2*vNDensity*exp(RateCoefficients[0]*O3*1e-9*vNDensity*dt*3600) //Estimated initial NO2 - molec/cm3 Endif InitialConcentrations[1] = 1e-9*O3*vNdensity //initial guess at O3 concentation at t= 0 - molec/cm3 //Iterate until NO2, O3, and N2O5 are simultaneous fit to observations (i.e. SumIter = 3) Do NUM_FitNO2_MultiStep(dt, dtpresunset, NO2, BkgO3, NOxForm) NUM_FitO3_MultiStep(dtpresunset, dt, O3, BkgO3) NUM_FitN2O5_MultiStep(dtpresunset, dt, N2O5, BkgO3) SumIter = Iterations[0] + Iterations[1] + Iterations[2] //Iterations Wave Assignments: NO2 = 0, O3 = 1, N2O5 = 2, ClNO2 = 3 //Error handling n += 1 While(SumIter > 3 && n < 50) If(n == 50) //if the NO2, O3, N2O5 loop has tried 50 times and it wont converge, report error Error =1 Else Error = 0 Endif Return Error End Function //************************************************************************************************ Function NUM_FitNO2_MultiStep(dt, dtpresunset, NO2_t, bkgO3, NOxForm) //Solve ODEs until final modeled NO2 is within fit accuracy of initially input (downwind) NO2 value Variable dt //simulation duration Variable dtpresunset //time between start of chemistry and sunset Variable NO2_t //measured NO2 that the model is fit to Variable bkgO3 //ppbv of bkg O3 Variable NOxForm //NOx initialized as NO (=0), NOx initialized as NO2 (=1) Variable dNO2, fNO2 //aboslute (dNO2) and fractional (fNO2) difference between measured and calculated values for final NO2 Variable Stop = 0 // stop model (=1), keep going ( = 0) Variable n = 1 //iteration number counter Variable gNO2 = 0.005 //Fit threshold (maximum fractional difference) Variable initNO, initNO2 //tracker variables for initial conc. Variable slast, snew //variables that indicate what direction (+ or -) initial concentrations need to be adjusted (i.e. -1= initial concentration needs to increase, 1= initial conc. needs to decrease) Variable StepSize //size of iterative step change Variable RoundNum = 0 //round number Wave model_NO2_ppbv, InitialConcentrations, Iterations NVAR vndensity Duplicate/o InitialConcentrations InitConcWave //wave for logging the original initial concentrations of all species (O3 and NO or NO2 should be the only populated species) InitNO2 = InitConcWave[0] //log initial NO2 from previous iteration InitNO = InitConcWave[7] //log initial NO from previous iteration //Initial Model Run to gauge current difference InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, bkgo3) //Run Model dNO2 = model_NO2_ppbv[49]-NO2_t //Determine initial difference between model and measurements StepSize = dNO2 //Set Initial Ste size If(dNO2 > 0) //check to see if the model has predicted too much or too little NO2 slast = 1 Else slast = -1 Endif //Solve ODEs iteratively until final NO2 is within the fit accuracy of the initially input (downwind) NO2 value Do fNO2 = dNO2/NO2_t //fractional NO2 difference between model and meausurement If (abs(fNO2) > gNO2) //stop iterating if calculated value is within fit accuracy Stop = 0 If(NOxForm == 0) InitNO -= (StepSize)*1e-9*vNdensity //adjust initial NO to new value InitConcWave[7] = InitNO InitNO2 = 0 InitConcWave[0] = InitNO2 Else InitNO2 -= (StepSize)*1e-9*vNdensity //adjust initial NO2 to new value InitConcWave[0] = InitNO2 InitNO = 0 InitConcWave[7] = InitNO Endif If(InitNO2 < 0) //check to make sure that the step size didnt push init concentration below zero InitNO2 = 0 InitConcWave[0] = InitNO2 Endif If(InitNO <0) InitNO = 0 InitConcWave[7] = InitNO Endif InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, bkgo3) //re-run the model with adjsted initial concentrations dNO2 = model_NO2_ppbv[49]-NO2_t //calculate new absolute NO2 difference If(dNO2 > 0) //check to see if the model needs to increase or decrease the initial concentrations snew = 1 Else snew = -1 Endif If(snew != slast) //if the model has overshot in either direction, cut the step size in half and reverse the direction StepSize /= -2 slast = snew //update direction to new value Endif n+=1 If(InitNO ==0 && InitNO2 ==0 && snew ==-1) //keep going one more round if the init conc is zero and set to increase in the next round Stop = 0 Elseif(InitNO ==0 && InitNO2 ==0 && snew ==1) //stop if the init conc is zero and set to decrease in the next round - and report a non-1 n value Stop = 1 n = 20 Endif Else Stop = 1 Endif While (Stop == 0 && n <= 20) //continue until the model converges or has iterated 20 times InitialConcentrations = InitConcWave //update Initial Concentrations wave Iterations[0] = n //log the number of iterations - used in the higher level function End Function //************************************************************************************************ Function NUM_FitO3_MultiStep(dtpresunset, dt, O3_t, BkgO3) //Solves Chemcial ODEs, changing the initial concentration each iteration, until final O3 is within the specifed fit accuracy of initially input (downwind) O3 Variable dtpresunset //time between start of chemistry and sunset Variable dt //nocturnal simulation duration Variable O3_t //downwind (measured O3) Variable BkgO3 //ppbv bkg O3 Variable dO3, fO3 //absolute (dNO2) and fractional (fNO2) difference between input and calculated values for final O3 Variable Stop = 0 //Stop iterations (=1), keep going (= 0) Variable n = 1 //iteration number counter Variable gO3 = 0.005 //Fit threshold (maximum fractional difference) Variable InitO3 //tracker variable for initial concentration Variable snew, slast //variables that indicate what direction (+ or -) initial concentrations need to be adjusted (i.e. -1= initial concentration needs to increase, 1= initial conc. needs to decrease) Variable StepSize //size of iterative step Wave model_o3_ppbv, initialconcentrations, iterations Wave/T Names nvar vndensity Make/o/n=(numpnts(names)) InitConcWave //wave for logging the original initial concentrations InitConcWave = InitialConcentrations InitO3 = InitConcWave[1] //log the initial O3 concnetrations (it will be updated later) InitialConcentrations = InitConcWave //determine the initial difference between model and measurement NUM_Model_MultiStep(dtpresunset, dt, bkgo3) //Intial Model Run dO3= model_o3_ppbv[49] - O3_t //Initial difference between model output and measured O3 StepSize = dO3 //set the initial step size If(dO3 < 0) //determine whether the init O3 concentration needs to increase or decrease slast = -1 Else slast = 1 Endif //Solve ODEs iteratively until final O3 is within the fit accuracy of the initially input (downwind/measured) O3 Do fO3 = dO3/O3_t //fractional O3 difference If(abs(fO3) > gO3) //stop iterating if calculated value is within fit accuracy Stop = 0 InitO3 -= (StepSize)*1e-9*vndensity //adjust initial O3 value (will go up or down as necessary due to dO3 being the absolute difference) If(InitO3 < 0) //check to see if Init O3 was pushed below zero InitO3 = 100 Endif InitConcWave[1] = InitO3 InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, bkgO3) //re-run the model with adjusted O3 concentration dO3 = model_O3_ppbv[49] -O3_t //re-calculate absolute O3 difference If(dO3 < 0) //determine if init O3 needs to increase or decrease snew = -1 Else snew = 1 Endif If(snew != slast) //if the model has overshot (i.e. O3 needed to increase and after previous iteration it now needs to decrease), then cut the step size in half and change direction StepSize /= -2 slast = snew Endif n +=1 //count iteraiton number Else Stop = 1 Endif While (Stop ==0 && n <= 20) InitialConcentrations = InitConcWave //update Initial Concentrations to derived values from this procedure Iterations[1] = n //log iteration number End Function //************************************************************************************************ Function NUM_FitN2O5_MultiStep(dtpresunset, dt, N2O5_t, BKgO3) //Solves Chemical ODEs, changing the N2O5 loss rate each iteration, until final N2O5 is within the specified fit accuracy of initially input (downwind, measured) N2O5 //assuming that at time t = 0, N2O5 = 0 Variable dtpresunset //time between start of chemistry and sunset Variable dt //nocturnal simulation duration Variable N2O5_t //measured N2O5 that is being fit to Variable BkgO3 //pbv Bkg O3 Variable dN2O5, fN2O5 //absolute and fractional N2O5 differences between calculated and measured downwind values Variable Stop = 0 //Stop iterations (=1), keep going (= 0) Variable n = 0 //iteration counter variable variable gN2O5 = 0.01 //Fit threshold (maximum fractional difference) Variable snew, slast //variables that indicate what direction (+ or -) kN2O5 need to be adjusted (i.e. kN2o5 down = -1, kN2O5 up = +1) Variable StepSize //Iteration step size for kN2O5 Variable order, olast //tacker variables for kN2o5 order of magnitude for current and previous iterations Variable remainder, interger //variables to help calculate correct order of magnitude Variable counter = 0 //count how many iterations through the order of magnitude loop (if counter =1, then the order of magnitude is the same as the initial value and start with that kN2O5 for the remaining procedre) Variable InitCoef //variable to log the initial kN2O5 Variable CalcMag = 1 // = 1 if order of magnitude needs to be calculated, =0 if order does not need to be calculated (e.g. if previous iteration is far from solution, the order of magnitude needs to be calculated - meant to speed up calculation time) Wave InitialConcentrations, CalTime, KK, RateCoefficients, Iterations, model_n2o5_ppbv Make/o/n=(numpnts(InitialConcentrations)) InitConcWave //make wave to keep track of the initial concentrations of each compound InitConcWave = InitialConcentrations //log the initial concentrations from the past calculation (e.g. O3 and NO2 functions) so they can be reset every iteration (since it is only the N2O5 loss rate that is changing) InitCoef = RateCoefficients[3] //log initial rate coefficient //Step 1. Calculate the initial difference between model and measurement, do nothing if already within fit accuracy InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, bkgO3) dN2O5 = model_N2O5_ppbv[49]-N2O5_t fN2O5 = dN2O5/N2O5_t If (abs(fN2O5) < gN2O5) Stop =1 Iterations[2] = 1 //log the number of iterations InitialConcentrations = InitConcWave //set initial concentrations back to original values Return Stop //"Return" stops the rest of the code from running if this condition is met Endif If (dN2O5 < 0) //determine whether initial kN2O5 needs to increase or decrease slast = -1 Else slast = 1 Endif //Step 2. Determine Initial Order of Magnitude of kN2O5 (e.g. kN2O5 = x*10^order) and test if the solution has the same order of magnitude order = -1*(round(log(RateCoefficients[3]))) // round because IGOR will sometimes record 1e-4 as 9.9999e-5 and order would be reported incorrectly interger = round(RateCoefficients[3]*10^order) //calculate the closest interger number (i.e. 1 for 9.999e-5*1e4) remainder = mod(RateCoefficients[3]*10^order, 1) //find the remainder of ratecoefficient * order (i.e. (4e-4*1e4)/1 = 4/1 = 4) so rounding can be done correctly If(interger == 1 || interger == 10) //this is the case when IGOR has turned 1e-X into 9.999e-X If(remainder >= 0.99) //if remainer is also .99 or greater, it is another indicator that IGOR incorrectly reported the original kN2O5 value order = -1*(ceil(log(RateCoefficients[3]))) //return the upper value to ensure 9.999e-4 has order 3 and not order 4 Else order = -1*(round(log(RateCoefficients[3]))) //if interger were 1, the order would have been initially calculated properly Endif Else //if the starting value is not 1e-X or 9.99e-X.... order = -1*(floor(log(RateCoefficients[3]))) //if order was reported by IGOR properly, round to the lowest value (i.e. to ensure 4e-4 has order 4 and not 3 as would be given by log(4e-4)) If(slast == -1) RateCoefficients[3] = 1*10^-order //if kN2O5 needs to go down from the original value, then go down (i.e. 4e-4 --> 1e-4) Else RateCoefficients[3] = 1*10^-(order-1) //if kN2O5 needs to go up...(i.e. 4e-4 --> 1e-3) Endif KK[3] = RateCoefficients[3] InitialConcentrations= InitConCWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run the model with adjusted kN2O5 dN2O5 = model_N2O5_ppbv[49]-N2O5_t //re-calculate the model/obs difference n+=1 If (dN2O5 < 0) //determine whether kN2O5 needs to increase or decrease snew = -1 Else snew = 1 Endif If(snew != slast) //if snew and slast are not the same, the original kN2O5 and the adjusted value are the bounds CalcMag = 0 //order of magnitude loop is not required Counter = 1 //start step 4 with initial kN2O5 adjusted up or down by appropriate step size Else Counter = 1 CalcMag = 1 //if snew and slast are the same, calculation of the order of magnitude is required Endif Endif olast = order //log current order of magnitude //ERROR Checking If(order >= 10) //error checking - make sure kN2O5 does not get too small Stop = 1 Iterations[2] = 6 RateCoefficients[3] = 1e-7 //SET RATE COEFFICIENT TO LOWER LIMIT so the NO2, O3, N2O5 loop can try again KK[3] = RateCoefficients[3] InitialConcentrations = InitConcWave Return Stop Endif If(order <= 0) //error checking - make sure kN2O5 does not get too large Stop = 1 Iterations[2] = 6 //SET RATE COEFFICIENT TO LOWER LIMIT so the NO2, O3, N2O5 loop can try again RateCoefficients[3] = 1 kk[3] = RateCoefficients[3] InitialConcentrations = InitConcWave Return Stop //Stop iteration Endif //Step 3. Find Order of Magnitude Upper and Lower Bounds - Start with the initial order of magnitude and change the magnitude up or down until the upper and lower magnitude bounds have been found If(CalcMag == 1) Do If (dN2O5 < 0) //check whether kN2O5 needs to increase of decrease (will be the same as slast during first iteration) snew = -1 Else snew = 1 Endif If (snew == slast) //if the model is still under/over-predicitng relative to the previous iteration, adjust the order of magnitude up/down by 1 (slast is initial direction) order -= snew //change order of magnitude by 1 (3--> 4 if kN2O5 needs to go down (snew = -1), 3-->2 if kN2O5 needs to go up (snew = 1)) RateCoefficients[3] = 1*10^-order //adjust kN2O5 with new order KK[3] = RateCoefficients[3] InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, bkgO3) //re-run the model with adjusted kN2O5 dN2O5 = model_N2O5_ppbv[49]-N2O5_t //re-calculate the model/obs difference olast = order //log the current order of magnitude counter +=1 //logs whether order of magnitude has changed n +=1 //counts to number of iterations Else Stop =1 //if the model is now over/underprediciting relative to under/overpredicting of the last iteration, the order of magnitude bounds have been found and recorded as order and olast order -= snew // record the new order as a bound, but don't change the rate coefficient Endif If(order >= 10) //error checking to make sure kN2O5 does not get too large or small - ends function if error occurs Stop = 1 Iterations[2] = 6 RateCoefficients[3] = 1e-7 kk[3] = RateCoefficients[3] InitialConcentrations = InitConcWave Return Stop Endif If(order <= 0) Stop = 1 Iterations[2] = 6 RateCoefficients[3] = 1 kk[3] = RateCoefficients[3] InitialConcentrations = InitConcWave Return Stop Endif While(Stop ==0) Endif //Step 4. Set Size Step Based on Magnitude Bounds and Re-Run Model - step equal to the lower order of magnitude bound (i.e. if range in 1e-4 - 1e-3, then Step Size = 1e-4) If(order < olast) //if order is upper bound and olast is lower bound, step size is equal to olast magnitude StepSize = 1*10^-(olast) RateCoefficients[3] = 2*StepSize Else //if olast is upper bound and order is lower bound, step size is equal to order magnitude StepSize = 1*10^-order RateCoefficients[3] = 2*StepSize Endif If(Counter == 1) //do this procedure if the order is the same as the initial order RateCoefficients[3] = InitCoef + slast*StepSize //start with the initial kN2O5 value, adjusted by the step size (since the initial value was determined in step 1 to not be within the fit error) KK[3] = RateCoefficients[3] InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run the model with new kN2O5 dN2O5 = model_N2O5_ppbv[49]-N2O5_t If (dN2O5 < 0) //determine whether kN2O5 needs to increase or decrease snew= -1 Else snew = 1 Endif If (snew != slast) //if the current direction is different than the previous direction, new upper and lower bounds have been found and the Step Size should be cut in half (slast is initial direction) StepSize /= 2 Endif slast = snew //reset slast for Step 5 Else //do this procedure if the order has changed from it's initial vlaue KK[3] = RateCoefficients[3] InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run model with new kN2O5 determined above dN2O5 = model_N2O5_ppbv[49]-N2O5_t If (dN2O5 < 0) //determine whether kN2O5 needs to increase or decrease slast= -1 //reset slast to new value Else slast = 1 Endif Endif n += 1 //count number of iterations (i.e. model runs) // Step 5. Iteratively solve for kN2O5 - within determined order of magnitude Do fN2O5 = dN2O5/N2O5_t If (abs(fN2O5) > gN2O5) //if model-obs. difference within accuracy, stop model, if not adjust kN2O5 rate coefficient Stop =0 RateCoefficients[3] += slast*StepSize //adjust the rate coefficient up or down according to current step size and direction (i.e. slast) KK[3] = RateCoefficients[3] InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run model with new kN2O5 dN2O5 = model_N2O5_ppbv[49]-N2O5_t If (dN2O5 < 0) //determine whether kN2O5 needs to increase or decrease snew = -1 Else snew = 1 Endif If ( snew != slast) //if model has overshot, then decrease the step size StepSize /= 2 slast = snew //change direction of step Endif n += 1 //counter iterations Else Stop = 1 Endif While (Stop == 0 && n < 20) Iterations[2] = n //log the number of iterations InitialConcentrations = InitConcWave //reset initial concentrations for the next procedure End Function //************************************************************************************************ Function NUM_Fit_ClNO2_MultiStep(dtpresunset, dt, ClNO2, ktot, BkgO3) //Macro iteratively solves ODEs while fitting measured ClNO2, using total N2O5 loss rate and adjusting ClNO2 production rate Variable dtpresunset //time between start of chemistry before sunset and sunset Variable dt //nocturnal simulation duration (hr) Variable ClNO2 //measured ClNO2 (ppbv) Variable ktot //total N2O5 loss rate constant from NO2, O3, N2O5 function Variable BkgO3 //ppbv bkg O3 Variable SumIter, Error Wave ErrorWave, Iterations, KK, RateCoefficients NUM_FitClNO2_MultistepFn(dtpresunset, dt, ClNO2, ktot, BkgO3) //Fit the modeled ClNO2 by iterating the ClNO2 production rate constant SumIter = Iterations[3] If(SumIter >= 30) //model probably wont converge after 30 tries Error = 1 Endif Return Error End Function //************************************************************************************************ Function NUM_FitClNO2_MultiStepFn(dtpresunset, dt, ClNO2_t, ktot, BkgO3) //Iterates to fit final ClNO2 concentration by adjusting the production rate of ClNO2 and the loss rate constant of N2O5 keeping the total N2O5 loss rate constant Variable dtpresunset //time between start of chemistry and sunset Variable dt //nocturnal simualtion duration Variable ClNO2_t //measured ClNO2 (ppbv) Variable ktot //total loss rate of kN2O5 (derived from FitN2O5 prodcedure) Variable BkgO3 //ppbv BkgO3 Variable dClNO2, fClNO2 //absolute and fractional HNO3 differences between calculated and measured downwind values Variable Stop = 0 // Stop (=1), keep going ( = 0) Variable n = 1 //Iteration counter variable gClNO2 = 0.01 //Minimum Fit threshold variable snew, slast //Tracker variables for direction. order down = -1, order up = +1 Variable StepSize //Size of iteration step Variable order, olast //tracker variables for order of magnitude Variable remainder, interger //variables to ensure correct order calculation Variable counter //counter for iterations in calculating the order of magnitude Variable InitCoef //Log of Initial kClNO2 Wave model_clno2_ppbv, ratecoefficients, iterations, KK, InitConcWave, InitialConcentrations InitConcWave = InitialConcentrations //log initial concentrations //Step 1. Set initial value to 1e-order, where order is the same as that of the derived ktot order = -1*(round(log(ktot))) interger = round(ktot*1*10^order) remainder = mod(ktot*1*10^order, 1) //find the remainder so rounding can be done correctly If(Interger == 1 || interger == 10) If(remainder >= 0.99) //if coefficient is order = -1*(ceil(log(ktot))) Endif Else order = -1*(floor(log(ktot))) //check the order of magnitude of the incoming kN2O5 Endif olast = order //Step 2. Determine if this initial value is within fit accuracy RateCoefficients[4] = 1*10^-order InitCoef = RateCoefficients[4] //log initial kClNO2 RateCoefficients[3] = ktot - ratecoefficients[4] //must keep the total N2O5 loss rate constant while changing the fit to ClNO2 If(RateCoefficients[3] <0) //if the prodction rate of ClNO2 becomes larger than the total N2O5 loss rate, set HNO3 production to zero and continue to calc. ClNO2 (will allow model to get yield over 100%) RateCoefficients[3] =0 Endif kk = ratecoefficients InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) dClNO2 = model_ClNO2_ppbv[49]-ClNO2_t fClNO2 = dClNO2/ClNO2_t //calc. fractional ClNO2 difference between calculated and measured downwind values If (abs(fClNO2) < gClNO2) //if this initial value is within fit accuracy, stop Stop = 1 InitialConcentrations = InitConcWave Iterations[3] = 1 Return Stop Endif If (dClNO2 < 0) //determine if kClNO2 needs to increase or decrease slast = -1 Else slast = 1 Endif //Step 3. Find order of magnitude of solution counter = 0 Do If (dClNO2 < 0) //check to see whether kClNO2 needs to increase or decrease snew = -1 //if kClNO2 needs to go up (order down), then it is most liekly not an order of magnitude larger than the magnitude of ktot and this loop is stopped Stop = 1 Else snew = 1 Stop =0 Endif If (snew == slast &&Stop ==0) //if the model is still under/over-predicitng relative to the previous step, change order of magnitude accordingly. If not, the solution is within the current order of magnitude order += snew //change order of magnitude ( up if over predicting, down if under predicting) RateCoefficients[4] = 1*10^-order RateCoefficients[3] = ktot - ratecoefficients[4] //must keep the total N2O5 loss rate constant while changing the fit to ClNO2 If(RateCoefficients[3] <0) //if the prodction rate of ClNO2 becomes larger than the total N2O5 loss rate, set HNO3 production to zero and continue to calc. ClNO2 (to get yeilf over 100%) RateCoefficients[3] =0 Endif kk = ratecoefficients InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) dClNO2 = model_ClNO2_ppbv[49]-ClNO2_t olast = order //if magnitude order is changing, record the previous value n+=1 //iteration counter counter +=1 //logs whether order of magnitude has changed Else Stop =1 //if slast and snew are different, the bounds are order and olast order += snew Endif If(order >= 9) //Error checking to make sure kClNO2 is within a reasonable range Stop = 1 Iterations[3] = 30 RateCoefficients[4] = 0 kk[4] = RateCoefficients[4] Return Stop //stop procedure if an error has occured Endif If(order < 0) //Error checking to make sure kClNO2 is within a reasonable range Stop = 1 Iterations[3] = 30 RateCoefficients[4] = 0 kk[4] = RateCoefficients[4] Return Stop //stop procedure if an error has occured Endif While(Stop ==0) //Step 4. Set Step Size If(order < olast) //if order is upper bound and olast is lower bound, step size is equal to olast magnitude StepSize = 1*10^-(olast) RateCoefficients[4] = 2*StepSize Else //if olast is upper bound and order is lower bound, step size is equal to order magnitude StepSize = 1*10^-order RateCoefficients[4] = 2*StepSize Endif If(Counter == 0) //if order hasn't changed, start with initial kClNO2 value RateCoefficients[4] = ktot/2 //slast is initial direction RateCoefficients[3] = ktot - ratecoefficients[4] //must keep the total N2O5 loss rate constant while changing the fit to ClNO2 If(RateCoefficients[3] <0) //if the prodction rate of ClNO2 becomes larger than the total N2O5 loss rate, set HNO3 production to zero and continue to calc. ClNO2 (to get yeilf over 100%) RateCoefficients[3] =0 Endif kk = ratecoefficients InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run model with adjust kClNO2 dClNO2 = model_ClNO2_ppbv[49]-ClNO2_t If (dClNO2 < 0) //determine if kClNO2 needs to increase or decrease snew = -1 Else snew = 1 Endif If ( snew != slast) //if new direction of first step is different than intitial direction, cut step size in half StepSize /= 2 slast = snew //reverse direction Endif Else RateCoefficients[3] = ktot - ratecoefficients[4] //If order of magnitude has changed... If(RateCoefficients[3] <0) RateCoefficients[3] =0 Endif kk = ratecoefficients InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run model dClNO2 = model_ClNO2_ppbv[49]-ClNO2_t If (dClNO2 < 0) //check whether kClNo2 needs to increase or decrease slast = -1 Else slast = 1 Endif Endif n +=1 //Step 5. Iteratively find kClNO2 (within order of magnitude) Do fClNO2 = dClNO2/ClNO2_t //calc. fractional ClNO2 difference between calculated and measured downwind values If (abs(fClNO2) > gClNO2) Stop = 0 RateCoefficients[4] -=slast*StepSize //adjust kClNo2 RateCoefficients[3] = ktot - ratecoefficients[4] //must keep the total N2O5 loss rate constant while changing the fit to ClNO2 If(RateCoefficients[3] <0) //if the prodction rate of ClNO2 becomes larger than the total N2O5 loss rate, set HNO3 production to zero and continue to calc. ClNO2 (to get yeilf over 100%) RateCoefficients[3] =0 Endif kk = ratecoefficients InitialConcentrations = InitConcWave NUM_Model_MultiStep(dtpresunset, dt, BkgO3) //re-run model dClNO2 = model_ClNO2_ppbv[49]-ClNO2_t If (dClNO2 < 0) //check if kClNo2 needs to increase or decrease snew = -1 Else snew = 1 Endif If ( snew != slast) //If model has overshot... StepSize /= 2 //cut step size in half slast = snew //reverse direction Endif n +=1 Else Stop = 1 Endif While (Stop == 0 && n < 30) Iterations[3] = n //count the number of iterations InitialConcentrations = InitConcWave Return Stop End Function //************************************************************************************************ Function NUM_Model_MultiStep(dtpresunset, dt, bkgO3) //Do the model calculation for 2 step integration. Step 1 = before sunset(3 steps), Step 2 = aftersunset - only adjusts the RateCoefficient and KK waves, not any initial concentrations Variable dtpresunset //time between start of chemistry and susnet Variable dt //nocturnal simulation duration variable bkgO3 //ppbv of background O3 Wave model_n2O5_ppbv, ratecoefficients, iterations, KK, InitConcWave, CalTime, InitialConcentrations string cmd, wname, model_wname Variable i Variable istep //current step before sunset Variable totalsteps //total number of time steps before sunset Variable timestep //the time for the current ODE iteration Wave/T Names //names of mechansim species, from Model SetUp NVAR vnDensity wave model_no2_ppbv, model_o3_ppbv, model_n2o5_ppbv, model_clno2_ppbv, model_no_ppbv //determine the number of pre-sunset steps to take, depending on the length of dtpresunset (based on an ideal 3 steps for a dtpresunset of 1.3 hours) If(dtpresunset > 0) //reduce unnecessary number of steps in some instantces If(dtpresunset > 1) totalsteps = 3 Endif If(dtpresunset < 1 && dtpresunset > 0.5) totalsteps = 2 Endif If(dtpresunset < 0.5) totalsteps = 1 Endif timestep = dtpresunset / totalsteps //resolution of each of the 3 time intervals before sunset CalTime = timestep*3600*(x/49) Variable jNO3_value Variable jN2O5_value Variable jClNO2_value Variable jNO2_value Variable jO3_Value //Do this loop for the calculated number of pre-sunset steps For (istep = 0; istep < totalsteps; istep += 1) jNO3_value= ((-0.29+0.3*exp(-0.26*-1*(dtpresunset-(timestep*istep))))+(-0.29+0.3*exp(-0.26*-1*(dtpresunset-(timestep*(istep+1))))))/2 //avg between two timesteps, fit to jNO3 data - derived from CalNEX jNO3 as a function of jNO2 and SZA (average of all day --> night flights) jN2O5_value =((-9e-6 +1e-5*exp(-0.47*-1*(dtpresunset-(timestep*istep))))+(-9e-6 +1e-5*exp(-0.47*-1*(dtpresunset-(timestep*(istep+1))))))/2 //avg between two timesteps,, fit to jN2O5 data - measured jClNO2_value = ((-1.4e-4 +1.5e-4*exp(-0.42*-1*(dtpresunset-(timestep*istep))))+(-1.4e-4 +1.5e-4*exp(-0.42*-1*(dtpresunset-(timestep*(istep+1))))))/2 //avg between two timesteps,, fit to jClNO2 data - measured jNO2_value = ((-7.4e-3 +7.5e-3*exp(-0.25*-1*(dtpresunset-(timestep*istep))))+(-7.4e-3 +7.5e-3*exp(-0.25*-1*(dtpresunset-(timestep*(istep+1))))))/2 //avg between two timesteps,, fit to jNO2 data - measured jO3_Value = ((-5.2e-7 +4.1e-7*exp(-1.08*-1*(dtpresunset-(timestep*istep))))+(-5.2e-7 +4.1e-7*exp(-1.08*-1*(dtpresunset-(timestep*(istep+1))))))/2 //avg between two timesteps,, fit to jO3 data - measured RateCoefficients[10] = jNO3_value //assign values to RateCoefficients wave RateCoefficients[11] = jN2O5_value RateCoefficients[12] = jClNO2_value RateCoefficients[9] = jNO2_value RateCoefficients[13] = jO3_Value If(RateCoefficients[9] < 0) RateCoefficients[9] = 0 Elseif(RateCoefficients[10] < 0) RateCoefficients[10] = 0 Elseif(RateCoefficients[11] < 0) RateCoefficients[11] = 0 Elseif(RateCoefficients[12] < 0 ) RateCoefficients[12] = 0 Elseif(RateCoefficients[13] < 0) RateCoefficients[13] = 0 Endif KK = RateCoefficients cmd = "NUM_Run_Model2()" Execute cmd //integrate ODEs with current initial concentrations for the duration of the current istep For(i = 0; i < (numpnts(Names)-2); i += 1) //adjust initial concentrations for the next step of model run wname = Names[i] model_wname = "Model_" + wname +"_ppbv" Wave modelwaveref = $model_wname InitialConcentrations[i] = (modelwaveref[49]*1e-9*vndensity) Endfor InitialConcentrations[9] = bkgo3*1e-9*vndensity Endfor Endif //post sunset RateCoefficients[10] = 0 //no more photolysis RateCoefficients[11] = 0 RateCoefficients[12] = 0 RateCoefficients[9] = 0 RateCoefficients[13] = 0 KK = RateCoefficients CalTime = dt*3600*(x/49) cmd = "NUM_Run_Model2()" //integrate ODEs for the nocturnal simulation suration Execute cmd End Function //// ************************************************************************************** //// ************************************************************************************** // The following Functions/Macros: ***************************************************** // // 1. NUM_IUPAC_NO3NO2() // 2. NUM_IUPAC_UniN2O5() // 3. NUM_ReplaceNaNZero() // 4. NUM_Calc_SZAFn() // 5. NUM_Time2Sunset() // 6. NUM_Time2SunsetFn() // 7. NUM_Time2Sunrise() // 8. NUM_Time2SunriseFn() //// 9. NUM_InterpNaNFn() // 10. NUM_Filter2WavesFn() // 11. NUM_FilterWaveFn() // 12. NUM_MakeAvgTimeWaveFn() // // are used to support the Model Setup and Model Model Macros // // Author:Steven Brown // Modified by Erin McDuffie (2018) //// ****************************************************************************************** //********************************************************************************************* Function NUM_IUPAC_NO3NO2(P, T) //calculates the rate equation for NO2+NO3 --> N2O5 - IUPAC (2012) Variable P // mbar Variable T // K Variable M = P*1e-4/(1.381e-23*T) //complex rate reaction from MCM 3.2 Variable K0 = 3.6e-30*M*((T/300)^-4.1) Variable K1 = 1.9e-12*((T/300)^0.2) Variable KR = K0/K1 Variable FC = 0.35 Variable NC = 0.75-1.27*log(FC) Variable F = 10^(log(FC)/(1+(log(KR)/NC)^2)) Variable KMT = (K0*K1)*F/(K0+K1) //print KMT Return KMT //3.6e-30*((T/300)^-4.1)*M End Function //// ************************************************************************************** Function NUM_IUPAC_UniN2O5(T, P) //calculates the rate constant for unimolecular dissociation of N2O5 --> NO3 + NO2 - IUPAC (2012) Variable T //K Variable P //mbar Variable M = P*1e-4/(1.381e-23*T) Variable K0 = 1.3e-3*M*(T/300)^-3.5*exp(-11000/T) Variable K1 = 9.7e14*(T/300)^0.1*exp(-11080/T) Variable KR = K0/K1 Variable FC = 0.35 Variable NC = 0.75-1.27*log(FC) Variable F = 10^(log(FC)/(1+(log(KR)/NC)^2)) Variable KMT = (K0*K1)*F/(K0+K1) Return KMT// //9.7e14*((T/300)^0.1)*exp(-11080/T) End Function //// ************************************************************************************** Function NUM_ReplaceNaNZeroFn(w1, w2) //Function takes a wave (w1) with real numbers and NaN values and makes a new wave (2) //with only real numbers. Previious NaN are set to zero Wave w1 //Original Wave String w2 //New Wave without NaNs Duplicate/O w1, $w2 Wave wave2 = $w2 variable ntot = numpnts(w1) variable n = 0 Silent 1 PauseUpdate Do If(numtype(w1[n]) == 2) wave2[n] = 0 Endif n += 1 While (n < ntot) End Function //// ************************************************************************************** Function NUM_Calc_SZAFn(lat, lon, Year, twave, offset) //Function returns the SZA in degrees Variable lat, lon, Year, twave, offset //latitude, longitude in degrees, calendar year, local time in seconds since 1/1/04, offset in hours to UTC (offset = 0 if timewave is in UTC) Variable Day, Hour Variable gam, b0, b1, b2, b3, b4, b5, b6, EQT,c0, c1, c2, c3, c4, delta, w Day = floor((twave +offset*3600 - date2secs(Year,1,1))/86400) //integer day number Hour = ((twave +offset*3600 - date2secs(Year, 1, 1) - Day*86400)/86400)*24 //decimal hour number //Calculate solar declination (+/- 23 degrees) in radians according to empirical formula (source: http://www.met.wau.nl/Courses/Micrometcourse/Modules/Shortwave/ModulePage3.html) If (mod(Year, 4) == 0) //Test for leap year gam = 2*pi*(Day -1)/366 Else gam = 2*pi*(Day -1)/365 Endif b0 = 0.006918 b1 = -0.399912 b2 = 0.070257 b3 = -0.008758 b4 = 0.000907 b5 = -0.002697 b6 = 0.00148 delta = b0 + b1*cos(gam) + b2*sin(gam) + b3*cos(2*gam) + b4*sin(2*gam) + b5*cos(3*gam) + b6*sin(3*gam) //radians //Calculate hour angle (radians) c0 = 7.5e-5 c1 = 1.868e-3 c2 = -0.032077 c3 = -0.014615 c4 = -0.040849 EQT = c0 + c1*cos(gam) + c2*sin(gam) + c3*cos(2*gam) + c4*sin(2*gam) w = pi*(Hour/12 - 1 + lon/180) + EQT //radians return (180/pi)*acos(sin(delta)*sin(lat*pi/180) + cos(delta)*cos(lat*pi/180)*cos(w)) //Solar Zenith Angle in degrees (SZA = sin(delta)*sin(local latitude) + cos(local latitude)*cos(delta)*cos(hour angle) - wikipedia) End Function //// ************************************************************************************** Macro NUM_Time2Sunset(Year, wtime, Lat, lon) //makes dtSunset wave and calculates the amount of time (hours) since sunset for each timepoint based on aircraft GPS location and time Variable Year = 2015 String wtime = "Start_UTC" String lat = "RAF_Latc" String lon = "RAF_Lonc" Prompt Year,"Enter Year" Prompt wtime, "Time Wave" Prompt lat, "Latitude Wave" Prompt lon, "Longitude Wave" NUM_Time2SunsetFn(Year, $wtime, $lat, $lon) //makes and calculates dtSunset wave End Macro //// ************************************************************************************** Function NUM_Time2SunsetFn(Year, wtime, wlat, wlon) //Makes dtSunset wave and calculates the time (hours) since sunet for each wave point based on aircraft GPS location Variable Year Wave wtime Wave wlat Wave wlon Variable SZA, i, j, jmax Duplicate/O/D wlat, dtSunset //time since sunset wave dtSunset = 0 jmax = numpnts(wlat) j = 0 i = 0 Do //calculate the time since sunset (based on SZA) at every wave point of the given flight SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) If (SZA < 90) Do //move forwards in time (+i) until calculated SZA > 90 i += 1 SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) While (SZA < 90) Else //move backwards in time (-i) until calculated SZA < 90 Do i -= 1 SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) While (SZA > 90) Endif dtSunSet[j] = -i/3600 //time since sunset in hours (negative if the sun is still up, positive if the sun has already set) j += 1 //loop through entire flight While (j < jmax) End Function //// ************************************************************************************** Macro NUM_Time2Sunrise(Year, wtime, Lat, lon) //makes dtSunset wave and calculates the amount of time (hours) since sunset for each timepoint based on aircraft GPS location and time Variable Year = 2015 String wtime = "Start_UTC" String lat = "RAF_Latc" String lon = "RAF_Lonc" Prompt Year, "Enter Year" Prompt wtime, "Time Wave" Prompt lat, "Latitude Wave" Prompt lon, "Longitude Wave" NUM_Time2SunriseFn(Year, $wtime, $lat, $lon) End Macro //// ************************************************************************************** Function NUM_Time2SunriseFn(Year, wtime, wlat, wlon) //Makes dtSunset wave and calculates the time (hours) since sunet for each wave point based on aircraft GPS location Variable Year Wave wtime Wave wlat Wave wlon Variable SZA, i, j, jmax Duplicate/O/D wlat, dtSunrise //time since sunset wave dtSunrise = 0 jmax = numpnts(wlat) j = 0 i = 0 Do //calculate the time since sunset (based on SZA) at every wave point of the given flight SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) If (SZA > 90) Do //move forwards in time (+i) until calculated SZA < 90 i += 1 SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) While (SZA > 90) Else Do //move backwards in time (-i) until calculated SZA > 90 i -= 1 SZA = NUM_Calc_SZAFn(wlat[j], wlon[j], Year, wtime[j]+i,0) While (SZA < 90) Endif dtSunRise[j] = i/3600 //time until sunrise in hours (negative if the sun has risen, positive if the sun is still down) j += 1 //loop through entire flight While (j < jmax) End Function //// ************************************************************************************** Function NUM_InterpNaNFn(w1, w2) //Function takes a wave (w1) with real numbers and NaN or Inf values and makes a new wave (2) //with only real numbers. The real numbers are an interpolation from the last real number //prior to the NaN or Inf sequence in w1 to the next real number afterwards. The interpolation //is linear and assumes equally spaced points. Wave w1 //wave to be interpreted String w2 //new wave with interpolated values Duplicate/o w1 $w2 Wave wave2 = $w2 Variable ntot = numpnts(w1) Variable n = 0 Variable n1, n2, y1, y2 Silent 1 PauseUpdate Do If(numtype(w1[n]) == 0) wave2[n] = w1[n] n += 1 Else n1 = n-1 If(n1 <0) n1 = 0 Endif y1 = w1[n1] Do n += 1 While((numtype(w1[n]) != 0) && (n < ntot-1)) If (n == ntot) n2 = n-1 wave2[n1,n2] = y1 Else n2 = n y2 = w1[n] wave2[n1,n2] = w1[n1] + ((y2-y1)/(n2-n1))*(x-n1) Endif Endif While (n < ntot-1) End Function //// ************************************************************************************** Function NUM_Filter2WavesFn(name1, name2, low,high) //filters name1 for low and high values in name2 Wave name1 //wave to be filtered Wave name2 //wave whose values will be the filter Variable low //lower threshold (<) Variable high //uperthreshol (>) Variable imax = numpnts(name1) Variable i = 0 Silent 1 PauseUpdate Do If (name2[i] > high) name1[i] = NaN Endif If (name2[i] < low) name1[i] = NaN Endif i += 1 While (i < imax) End Function //// ************************************************************************************** Function NUM_FilterWaveFn(name,low,high) //filters a wave for values less than and greater than the input limits Wave name //wave to be filtered Variable low, high //low and high filter limits (< and >) Variable imax = numpnts(name) Variable i = 0 Silent 1 PauseUpdate Do If (name[i] > high) name[i] = NaN Endif If (name[i] < low) name[i] = NaN Endif i += 1 While (i < imax) End Function //// ************************************************************************************** Function NUM_MakeAvgTimeWaveFn(timewave, interval) //Makes waves named "Start_interval+"s" and "Stop_interval+"s" from timewave Wave timewave //timewave to average Variable Interval //average interval Variable n = numpnts(timewave)/Interval String IntStr = num2str(interval) Make/n=(n)/D/O Temp String Start = "Start_" + IntStr + "s" String Stop = "Stop_" + IntStr + "s" String TimeStr = "Time_" + IntStr + "s" Duplicate/o Temp $Start, $Stop, $TimeStr Wave Temp1 = $Start Wave Temp2 = $Stop Wave Temp3 = $TimeStr Temp1 = timewave[0] + x*Interval Temp2 = (Temp1) + (Interval - 1) Temp3 = (Temp1) + 0.5*Interval SetScale d 0,0,"dat", Temp1, Temp2, Temp3 Edit $Start, $Stop, $TimeStr ModifyTable format($Start)=8,format($Stop)=8,format($TimeStr)=8 Killwaves/z temp, Temp1, Temp2, Temp3 End Function //// ************************************************************************************** //// ************************************************************************************** // The following Functions/Macros: ***************************************************** // //// 1. NUM_DoAveragUsingStartStopFn() // 2. NUM_FindFirst() // 3. NUM_FindFirstGE() // 4. NUM_FindRightMostLE() // 5. NUM_ExpandItFillFn() //// // are used to support the Model Setup and Main Model macros // // Author:Donna Sueper // Modified by Ken Aikin (2015) and Erin McDuffie (2018) //// ************************************************************************************** //// ************************************************************************************** Function NUM_DoAveragUsingStartStopFn(timeWave, wave2avg, startWave, stopWave) // Calculates the average, standard deviation and number of non-nan points between time intervals start and stop times waves Wave timeWave //Data timewave Wave wave2avg //Data to average Wave startWave //Wave with average start times Wave stopWave //Wave with average stop times Make/o/n=(numpnts(startWave)) waveAvg, wave_std, wave_numpnts Variable start, stop, n=numpnts(timeWave) // start and stop range over the values 0 through numpnts(timeWave) Variable i=0, numpts_start=numpnts(startWave) // i ranges over the values 0 through numpnts(startWave) WaveStats/q timeWave If (V_numNans !=0) Abort "The timewave "+NameofWave (timeWave) + " has nans - Aborting from DoAveragUsingStartStop" Endif WaveStats/q startWave If (V_numNans !=0) Abort "The timewave "+NameofWave (startWave) + " has nans - Aborting from DoAveragUsingStartStop" Endif WaveStats/q stopWave If (V_numNans !=0) Abort "The timewave "+NameofWave (stopWave) + " has nans - Aborting from DoAveragUsingStartStop" Endif waveAvg=nan wave_std=nan wave_numpnts=nan start = 0 stop = 0 For(i=0;i=startWave[i])&&(timeWave[stop]<=stopWave[i]) ) If ( (stop - start)>=1) WaveStats/q/r=[start,stop] wave2avg waveAvg[i]=V_avg wave_std[i]=V_sdev wave_numpnts[i]=V_npnts Else If ( (stop - start)==0) waveAvg[i]=wave2avg[start] wave_std[i]=0 wave_numpnts[i]=1 Endif Endif Endif Endfor String AvgStr = NameofWave(wave2avg)+"_avg" Duplicate/o waveAvg $AvgStr String stdStr = NameofWave(wave_std)+"_Std" Duplicate/o wave_std $stdstr String npntsstr = NameofWave(wave_numpnts)+"_npnts" Duplicate/o wave_numpnts $npntsstr killwaves/Z waveavg, wave_std, wave_numpnts End //// ************************************************************************************** // Returns the first non-nan value after and including beginPt in a wave. // If all values after beginPt are nans, Findfirst will return number of points. Function NUM_FindFirst(w, beginPt) Wave w Variable beginPt variable i=beginPt, n=numpnts(w) If (i >= n) Return n Endif If (i < 0) Return -1 Endif If (numtype(w[i])== 0 ) // is w[beginPt] a number? Return i Endif For (i=beginPt; i=value. // It returns n if all the points in w are less than value // It returns -1 if beginPt < 0 // Typically used with waves that increase with increasing point number, like timewaves. Function NUM_FindFirstGE(w, value, beginPt) Wave w Variable value, beginPt variable i=beginPt, n = numpnts(w) If (i >= n) Return n Endif If (i < 0) Return -1 Endif If (w[i] >= value) Return i Endif For(i=beginPt; i=value) Break Endif Endfor Return i End //// ************************************************************************************** // Returns the right-most point less than or equal to value // in the wave w. It begin searching at beginPt and continues in increasing point number. // It returns beginPt if w[beginPt] > value. // This function assumes that w is sorted in increasing order. Function NUM_FindRightMostLE(w, value, beginPt) Wave w Variable value, beginPt variable i=beginPt, n= numpnts(w) If (i >= n) Return n Endif If (i < 0) Return -1 Endif If (w[i] > value) Return beginPt //=i Endif If (w[n-1] <= value) // the last point in w is less than or equal to value Return (n-1) Endif Do If (w[i+1] <=value) i+=1 Endif While( (i<(n-1)) && (w[i] <=value) && ( w[i+1] <= value ) ) Return i End //// ************************************************************************************** // Expands data which is on StartTime/StopTime, onto another wave such as AOCTimewave, and fills in the values between // the start/stop time with the value that was in the data wave. // Useful for example to expand CH2O data from Start/Stop time onto AOCtimewave. Function NUM_ExpandItFillFn(start1, stop1, w1, t2, w2) Wave start1 //Wave with start time values Wave stop1 //Wave with stop time values Wave w1 //Wave you want expanded Wave t2 //Time wave you want to expand to Wave w2 //Expanded Data wave (must already exist) If (numpnts(Start1)!=numpnts(w1) ) Print "// ", Start1, w1, numpnts(Start1), numpnts(w1) Abort "The number of points must be the same; see the history window to help debug. Aborting from TimeExpandCreateWave" Endif Variable x1=0, x2=0 Variable n2=numpnts(t2), n1=numpnts(start1) w2=nan Do If ((t2[x2] >= start1[x1]) * (t2[x2] < stop1[x1])) //need to fill in between values Do w2[x2] = w1[x1] //set to same value between start and stoptime x2+=1 If (x2>n2) //exceeded length of wave t2, so there is no more time to expand onto or fill Break Endif While(( t2[x2] >= start1[x1] ) * ( t2[x2] < stop1[x1] )) w2[x2] = w1[x1] x1+=1 x2+=1 Else If ( t2[x2] < start1[x1] ) x2+=1 Else x1+=1 Endif Endif While ( ( x2 < n2 )&&(x1< n1 ) ) End //// ************************************************************************************** //// ************************************************************************************** // The following Functions/Macros: ***************************************************** // // 1. NUM_ChemMechanism_Setup() // 2. NUM_IncludeChemKineticFuncs() // 3. NUM_RemoveChemKineticFuncs() // 4. NUM_MakeChemKineticFuncs() // 5. NUM_Mechanism_Compile() // 5. NUM_Run_Model2() //// // are used to setup the reaction mechanism table, compile the mechanism, and compute the ODE integration // // Author:Unknown (from KinSim.ipf) // Modified by Erin McDuffie (2018) //// ************************************************************************************** //// ************************************************************************************** Macro NUM_ChemMechanism_Setup() // This macro will setup the input waves needed to generate the reaction mechanism // After running this macro, enter the reaction mechanism // Reactant and Product names need to be less than 11 characters in length // Do not use Reactant and Product as species names. // Enter the Rate coefficients for the reactions in the RateCoefficients wave // The time base is set for equally spaced points (100) over a period of 0.005 s. // This can be changed in the command window as desired later: i.e. caltime = 1e-3*exp(x/5) (unequally spaced points) Make/O/T/N=100 Reactant1, Reactant2, Product1, Product2, Product3, Names Reactant1 = "" Reactant2 = "" Product1 = "" Product2 = "" Product3 = "" Names = "" Make/O/N=50 RateCoefficients RateCoefficients=0 Make/O/N=50 CalTime=x*1e-5 Edit/W=(5.25,42.5,1122,236) reactant1,reactant2,product1,product2, product3, RateCoefficients DoWindow/K ReactionMechanism DoWindow/C/T ReactionMechanism,"Reaction Mechanism" ModifyTable size=11,width=90 ModifyTable width(ratecoefficients)=100 //Variable /G dilution = 0 NUM_IncludeChemKineticFuncs() // load subroutines for reaction mechanism from experiment folder End Function //*************************************************************************************************** Function NUM_IncludeChemKineticFuncs() // adds include statement to procedure window // note that ChemKineticsFuncs.ipf must be present in experiment folder // if not present, a blank one will be generated Pathinfo home // check if experiment was saved If (StrLen(S_Path)<1) abort "Save experiment first!" Endif String SelectStr=S_Path + "ChemKineticFuncs.ipf" String homepath = S_Path GetFileFolderInfo/Z/Q SelectStr // check if ChemKineticFuncs.ipf exists in same folder If (( V_Flag != 0) || (V_isFile==0)) // ChemKineticFuncs.ipf doesn't exist DoAlert 1, "ChemKineticFuncs.ipf does not exist, should I generate a blank one?" If (V_flag==1) // Yes clicked NUM_MakeChemKineticFuncs() Else abort "could not include ChemKineticFuncs.ipf" Endif Endif String cmd = "INSERTINCLUDE "+"\""+homepath+"ChemKineticFuncs\"" Execute /P cmd Execute /P "COMPILEPROCEDURES " End Function //*************************************************************************************************** Function NUM_RemoveChemKineticFuncs() // removes include statement from procedure window // also removes ChemKineticsFuncs.ipf from experiment Pathinfo home String cmd = "DELETEINCLUDE "+"\""+S_path+"ChemKineticFuncs\"" Execute /P cmd Execute /P "COMPILEPROCEDURES " End Function //*************************************************************************************************** Function NUM_MakeChemKineticFuncs() // makes new ChemKineticFuncs.ipf in experiment folder String noteBookName = "ChemKineticsFuncsText" NewNotebook /F=0 /N=$noteBookName as "ChemKineticFuncs.ipf" // make notebook NoteBook $notebookName text = "Function ChemKinetic(pw, tt, yw, dydt)"+"\r" // add hard-coded text into notebook NoteBook $notebookName text = "\tWave pw // rate coefficients"+"\r" NoteBook $notebookName text = "\tVariable tt // time value at which to calculate derivatives"+"\r" NoteBook $notebookName text = "\tWave yw // containing concentrations of species"+"\r" NoteBook $notebookName text = "\tWave dydt // wave to receive dA/dt, dB/dt etc. (output)"+"\r" NoteBook $notebookName text = ""+"\r" NoteBook $notebookName text = ""+"\tif (exists(\"dilution\"))\r" NoteBook $notebookName text = ""+"\t\tNVAR d = dilution\r" NoteBook $notebookName text = ""+"\tendif\r" NoteBook $notebookName text = "// begin dynamic code, don't modify this line!"+"\r" NoteBook $notebookName text = ""+"\r" NoteBook $notebookName text = "// end dynamic code, don't modify this line!"+"\r" NoteBook $notebookName text = ""+"\r" NoteBook $notebookName text = "End"+"\r" SaveNotebook/O /P=home $noteBookName as "ChemKineticFuncs.ipf" // Save the notebook DoWindow/K $noteBookName // Kill the notebook End Function //*************************************************************************************************** Macro NUM_Mechanism_Compile(fctname, PrintMech, Update_ChemKineticFuncs) //Compiles mechanism and makes ipf for use in "IntegrateODE function in Run_Model2() String fctname = "ChemKinetic" Variable printMech = 1 //print mechanism in command window (1 = yes) Variable Update_ChemKineticFuncs = 1 //update ChemKineticFuncs.ipf ( 1 = yes) (only required if the mechanism has changed, not the numeric values) Prompt fctname, "Select Kinetic Function", popup FunctionList("ChemKinetic*",";","KIND:2") Prompt Update_ChemKineticFuncs, "Update Chem Kinetic Functions Notebook? (1 = yes)" String/G CurrChemKineticFct = fctname Variable Rnumber, SpeciesNumber, SpeciesDiff Variable flag,n,m, nn, position, position2 Variable na,nb,nc,nd, ne Variable TotalLen, breakpos String temp, temp1, temp2, temp3, temp4, temp5 String TotalTemp, totaltemp1, totaltemp2, totaltemp3, totaltemp4, totaltemp5 //String dilutiontemp String sa,sb,sc,sd, se String s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10 String Plus, Equal String MyList, tempNote, tempNote2 String/G dydtstr = "" // string to be passed to ChemKineticsFuncs Silent 1 PauseUpdate Plus = "+" Equal = "=" s0 = "" s1 = " " s2 = S1 + " " s3 = s2 + " " s4 = s3 + " " s5 = s4 + " " s6 = s5 + " " s7 = s6 + " " s8 = s7 + " " s9 = s8 + " " s10 = s9 + " " n=0 If(PrintMech == 1) printf "\r\t\t\t\tReaction Mechanism \t\t\t\t\t\t\t\t\t\t Rate Coefficient\r" //print to command window Endif Do //Sets variables for each reactant and product for each reaction, prints to command window if specified temp = Reactant1[n] + Reactant2[n] + Product1[n] + Product2[n] + Product3[n] If (stringmatch(temp,"")) //Empty line in rxn. mechanism, have reached the end, jump to the end of the loop Rnumber = n // This is the number of reactions n += 500 Else //There is at least one reactant or product na = 10 - strlen(Reactant1[n]) temp = "s" + num2str(na) sa = Reactant1[n] + ($temp) If(strlen(Reactant2[n]) == 0) sa =Reactant1[n]+ ($temp) Endif If (strlen(Reactant1[n]) == 0) sa = " " Endif nb = 10 - strlen(Reactant2[n]) temp = "s" + num2str(nb) sb = plus + " " + Reactant2[n] +($temp)//+ equal If (strlen(Reactant2[n]) == 0) sb = " " //+ //equal Endif nc = 10 - strlen(Product1[n]) temp = "s" + num2str(nc) sc = equal + " " + Product1[n] +($temp) If (strlen(Product1[n]) == 0) sc = " " Endif nd = 10 - strlen(Product2[n]) temp = "s" + num2str(nd) sd = plus + " "+ Product2[n] + ($temp) If (strlen(Product2[n]) == 0) sd = " " Endif ne = 10 - strlen(Product3[n]) temp = "s" + num2str(ne) se = plus + Product3[n]+ ($temp) If (strlen(Product3[n]) == 0) se = " " Endif If(PrintMech == 1) print "\t\t\t\t" + sa + " "+ sb +" "+ sc + " " +sd+ " "+ se + "\t\t" + num2str(ratecoefficients[n]) //print each reaction as a line in the command window Endif Endif n += 1 While (n < 100) //Loop through the number of reaction If(PrintMech == 1) //Next, print the ODEs in compound format in command window if specified print "\r\r" Print " ODEs in Compound Format" Endif Make/O/N=(Rnumber) KK // Transfer Rate Coefficient Data to KK wave KK = RateCoefficients[p] n = 0 // Capture species names and write to Names[n] wave Make/O/N=600/T Names = "" Do nn = n*5 Names[nn] = Reactant1[n] Names[nn+1] = Reactant2[n] Names[nn+2] = Product1[n] Names[nn+3] = Product2[n] Names[nn+4] = Product3[n] n += 1 While (n < rnumber) n=numpnts(Names)-1 // last point in Names // Remove Blanks from Names List, changed to start at the end, since double empty names were not properly removed before Do If(strlen(Names[n]) == 0) DeletePoints n,1,Names Endif n -=1 While(n > -1) Silent 1 // Remove Duplicate Names from Names List Do flag=0 n=0 Do m=n+1 Do If(stringmatch(Names[m], Names[n])) DeletePoints m,1,Names flag=1 Else m += 1 Endif While(m < numpnts(Names)) n += 1 While(n< numpnts(Names) - 1) While (flag==1) n=0 // Determine number of Species in Reaction Names List Do If(stringmatch(Names[n],"")) SpeciesNumber = n n = 500 Endif n += 1 While (n < numpnts(Names)) SpeciesNumber = numpnts(Names) Make/O/N=(SpeciesNumber) InitialConcentrations n =0 // Make Output waves for Species with wave Names from Names[n] Do Make/O/D/N=(numpnts(CalTIME)) $Names[n] = 0 n +=1 While (n < SpeciesNumber) DoWindow ReactionMechanism //Append mechanism to Tbale "Reaction Mechansim" AppendToTable Names AppendToTable KK AppendToTable InitialConcentrations AppendToTable CalTime ModifyTable size(Names)=14, width(Names)=100 ModifyTable size(KK)=14,rgb(KK)=(65535,0,0), width(KK) = 100 ModifyTable size(initialconcentrations)=14, rgb(InitialConcentrations)=(65535,0,0), width(initialconcentrations)=150 ModifyTable size(CalTime)=14, rgb(CalTime)=(65535,0,0) If(printmech == 1) // Print ODEs for Reaction Mechanism in Chemical Format n = 0 Do m=0 Totaltemp = "" Totaltemp1 = "" Totaltemp2 = "" Totaltemp3 = "" Totaltemp4 = "" Totaltemp5 = "" Do // Find Reactant 1 in Names temp1 = "" If(stringmatch(Names[n], Reactant1[m])) temp1 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "] * ["+Reactant2[m] + "]" If(strlen(Reactant1[m]) == 0) temp1 = "-" + num2str(ratecoefficients[m]) + " * " + "["+Reactant2[m] + "]" Endif If(strlen(Reactant2[m]) == 0) temp1 = "-" + num2str(ratecoefficients[m]) + " * [" + Reactant1[m] + "]" Endif Endif Totaltemp1 += temp1 m += 1 While(m 0) If (TotalLen>390) // string too long for one line in function breakpos = StrSearch(Totaltemp, "+", 300) // find first "+" past 300 chars If (breakpos<0 || breakpos> 390) breakpos = StrSearch(Totaltemp, "-", 300)) // Check if there's maybe a "-" Endif If (breakpos>300) //Print "dydt["+ num2str(n) + "]=" + Totaltemp[0,breakpos-1] //Print "dydt["+ num2str(n) + "]+=" + Totaltemp[breakpos, StrLen(Totaltemp)-1] dydtstr += "dydt["+ num2str(n) + "]="+Totaltemp[0,breakpos-1] + "\r" dydtstr += "dydt["+ num2str(n) + "]+="+ Totaltemp[breakpos, StrLen(Totaltemp)-1] + "\r" Else print Totaltemp abort "something is strange with this string" Endif Else //Print "dydt["+ num2str(n) + "]=" + Totaltemp dydtstr += "dydt["+ num2str(n) + "]=" + Totaltemp+ "\r" Endif Endif n +=1 While (n < speciesnumber) If(Update_ChemKineticFuncs == 1) NUM_RemoveChemKineticFuncs() Execute/Q /P "EditChemKineticFuncs()" NUM_IncludeChemKineticFuncs() Endif End Function //************************************************************************************************ Macro NUM_Run_Model2() //Integrate ODEs - defined in "ChemKineticFct".ipf and output to individual Specie waves (waves made for each chemical compound in Mechanism_Complie) Silent 1 If (!exists("CurrChemKineticFct")) Abort "Must compile ODE first!" Endif String cmd, SpeciesList, CurrSpecies Variable n, speciesnumber Variable npnts //SetDataFolder root: Silent 1 PauseUpdate Variable starttime = DateTime //make species list to be used as list of input/output waves n=0 SpeciesList = "" SpeciesNumber = numpnts(Names) Do if(n == 0) SpeciesList += Names[n] else SpeciesList = SpeciesList + "," + Names[n] endif n += 1 while (n < numpnts(Names)) // set each specie wave (n = 0) to its initial concentration n=0 Do CurrSpecies = Names[n] $CurrSpecies[0] = InitialConcentrations[n] n += 1 while(n < SpeciesNumber) cmd= "IntegrateODE/M=3/F=(1+2)/U=10000/X=CalTime "+CurrChemKineticFct+","+" KK, {"+ SpeciesList + "}" //derivative functions in CurrChemKineticFct //parameter/constant coefficients wave = KK (rate constants) //xx = cal time ( x values at which to report results) //yw = NO2, O3, N2O5, etc. (individual specie waves, species specified by "SpeciciesList") - waves containing initial y[i] values //dydt = NO2, O3, N2O5, etc. - waves to receive solutions Execute cmd End Macro