bug fix and new spin up option

in ui_check, tiling of seedx was leading to the first dimension being transformed into the second dimension, so this was fixed. Users can now set a new spin up option, where the spin-up loops over the total simulation time until total RO2 concentrations at the end of the simulation vary by less than 10%.

PyCHAM/dydt_rec.py

@@ -24,7 +24,7 @@
 # changes due to gas-phase photochemistry and partitioning are included; 
 # generated in init_conc and treats loss from gas-phase as negative
 
-# File Created at 2024-11-12 16:11:40.155400
+# File Created at 2024-11-15 17:34:41.360613
 
 import numpy as np 
 

PyCHAM/eqn_interr.py

@@ -29,7 +29,7 @@
 import openbabel.pybel as pybel
 import sys
 
-def eqn_interr(num_sb, self):
+def eqn_interr(num_sb, erf, err_mess, self):
 
 	# inputs: ----------------------------------------------------------------------------
 	# self.eqn_num - number of equations
@@ -554,8 +554,14 @@ def eqn_interr(num_sb, self):
 			# account for this component
 			comp_num += 1
 			# convert MCM chemical names to SMILES
-			# index where xml file name matches reaction component name
-			name_indx = self.comp_xmlname.index(sname)
+			# index where xml file name matches component name
+			try:
+				name_indx = self.comp_xmlname.index(sname)
+			except:
+				erf = 1 # raise error
+				err_mess = str('Error: inside eqn_interr, seed component ' + 
+				str(sname) + ' not found in xml file.')
+				return(comp_list, Pybel_objects, comp_num, erf, err_mess, self)
 
 			# get SMILES string
 			name_SMILE = self.comp_smil[name_indx]
@@ -995,4 +1001,4 @@ def eqn_interr(num_sb, self):
 			self.obs_comp_i[i] = self.comp_namelist.index(
 			self.obs_comp[i])
 
-	return(comp_list, Pybel_objects, comp_num, self)
+	return(comp_list, Pybel_objects, comp_num, erf, err_mess, self)

PyCHAM/eqn_pars.py

@@ -96,8 +96,11 @@ def extr_mech(int_tol, num_sb, drh_str, erh_str, self):
 	[err_mess_new, self] = xml_interr.xml_interr(self)
 
 	# get equation information for chemical reactions
-	[comp_list, Pybel_objects, comp_num, self] = eqn_interr.eqn_interr(
-		num_sb, self)                                   
+	[comp_list, Pybel_objects, comp_num, erf, err_mess, self] = eqn_interr.eqn_interr(
+		num_sb, erf, err_mess, self)    
+
+	if (erf != 0): # exit and display error, if it exists
+		return([], [], [], [], erf, err_mess, self)                               
 
 	# prepare aqueous-phase and surface (e.g. wall) reaction matrices for applying 
 	# to reaction rate calculation
@@ -110,7 +113,11 @@ def extr_mech(int_tol, num_sb, drh_str, erh_str, self):
 	# not turned on then raise an error
 	if (self.eqn_num[1] > 0 and num_sb-self.wall_on == 0):
 		erf = 1 # raise error
-		err_mess = str('Error: ' +str(self.eqn_num[1]) + ' particle-phase reactions were registered (from the chemical scheme input file), but no particle size bins have been invoked (number_size_bins variable in the model variables input file). Please ensure consistency. (message generated by eqn_pars.py module)')
+		err_mess = str('Error: ' + str(self.eqn_num[1]) + ' particle-phase ' +
+		'reactions were registered (from the chemical scheme input file), but ' +
+		'no particle size bins have been invoked (number_size_bins variable in ' +
+		'the model variables input file). Please ensure consistency. (message ' +
+		'generated by eqn_pars.py module)')
 
 	# get index of components with continuous influx/concentration -----------
 	# empty array for storing index of components with constant influx

PyCHAM/err_log.txt

@@ -1 +1 @@
-/Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_Exp2p1
+/Users/user/Documents/GitHub/PyCHAM/PyCHAM/output/mcmAP_scheme/20Nov_0900_0959_Base

