Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add new aci and aerosol activation metrics from ARM Diags v3 #679

Merged
merged 11 commits into from
May 23, 2023
Merged

Add new aci and aerosol activation metrics from ARM Diags v3 #679

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
d4df19f
add new variables needed for aci and aerosol activation metrics
chengzhuzhang Mar 17, 2023
b005733
update time_range
chengzhuzhang Apr 6, 2023
b0b0f90
code change for arm diags v3 data
chengzhuzhang Apr 13, 2023
25a9402
fix viwer; fix format
chengzhuzhang Apr 14, 2023
c7cbd32
add ENA site
chengzhuzhang Apr 14, 2023
4fc8b3c
add aerosol activation driver and plotting
chengzhuzhang Apr 17, 2023
a438eca
add aerosol activation viewer
chengzhuzhang Apr 17, 2023
f5128a8
add aerosol annual cycle
chengzhuzhang Apr 18, 2023
227561f
fix units
chengzhuzhang Apr 18, 2023
dde1ef7
Merge branch 'main' into arm_diags_v3
chengzhuzhang May 23, 2023
ab041b4
fix integration tests
chengzhuzhang May 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 68 additions & 1 deletion e3sm_diags/derivations/acme.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
from cdms2.fvariable import FileVariable

AVOGADOR_CONS = 6.022e23
AIR_DENS = 1.225 # standard air density 1.225kg/m3


def rename(new_name):
Expand Down Expand Up @@ -147,7 +148,23 @@ def molec_convert_units(var, molar_weight):
# Convert molec/cm2/s to kg/m2/s
if var.units == "molec/cm2/s":
var = var / AVOGADOR_CONS * molar_weight * 10.0
var.units == "kg/m2/s"
var.units = "kg/m2/s"
return var


def a_num_sum(var):
# Calculate: total aerosol number concentration (#/cm3)
var = var * AIR_DENS / 1e6
var.units = "/cm3"
var.long_name = "aerosol number concentration"
return var


def so4_mass_sum(var):
# Calculate: SO4 mass conc. (ng/m3) (< 1um)
var = var * AIR_DENS * 1e9
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lannyzxj Xiaojian,this formula converts so4 mass conc. from (kg/kg) to ug/m3 which gives comparable results to models. But in your version, the results are in ug/cm3. I might get something wrong here..would you help me double check.

Copy link

@lannyzxj lannyzxj Apr 18, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Dear Dr. Zhang @chengzhuzhang ,

Thanks for letting me know. I have cross-compared your plots and mine offline, the plots look good.
My units are also ug/m3 for the sulfate output in my plot, which looks comparable to yours.

I double-checked my processing code for the model so4,

it was first converted to ng/m3 after reading the data:
"
iso4_Acc_tmp= iso4_Acc_tmp * 1e12 * 1.225 #kg/kg to ng/kg then to ng/m3
"

and then further divided by 1000 in order to convert to ug/m3 in the data output part of the code:
"
#SO4
ivar=np.reshape(so4_all,(numday*24),order='C')
ivar=ivar/1000. #convert to ug/m3
Mvar=MV2.masked_array(ivar,mask=[np.isnan(ivar)],fill_value=-9999.)
Mvar.id="sulfate"
Mvar.long_name='Sulfate (<1um) Mass Conc.'
Mvar.units='ug/m3'
Mvar.missing_value=-9999.
Mvar.setAxis(0,hrtime_axes)
fout.write(Mvar)
"
The above processes should be equivalent to your "var = var * AIR_DENS * 1e9".

Would you point me to where my version states the unit as ug/cm3? That could be a typo in stating the unit, the actual calculation should indeed yield ug/m3.

Best,
Xiaojian

Copy link
Contributor Author

@chengzhuzhang chengzhuzhang Apr 18, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for your prompt reply! I missed the second units conversion step in your script when writing out files! It is good to know that our results are consistent.
I was looking at a version of ppt for ARM diags v3(for verification purpose), which includes aerosol diurnal cycle plots that has units in ug/cm3, which could be a typo. And I also confirmed that the ARM Diags standalone run provide the correct units!

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good to know that! Thank you!

