What version of the code are you using?
CESM 2.1.5
Have you made any changes to files in the source tree? User_nl_cam:
! Users should add all user specific namelist changes below in the form of
! namelist_var = new_namelist_value
empty_htapes = .true.
nhtfrq(2) = -24
mfilt(2) = 365
fincl2 = 'bc_a4', 'num_a4', 'bc_a1', 'num_a1', 'T', 'Z3'
fincl3 = ' '
fincl4 = ' '
fincl5 = ' '
fincl6 = ' '
fincl7 = ' '
fincl8 = ' '
fincl9 = ' '
nhtfrq = 0, -24
ext_frc_specifier = 'bc_a4 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_bc_a4_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'bc_a4 -> /home/jjs1u25/code/spacecraft_emissions/a4_bc_spacecraft_south_pole.nc',
'so4_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a1_anthro-ene_vertical_2000climo_0.9x1.25_c20170616.nc',
'so4_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a1_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'so4_a2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a2_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_so4_a1_anthro-ene_vertical_2000climo_0.9x1.25_c20170616.nc',
'num_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_a1_so4_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_a2_so4_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a4 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_bc_a4_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'num_a4 -> /home/jjs1u25/code/spacecraft_emissions/num_a4_bc_spacecraft_south_pole.nc',
'NO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_NO2_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'SVOC -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SVOC_bb_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SO2_contvolcano_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SO2_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/stratvolc/VolcanEESMv3.10_SO2_1995-2005average_1deg_ZeroTrop_c180912.nc'
Describe every step you took leading up to the problem:
import xarray as xr
import numpy as np
import math
# Define constants
# rho = 1700.0 # kg/m3 for black carbon
rho = 2500.0 # kg/m3 for dust
# d = 0.134e-6 # m for black carbon in BC mode a4
d = 8.874e-8 # m for aitken mode dust
# mw = 12. # g/mol molecular weight for carbon (mw for C)
mw = 135.064 # molecular weight for dust (AlSiO5)
N_A = 6.02214076e23 # Avogadro's constant (matches CESM internal code)
total_mass = 10e9 # grams
molecular_mass = 12 # black carbon
avogadros_number = 6.022e23
seconds_per_year = 365.25 * 24 * 3600
radius_earth = 6.371e8
# Returns volume between latitudes phi_max and phi_min and altitudes h_min and h_max
def sector_volume(phi_max, phi_min, h_max, h_min):
phi_max = phi_max * np.pi / 180
phi_min = phi_min * np.pi / 180
a = 2*np.pi * (np.sin(phi_max) - np.sin(phi_min))
b = radius_earth**2 * (h_max-h_min) + radius_earth * (h_max**2-h_min**2) + 1/3 * (h_max**3 - h_min**3)
return a*b
# Converts molecular emission to particular emission for use in CESM
def convert(molecule_file):
M_particle = (math.pi/6 * rho * d**3)
Pre_factor = mw / M_particle
particle_file = np.asarray(molecule_file) * Pre_factor
return particle_file
ds = xr.open_dataset("/iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo_2deg/emissions-cmip6_bc_a4_aircraft_vertical_2000climo_1.9x2.5_c20200422.nc")
new_altitude_int = np.linspace(50,90,26)
new_altitude = np.zeros(len(new_altitude_int)-1)
for i in range(len(new_altitude_int)-1):
new_altitude = (new_altitude_int+new_altitude_int[i+1])/2
ds = ds.assign_coords(altitude=new_altitude)
ds = ds.assign_coords(altitude_int=new_altitude_int)
molecule_num = total_mass / molecular_mass * avogadros_number
# phi_max = 90 # polar
# phi_min = 67.26315789 # polar
phi_max = 23.6842105 # tropical
phi_min = -23.6842105 # tropical
h_max = 80.4 * 1000 * 100 # cm
h_min = 59.6 * 1000 * 100 # cm
flux = molecule_num / (seconds_per_year * sector_volume(phi_max, phi_min, h_max, h_min)) # molecules per cm3 per second
ds['emiss_aircraft'].values[:] = 0
ds['emiss_aircraft'].values[:, 6:19, 0:13, :] = flux
ds_num = ds.copy(deep=True)
ds_num['emiss_aircraft'].values[:] = 0
ds_num['emiss_aircraft'].values[:, 6:19, 0:13, :] = convert(flux)
ds.to_netcdf("dst_a2_spacecraft_south_pole.nc")
ds_num.to_netcdf("num_dst_a2_spacecraft_south_pole.nc")
Describe your problem or question:
I am trying to experiment with mesospheric emissions in WACCM. For now, I have been using the default black carbon emissions.
I would like to ask:
1) Have I introduced the black carbon emissions correctly?
2) How can I modify the size of my emitted particles? Does it suffice to simply change the diameter d in my python script? I'm not entirely sure what this does since I know CESM uses the Modal Aerosol Model, which only has 4 modes so I'm a little confused how changing particle sizes works.
3) How do I replace my black carbon with dust in my simulations? My understanding is that it's not as simple as replacing bc_a4 with dst_a2 or something.
4) How do I create a unique tracer for my newly emitted particles? Is this even possible? Becuase I know that as these particles age they will coagulate into different modes and so forth.
Thank you so much for any help in advance
CESM 2.1.5
Have you made any changes to files in the source tree? User_nl_cam:
! Users should add all user specific namelist changes below in the form of
! namelist_var = new_namelist_value
empty_htapes = .true.
nhtfrq(2) = -24
mfilt(2) = 365
fincl2 = 'bc_a4', 'num_a4', 'bc_a1', 'num_a1', 'T', 'Z3'
fincl3 = ' '
fincl4 = ' '
fincl5 = ' '
fincl6 = ' '
fincl7 = ' '
fincl8 = ' '
fincl9 = ' '
nhtfrq = 0, -24
ext_frc_specifier = 'bc_a4 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_bc_a4_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'bc_a4 -> /home/jjs1u25/code/spacecraft_emissions/a4_bc_spacecraft_south_pole.nc',
'so4_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a1_anthro-ene_vertical_2000climo_0.9x1.25_c20170616.nc',
'so4_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a1_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'so4_a2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_so4_a2_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_so4_a1_anthro-ene_vertical_2000climo_0.9x1.25_c20170616.nc',
'num_a1 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_a1_so4_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_a2_so4_contvolcano_vertical_2000climo_0.9x1.25_c20170724.nc',
'num_a4 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_num_bc_a4_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'num_a4 -> /home/jjs1u25/code/spacecraft_emissions/num_a4_bc_spacecraft_south_pole.nc',
'NO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_NO2_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'SVOC -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SVOC_bb_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SO2_contvolcano_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo/emissions-cmip6_SO2_aircraft_vertical_2000climo_0.9x1.25_c20170322.nc',
'SO2 -> /iridisfs/cesm_input_data/atm/cam/chem/stratvolc/VolcanEESMv3.10_SO2_1995-2005average_1deg_ZeroTrop_c180912.nc'
Describe every step you took leading up to the problem:
import xarray as xr
import numpy as np
import math
# Define constants
# rho = 1700.0 # kg/m3 for black carbon
rho = 2500.0 # kg/m3 for dust
# d = 0.134e-6 # m for black carbon in BC mode a4
d = 8.874e-8 # m for aitken mode dust
# mw = 12. # g/mol molecular weight for carbon (mw for C)
mw = 135.064 # molecular weight for dust (AlSiO5)
N_A = 6.02214076e23 # Avogadro's constant (matches CESM internal code)
total_mass = 10e9 # grams
molecular_mass = 12 # black carbon
avogadros_number = 6.022e23
seconds_per_year = 365.25 * 24 * 3600
radius_earth = 6.371e8
# Returns volume between latitudes phi_max and phi_min and altitudes h_min and h_max
def sector_volume(phi_max, phi_min, h_max, h_min):
phi_max = phi_max * np.pi / 180
phi_min = phi_min * np.pi / 180
a = 2*np.pi * (np.sin(phi_max) - np.sin(phi_min))
b = radius_earth**2 * (h_max-h_min) + radius_earth * (h_max**2-h_min**2) + 1/3 * (h_max**3 - h_min**3)
return a*b
# Converts molecular emission to particular emission for use in CESM
def convert(molecule_file):
M_particle = (math.pi/6 * rho * d**3)
Pre_factor = mw / M_particle
particle_file = np.asarray(molecule_file) * Pre_factor
return particle_file
ds = xr.open_dataset("/iridisfs/cesm_input_data/atm/cam/chem/emis/CMIP6_emissions_2000climo_2deg/emissions-cmip6_bc_a4_aircraft_vertical_2000climo_1.9x2.5_c20200422.nc")
new_altitude_int = np.linspace(50,90,26)
new_altitude = np.zeros(len(new_altitude_int)-1)
for i in range(len(new_altitude_int)-1):
new_altitude = (new_altitude_int+new_altitude_int[i+1])/2
ds = ds.assign_coords(altitude=new_altitude)
ds = ds.assign_coords(altitude_int=new_altitude_int)
molecule_num = total_mass / molecular_mass * avogadros_number
# phi_max = 90 # polar
# phi_min = 67.26315789 # polar
phi_max = 23.6842105 # tropical
phi_min = -23.6842105 # tropical
h_max = 80.4 * 1000 * 100 # cm
h_min = 59.6 * 1000 * 100 # cm
flux = molecule_num / (seconds_per_year * sector_volume(phi_max, phi_min, h_max, h_min)) # molecules per cm3 per second
ds['emiss_aircraft'].values[:] = 0
ds['emiss_aircraft'].values[:, 6:19, 0:13, :] = flux
ds_num = ds.copy(deep=True)
ds_num['emiss_aircraft'].values[:] = 0
ds_num['emiss_aircraft'].values[:, 6:19, 0:13, :] = convert(flux)
ds.to_netcdf("dst_a2_spacecraft_south_pole.nc")
ds_num.to_netcdf("num_dst_a2_spacecraft_south_pole.nc")
Describe your problem or question:
I am trying to experiment with mesospheric emissions in WACCM. For now, I have been using the default black carbon emissions.
I would like to ask:
1) Have I introduced the black carbon emissions correctly?
2) How can I modify the size of my emitted particles? Does it suffice to simply change the diameter d in my python script? I'm not entirely sure what this does since I know CESM uses the Modal Aerosol Model, which only has 4 modes so I'm a little confused how changing particle sizes works.
3) How do I replace my black carbon with dust in my simulations? My understanding is that it's not as simple as replacing bc_a4 with dst_a2 or something.
4) How do I create a unique tracer for my newly emitted particles? Is this even possible? Becuase I know that as these particles age they will coagulate into different modes and so forth.
Thank you so much for any help in advance