PyCHAM/hyst_eq.py

@@ -21,7 +21,7 @@
 ##########################################################################################
 '''solution of deliquescence and efflorescence RH, generated by eqn_pars.py in fully functioning mode, or by ui_check.py in testing mode'''
 # module to estimate deliquescence and efflorescence relative humidity as a function of temperature
-# File Created at 2024-11-12 16:26:35.051795
+# File Created at 2024-11-15 17:34:41.373159
 
 # function for deliquescence
 def drh(TEMP):

PyCHAM/input/Baker2024/mcmAP_PRAMAP_scheme.kpp

@@ -1426,6 +1426,10 @@ J(24) = 4.365E-05*(cos(zenith)**1.089)*exp(-0.323*(1./cos(zenith))) {MCM J=56.}
  {881 } HCOCO = HCOCO3 : 5.00E-12*O2*3.2*EXP(-550./TEMP) ;
  {882 } HCOCO = CO + OH : 5.00E-12*O2*3.2*(1.-EXP(-550./TEMP)) ;
  {883 } HMVKANO3 + hv = HMVKAO + NO2 : J(24)*0.91 ;
+# background reactivity used for OH in supplemntary material of Baker et al. 2024 
+# (10.5194/acp-24-4789-2024-supplement
+ {884 } OH = : 5. ;
+
 # PRAM (Peroxy Radical Autoxidation Mechanism) in Kinetic Pre Processor (KPP) format.
 # This file include the complete PRAM version used by Roldin et al., (NCOMMS 2019), 208 species and 1773 reactions.
 # Combine all PRAM reactions and species with a suitable Master Chemical Mechanism (MCM) KPP file.

PyCHAM/input/Baker2024/mod_var_Exp2p1.txt

@@ -0,0 +1,76 @@
+res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_PRAMAP_scheme/Baker_2024_Exp2p1_zpc1e-4_it1e-3
+chem_sch_name = mcmAP_PRAMAP_scheme.kpp
+xml_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_input/MCM_PRAM_xml.xml
+total_model_time = 8.64e4
+update_step = 3.6e2
+recording_time_step = 3.6e3
+temperature = 293.15
+tempt = 0.
+p_init = 101325.0
+rh = 0.52
+rht = 0
+# UVC lights on
+light_status = J_values.xlsx
+# note that the photon intensity for 254 nm in this file
+# has been tuned to give a photolysis rate of J(O1D) 
+# i.e. O3 + hv = O1D of 
+# 8x10-4 /s when the exposed light is at a minimum of 35 %,
+# which Baker et al. 2024 section 3.1 says is the minimum
+# J(O1D) rate. 
+#act_flux_file = SAPHIR_STAR_UVC_Actinic_Flux_Spectrum.csv
+const_comp = H2O
+# initial concentrations of components (ppb)
+Comp0 = O3, CO, APINENE, OH
+C0 = 29., 0., 7.5, 3.4e-5
+# continuously influxed components
+cont_infl = O3, APINENE, CO
+cont_infl_t = 0., 4.32e4
+# based on Figure 2 of Baker et al. 2024,
+# we need to tune the Cinfl values so that
+# in both the with and without CO stages of
+# experiment O3 has a mixing ratio of 25 ppb,
+# alpha-pinene has a mixing ratio of 5 ppb and
+# CO has a mixing ratio of 0 ppm in the first
+# stage and a mixing ratio of 2.5 ppm
+# these give correct results for APINENE, O3 and CO
+# when light intensity
+# in photofile set to 2.3e14 photons/cm2/s, but
+# OH comes out around a factor of two higher than Figure 2
+# observations
+# lights tuned down for O3 and OH agreement
+#Cinfl = 0.0076, 0.0089; 0.0028, 0.0028; 0., 0.69
+# lights set to give same J(O1D) as in Baker paper
+#Cinfl = 0.0083, 0.011; 0.0033, 0.0033; 0., 0.67
+Cinfl = 0.011, 0.016; 0.0033, 0.0033; 0., 0.67
+# residence time of 63 minutes
+dil_fac = 0.000265
+wall_on = 1
+# gas-wall partitioning
+# 1.1e-2 /s for all species from Ehn et al. 2013, 
+# from S3 of supplement of Baker et al. 2024:
+# 5.88e-3 /s for RO2, OH loss rate
+# (5 /s) is actually due to a 'background reactivity', so
+# it is stated in the chemical scheme
+# and HO2 wall loss rate of 
+mass_trans_coeff = RO2_wall1_5.88e-3; HO2_wall1_2e-2
+# once these components go to wall do not allow evaporation from wall
+vol_Comp = RO2_wall1, HO2_wall1
+volP = 0., 0.
+eff_abs_wall_massC = 1.e0
+chem_scheme_markers = {, RO2, +, C(ind_, ), , &, , , :, }, ;,
+# set the vapour pressure estimation method for HOMs
+HOMs_vp_method = Mohr2019
+# only need to use the below variables for the seeded case
+number_size_bins = 1
+pconc = 2.3e1
+pconct = 0.
+pcont = 1
+mean_rad = 3.0e-2
+std = 1.5
+Vwat_inc = 2
+# 125 nm radius
+upper_part_size = 1.25e-1
+z_prt_coeff = 1.e-4
+int_tol = 1.e-3, 1.e-3
+tracked_comp = O3, OH, HO2, APINENE
+pars_skip = 0