var.units = "\u03bcg/m3"
var.long_name = "SO4 mass conc."
return var


Expand Down Expand Up @@ -1742,6 +1759,56 @@ def cosp_histogram_standardize(cld: "FileVariable"):
(("Mass_pom",), rename),
]
),
# total aerosol number concentration (#/CC)
"a_num": OrderedDict(
[
(("cpc",), rename),
# Aerosol concentration from Aitken, Accumu., and Coarse mode
(
(
"num_a1",
"num_a2",
"num_a3",
),
lambda a1, a2, a3: a_num_sum(a1 + a2 + a3),
),
]
),
# total so4 mass concentration (ng/m3)
"so4_mass": OrderedDict(
[
(("sulfate",), rename),
# Aerosol concentration from Aitken, Accumu., and Coarse mode
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note: this only takes Aitken and Accumu. mode. remove Coarse in the comment

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lannyzxj Cheng and I are working on clearing up some documentation in preparation for E3SM Tutorial. When working on ARM Diags, I noticed that when computing so4 mass concentration, only Aitken and Accumu. mode are taking account. Could you remind me why is the case? Is it to make the model data comparable to obs?

(
(
"so4_a1",
"so4_a2",
),
lambda a1, a2: so4_mass_sum(a1 + a2),
),
]
),
# CCN 0.1%SS concentration (1/CC)
"ccn01": OrderedDict(
[
(("ccn01",), rename),
(("CCN3",), rename),
]
),
# CCN 0.2%SS concentration (1/CC)
"ccn02": OrderedDict(
[
(("ccn02",), rename),
(("CCN4",), rename),
]
),
# CCN 0.5%SS concentration (1/CC)
"ccn05": OrderedDict(
[
(("ccn05",), rename),
(("CCN5",), rename),
]
),
# Land variables
"SOILWATER_10CM": OrderedDict([(("mrsos",), rename)]),
"SOILWATER_SUM": OrderedDict([(("mrso",), rename)]),
Expand Down
5 changes: 3 additions & 2 deletions e3sm_diags/derivations/default_regions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,9 +174,10 @@
points_specs = {
# ARM sites coordinates, select nearest grid poit to ARM site coordinates
# Each point is supplied with [latitude, longitude ,select method, description of the point]
"sgp": [36.4, -97.5, "cob", "97.5W 36.4N Oklahoma ARM"],
"nsa": [71.3, -156.6, "cob", "156.6W 71.3N Barrow ARM"],
"sgpc1": [36.4, -97.5, "cob", "97.5W 36.4N Oklahoma ARM"],
"nsac1": [71.3, -156.6, "cob", "156.6W 71.3N Barrow ARM"],
"twpc1": [-2.1, 147.4, "cob", "147.4E 2.1S Manus ARM"],
"twpc2": [-0.5, 166.9, "cob", "166.9E 0.5S Nauru ARM"],
"twpc3": [-12.4, 130.9, "cob", "130.9E 12.4S Darwin ARM"],
"enac1": [39.1, -28.0, "cob", "28.0E 39.1N Graciosa Island ARM"],
}
214 changes: 184 additions & 30 deletions e3sm_diags/driver/arm_diags_driver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ def run_diag_diurnal_cycle(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
refs = []

if "armdiags" in ref_name:
if region != "sgp":
if region != "sgpc1":
msg = "Diurnal cycle of {} at Site: {} is not supported yet".format(
region, var
)
Expand All @@ -114,7 +114,9 @@ def run_diag_diurnal_cycle(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
) = utils.diurnal_cycle.composite_diurnal_cycle(
ref, season, fft=False
)
ref.long_name = ref.standard_name
if hasattr(ref, "standard_name"):
ref.long_name = ref.standard_name

ref = ref_diurnal

else:
Expand Down Expand Up @@ -195,19 +197,12 @@ def run_diag_diurnal_cycle_zt(parameter: ARMDiagsParameter) -> ARMDiagsParameter
refs = []

if "armdiags" in ref_name:
if region == "sgp":
ref_file_name = "sgparmdiagsmondiurnalclimC1.c1.nc"
elif region == "nsa":
ref_file_name = (
region[:3] + "armdiagsmondiurnalclim" + "C1.c1.nc"
)
else:
ref_file_name = (
region[:3]
+ "armdiagsmondiurnalclim"
+ region[3:5].upper()
+ ".c1.nc"
)
ref_file_name = (
region[:3]
+ "armdiagsmondiurnalclim"
+ region[3:5].upper()
+ ".c1.nc"
)

ref_file = os.path.join(ref_path, ref_file_name)
ref_data = cdms2.open(ref_file)
Expand Down Expand Up @@ -303,25 +298,27 @@ def run_diag_annual_cycle(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
refs = []

if "armdiags" in ref_name:
if region == "sgp":
ref_file = os.path.join(ref_path, "sgparmdiagsmonC1.c1.nc")
elif region == "nsa":
ref_file = os.path.join(
ref_path, region[:3] + "armdiagsmonclim" + "C1.c1.nc"
)
else:
ref_file = os.path.join(
ref_path,
region[:3]
+ "armdiagsmonclim"
+ region[3:5].upper()
+ ".c1.nc",
)
# in ARM Diags v2 only sgp site has monthly time series other sites have annual cycle , i.e. 12 time points.
# if "sgp" in region:
# ref_file = os.path.join(ref_path, "sgparmdiagsmonC1.c1.nc")
# else:
# ref_file = os.path.join(
# ref_path,
# region[:3]
# + "armdiagsmonclim"
# + region[3:5].upper()
# + ".c1.nc",
# )
ref_file = os.path.join(
ref_path,
region[:3] + "armdiagsmon" + region[3:5].upper() + ".c1.nc",
)
ref_data = cdms2.open(ref_file)
vars_funcs = get_vars_funcs_for_derived_var(ref_data, var)
target_var = list(vars_funcs.keys())[0][0]
ref_var = ref_data(target_var)
ref_var.long_name = ref_var.standard_name
if hasattr(ref_var, "standard_name"):
ref_var.long_name = ref_var.standard_name
ref = vars_funcs[(target_var,)](utils.climo.climo(ref_var, season))

else:
Expand Down Expand Up @@ -413,6 +410,159 @@ def run_diag_convection_onset(parameter: ARMDiagsParameter) -> ARMDiagsParameter
return parameter


def run_diag_aerosol_activation(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
regions = parameter.regions
ref_name = parameter.ref_name
ref_path = parameter.reference_data_path
variables = parameter.variables
# Read in observation data

for region in regions:
# The regions that are supported are in e3sm_diags/derivations/default_regions.py
# You can add your own if it's not in there.
logger.info("Selected region: {}".format(region))
# Possible variables are ccn01, ccn02, ccn05
for variable in variables:

test_data = utils.dataset.Dataset(parameter, test=True)

test_a_num = test_data.get_timeseries_variable("a_num", single_point=True)[
:,
-1,
].filled(fill_value=np.nan)
test_ccn = test_data.get_timeseries_variable(variable, single_point=True)[
:,
-1,
].filled(fill_value=np.nan)

# Get the name of the data, appended with the years averaged.
parameter.test_name_yrs = utils.general.get_name_and_yrs(
parameter, test_data
)

if "armdiags" in ref_name:

ref_file = os.path.join(
ref_path,
region[:3] + "armdiagsaciactivate" + region[3:5].upper() + ".c1.nc",
)
ref_data = cdms2.open(ref_file)
ref_a_num = ref_data("cpc_bulk").filled(fill_value=np.nan)
ref_ccn = ref_data(f"{variable}_bulk").filled(fill_value=np.nan)

else:
ref_data = utils.dataset.Dataset(parameter, ref=True)
ref_a_num = test_data.get_timeseries_variable(
"a_num", single_point=True
)[
:,
-1,
].filled(
fill_value=np.nan
)
ref_ccn = test_data.get_timeseries_variable(
variable, single_point=True
)[
:,
-1,
].filled(
fill_value=np.nan
)

parameter.output_file = "-".join(
[ref_name, "aerosol-activation", region, variable]
)
arm_diags_plot.plot_aerosol_activation(
test_a_num, test_ccn, ref_a_num, ref_ccn, parameter, region, variable
)

return parameter


def run_diag_annual_cycle_aerosol(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
variables = parameter.variables
regions = parameter.regions
ref_name = parameter.ref_name
ref_path = parameter.reference_data_path

seasons = ["ANNUALCYCLE"]

for region in regions:
# The regions that are supported are in e3sm_diags/derivations/default_regions.py
# You can add your own if it's not in there.
logger.info("Selected region: {}".format(region))
vars_to_data = collections.OrderedDict()

for season in seasons:
logger.info("Season: {}".format(season))
for var in variables:
logger.info("Variable: {}".format(var))
test_data = utils.dataset.Dataset(parameter, test=True)
test = test_data.get_climo_variable(var, season)[:, -1]

parameter.viewer_descr[var] = getattr(test, "long_name", var)
# Get the name of the data, appended with the years averaged.
parameter.test_name_yrs = utils.general.get_name_and_yrs(
parameter, test_data
)
parameter.var_name = getattr(test, "long_name", var)
parameter.var_units = getattr(test, "units", var)

refs = []

if "armdiags" in ref_name:
ref_file = os.path.join(
ref_path,
region[:3] + "armdiagsaciclim" + region[3:5].upper() + ".c1.nc",
)
ref_data = cdms2.open(ref_file)
vars_funcs = get_vars_funcs_for_derived_var(ref_data, var)
target_var = list(vars_funcs.keys())[0][0]
ref_var = ref_data(target_var)[:, 0] # 0 mean; 1 standard devation
if hasattr(ref_var, "standard_name"):
ref_var.long_name = ref_var.standard_name
ref = vars_funcs[(target_var,)](utils.climo.climo(ref_var, season))

else:
ref_data = utils.dataset.Dataset(parameter, ref=True)
ref = ref_data.get_climo_variable(var, season)[:, -1]
ref_domain = utils.general.select_point(region, ref)
ref.ref_name = ref_name
refs.append(ref_domain)

test_domain = utils.general.select_point(region, test)

metrics_dict = create_metrics(test_domain, ref_domain)

result = RefsTestMetrics(
test=test_domain, refs=refs, metrics=metrics_dict, misc=None
)
vars_to_data[season] = result
# Saving the metrics as a json.
metrics_dict["unit"] = test.units
print(parameter.var_units, test.units)
parameter.output_file = "-".join(
[ref_name, var, season, "aerosol", region]
)
fnm = os.path.join(
utils.general.get_output_dir(parameter.current_set, parameter),
parameter.output_file + ".json",
)
with open(fnm, "w") as outfile:
json.dump(metrics_dict, outfile)
# Get the filename that the user has passed in and display that.
fnm = os.path.join(
utils.general.get_output_dir(parameter.current_set, parameter),
parameter.output_file + ".json",
)
logger.info(f"Metrics saved in: {fnm}")

if season == "ANNUALCYCLE":
arm_diags_plot.plot_annual_cycle(var, vars_to_data[season], parameter)

return parameter


def run_diag_pdf_daily(parameter: ARMDiagsParameter):
logger.info("'run_diag_pdf_daily' is not yet implemented.")

Expand All @@ -429,5 +579,9 @@ def run_diag(parameter: ARMDiagsParameter) -> ARMDiagsParameter:
return run_diag_pdf_daily(parameter)
elif parameter.diags_set == "convection_onset":
return run_diag_convection_onset(parameter)
elif parameter.diags_set == "aerosol_activation":
return run_diag_aerosol_activation(parameter)
if parameter.diags_set == "annual_cycle_aerosol":
return run_diag_annual_cycle_aerosol(parameter)
else:
raise Exception("Invalid diags_set={}".format(parameter.diags_set))
6 changes: 6 additions & 0 deletions e3sm_diags/driver/default_diags/arm_diags_model_vs_model.cfg
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,3 +20,9 @@ variables = ["CLOUD"]
sets = ["arm_diags"]
diags_set = "convection_onset"
regions = ["sgp", "twpc1", "twpc2", "twpc3"]

[#]
sets = ["arm_diags"]
diags_set = "aerosol_activation"
regions = ["sgpc1", "enac1"]
variables = ["ccn02", "ccn05"]
Loading