PyCHAM/input/Baker2024/mv_Exp1.txt

@@ -1,5 +1,5 @@
-res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_unseeded_JO1Dfixed
-chem_sch_name = mcmAP_scheme.kpp
+res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_PRAMAP_scheme/Baker_2024_Exp1_defit
+chem_sch_name = mcmAP_PRAMAP_scheme.kpp
 xml_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_input/MCM_PRAM_xml.xml
 total_model_time = 8.64e4
 update_step = 3.6e2

PyCHAM/input/Baker2024/mv_Exp2p2.txt

@@ -0,0 +1,76 @@
+res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_Exp2p2_it1e-4
+chem_sch_name = mcmAP_scheme.kpp
+xml_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_input/MCM_PRAM_xml.xml
+total_model_time = 8.64e4
+update_step = 3.6e2
+recording_time_step = 3.6e3
+temperature = 293.15
+tempt = 0.
+p_init = 101325.0
+rh = 0.52
+rht = 0
+# UVC lights on
+light_status = J_values.xlsx
+# note that the photon intensity for 254 nm in this file
+# has been tuned to give a photolysis rate of J(O1D) 
+# i.e. O3 + hv = O1D of 
+# 8x10-4 /s when the exposed light is at a minimum of 35 %,
+# which Baker et al. 2024 section 3.1 says is the minimum
+# J(O1D) rate. 
+#act_flux_file = SAPHIR_STAR_UVC_Actinic_Flux_Spectrum.csv
+const_comp = H2O
+# initial concentrations of components (ppb)
+Comp0 = O3, CO, APINENE, OH
+C0 = 29., 0., 7.5, 3.4e-5
+# continuously influxed components
+cont_infl = O3, APINENE, CO
+cont_infl_t = 0., 4.32e4
+# based on Figure 2 of Baker et al. 2024,
+# we need to tune the Cinfl values so that
+# in both the with and without CO stages of
+# experiment O3 has a mixing ratio of 25 ppb,
+# alpha-pinene has a mixing ratio of 5 ppb and
+# CO has a mixing ratio of 0 ppm in the first
+# stage and a mixing ratio of 2.5 ppm
+# these give correct results for APINENE, O3 and CO
+# when light intensity
+# in photofile set to 2.3e14 photons/cm2/s, but
+# OH comes out around a factor of two higher than Figure 2
+# observations
+# lights tuned down for O3 and OH agreement
+#Cinfl = 0.0076, 0.0089; 0.0028, 0.0028; 0., 0.69
+# lights set to give same J(O1D) as in Baker paper
+#Cinfl = 0.0083, 0.011; 0.0033, 0.0033; 0., 0.67
+Cinfl = 0.011, 0.015; 0.0033, 0.0033; 0., 0.67
+# residence time of 63 minutes
+dil_fac = 0.000265
+wall_on = 1
+# gas-wall partitioning
+# 1.1e-2 /s for all species from Ehn et al. 2013, 
+# from S3 of supplement of Baker et al. 2024:
+# 5.88e-3 /s for RO2, OH loss rate
+# (5 /s) is actually due to a 'background reactivity', so
+# it is stated in the chemical scheme
+# and HO2 wall loss rate of 
+mass_trans_coeff = RO2_wall1_5.88e-3; HO2_wall1_2e-2
+# once these components go to wall do not allow evaporation from wall
+vol_Comp = RO2_wall1, HO2_wall1
+volP = 0., 0.
+eff_abs_wall_massC = 1.e0
+chem_scheme_markers = {, RO2, +, C(ind_, ), , &, , , :, }, ;,
+# set the vapour pressure estimation method for HOMs
+HOMs_vp_method = Mohr2019
+# only need to use the below variables for the seeded case
+number_size_bins = 0
+#pconc = 1.8e1
+#pconct = 0.
+#pcont = 1
+#mean_rad = 3.0e-2
+#std = 1.5
+#Vwat_inc = 2
+# 125 nm radius
+#upper_part_size = 1.25e-1
+#z_prt_coeff = 1.e-4
+int_tol = 1.e-4, 1.e-4
+tracked_comp = O3, OH, HO2, APINENE
+pars_skip = 1

PyCHAM/input/Baker2024/model_var_Exp2p1.txt → PyCHAM/input/Baker2024/mv_Exp3.txt

@@ -1,4 +1,4 @@
-res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_Exp2p1
+res_file_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_Exp3_zpc1e-4_it1e-3
 chem_sch_name = mcmAP_scheme.kpp
 xml_name = /Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_input/MCM_PRAM_xml.xml
 total_model_time = 8.64e4
@@ -39,8 +39,9 @@ cont_infl_t = 0., 4.32e4
 # observations
 # lights tuned down for O3 and OH agreement
 #Cinfl = 0.0076, 0.0089; 0.0028, 0.0028; 0., 0.69
-# lights set to give same J(O1D) as in Baker paper 
-Cinfl = 0.0083, 0.011; 0.0033, 0.0033; 0., 0.67
+# lights set to give same J(O1D) as in Baker paper
+#Cinfl = 0.0083, 0.011; 0.0033, 0.0033; 0., 0.67
+Cinfl = 0.010, 0.014; 0.0033, 0.0033; 0., 0.67
 # residence time of 63 minutes
 dil_fac = 0.000265
 wall_on = 1
@@ -56,17 +57,21 @@ mass_trans_coeff = RO2_wall1_5.88e-3; HO2_wall1_2e-2
 vol_Comp = RO2_wall1, HO2_wall1
 volP = 0., 0.
 eff_abs_wall_massC = 1.e0
-number_size_bins = 1
 chem_scheme_markers = {, RO2, +, C(ind_, ), , &, , , :, }, ;,
 # set the vapour pressure estimation method for HOMs
 HOMs_vp_method = Mohr2019
 # only need to use the below variables for the seeded case
-pconc = 4.e1
+number_size_bins = 1
+# 2.3e1
+pconc = 1.8e1
 pconct = 0.
-pcont = 1 
-mean_rad = 3.e-2
+pcont = 1
+mean_rad = 3.0e-2
 std = 1.5
-Vwat_inc = 0
+Vwat_inc = 2
 # 125 nm radius
 upper_part_size = 1.25e-1
-tracked_comp = O3, OH, HO2, APINENE
+z_prt_coeff = 1.e-4
+int_tol = 1.e-3, 1.e-3
+tracked_comp = O3, OH, HO2, APINENE
+pars_skip = 1