1+ # NOTE: This module uses the deprecated e3sm_diags.driver.utils.dataset.Dataset
2+ # class, which was replaced by e3sm_diags.driver.utils.dataset_xr.Dataset.
13import e3sm_diags
24from e3sm_diags .driver import utils
35import cdms2
1214
1315
1416def global_integral (var , area_m2 ):
15- """ Compute global integral of 2 dimentional properties"""
16- return numpy .sum (numpy .sum (abs (var )* area_m2 ,axis = 0 ), axis = 0 )
17+ """Compute global integral of 2 dimentional properties"""
18+ return numpy .sum (numpy .sum (abs (var ) * area_m2 , axis = 0 ), axis = 0 )
19+
1720
1821def calc_column_integral (data , aerosol , season ):
19- """ Calculate column integrated mass """
22+ """Calculate column integrated mass"""
2023
2124 # take aerosol and change it to the appropriate string
2225 # ncl -> SEASALT, dst -> DUST, rest1 -> REST1
@@ -32,129 +35,129 @@ def calc_column_integral(data, aerosol, season):
3235 burden = data .get_climo_variable (f"ABURDEN{ aerosol_name } " , season )
3336 except RuntimeError :
3437 # if not, use the Mass_ terms and integrate over the column
35- mass = data .get_climo_variable (f' Mass_{ aerosol } ' , season )
38+ mass = data .get_climo_variable (f" Mass_{ aerosol } " , season )
3639 hyai , hybi , ps = data .get_extra_variables_only (
37- f' Mass_{ aerosol } ' , season , extra_vars = ["hyai" , "hybi" , "PS" ]
40+ f" Mass_{ aerosol } " , season , extra_vars = ["hyai" , "hybi" , "PS" ]
3841 )
3942
4043 p0 = 100000.0 # Pa
41- ps = ps # Pa
42- pressure_levs = cdutil .vertical .reconstructPressureFromHybrid (ps , hyai , hybi , p0 )
44+ ps = ps # Pa
45+ pressure_levs = cdutil .vertical .reconstructPressureFromHybrid (
46+ ps , hyai , hybi , p0
47+ )
4348
44- #(72,lat,lon)
45- delta_p = numpy .diff (pressure_levs ,axis = 0 )
46- mass_3d = mass * delta_p / 9.8 # mass density * mass air kg/m2
47- burden = numpy .nansum (mass_3d ,axis = 0 ) # kg/m2
49+ # (72,lat,lon)
50+ delta_p = numpy .diff (pressure_levs , axis = 0 )
51+ mass_3d = mass * delta_p / 9.8 # mass density * mass air kg/m2
52+ burden = numpy .nansum (mass_3d , axis = 0 ) # kg/m2
4853 return burden
49-
54+
55+
5056def generate_metrics_dic (data , aerosol , season ):
5157 metrics_dict = {}
52- wetdep = data .get_climo_variable (f'{ aerosol } _SFWET' , season )
53- drydep = data .get_climo_variable (f'{ aerosol } _DDF' , season )
54- srfemis = data .get_climo_variable (f'SF{ aerosol } ' , season )
55- area = data .get_extra_variables_only (
56- f'{ aerosol } _DDF' , season , extra_vars = ["area" ]
57- )
58+ wetdep = data .get_climo_variable (f"{ aerosol } _SFWET" , season )
59+ drydep = data .get_climo_variable (f"{ aerosol } _DDF" , season )
60+ srfemis = data .get_climo_variable (f"SF{ aerosol } " , season )
61+ area = data .get_extra_variables_only (f"{ aerosol } _DDF" , season , extra_vars = ["area" ])
5862 area_m2 = area * REARTH ** 2
5963
6064 burden = calc_column_integral (data , aerosol , season )
61- burden_total = global_integral (burden , area_m2 )* 1e-9 # kg to Tg
62- print (f' { aerosol } Burden (Tg): ' , f' { burden_total :.3f} ' )
63- sink = global_integral ((drydep - wetdep ),area_m2 )* UNITS_CONV
64- drydep = global_integral (drydep ,area_m2 )* UNITS_CONV
65- wetdep = global_integral (wetdep ,area_m2 )* UNITS_CONV
66- srfemis = global_integral (srfemis ,area_m2 )* UNITS_CONV
67- print (f' { aerosol } Sink (Tg/year): ' , f' { sink :.3f} ' )
68- print (f' { aerosol } Lifetime (days): ' , f' { burden_total / sink * 365 :.3f} ' )
65+ burden_total = global_integral (burden , area_m2 ) * 1e-9 # kg to Tg
66+ print (f" { aerosol } Burden (Tg): " , f" { burden_total :.3f} " )
67+ sink = global_integral ((drydep - wetdep ), area_m2 ) * UNITS_CONV
68+ drydep = global_integral (drydep , area_m2 ) * UNITS_CONV
69+ wetdep = global_integral (wetdep , area_m2 ) * UNITS_CONV
70+ srfemis = global_integral (srfemis , area_m2 ) * UNITS_CONV
71+ print (f" { aerosol } Sink (Tg/year): " , f" { sink :.3f} " )
72+ print (f" { aerosol } Lifetime (days): " , f" { burden_total / sink * 365 :.3f} " )
6973 metrics_dict = {
70- "Surface Emission (Tg/yr)" : f' { srfemis :.3f} ' ,
71- "Sink (Tg/yr)" : f' { sink :.3f} ' ,
72- "Dry Deposition (Tg/yr)" : f' { drydep :.3f} ' ,
73- "Wet Deposition (Tg/yr)" : f' { wetdep :.3f} ' ,
74- "Burden (Tg)" : f' { burden_total :.3f} ' ,
75- "Lifetime (Days)" : f' { burden_total / sink * 365 :.3f} ' ,
74+ "Surface Emission (Tg/yr)" : f" { srfemis :.3f} " ,
75+ "Sink (Tg/yr)" : f" { sink :.3f} " ,
76+ "Dry Deposition (Tg/yr)" : f" { drydep :.3f} " ,
77+ "Wet Deposition (Tg/yr)" : f" { wetdep :.3f} " ,
78+ "Burden (Tg)" : f" { burden_total :.3f} " ,
79+ "Lifetime (Days)" : f" { burden_total / sink * 365 :.3f} " ,
7680 }
7781 return metrics_dict
7882
83+
7984param = CoreParameter ()
80- param .test_name = ' v2.LR.historical_0101'
81- param .test_name = ' F2010.PD.NGD_v3atm.0096484.compy'
82- param .test_data_path = ' /Users/zhang40/Documents/ACME_simulations/'
83- param .test_data_path = ' /compyfs/mahf708/E3SMv3_dev/F2010.PD.NGD_v3atm.0096484.compy/post/atm/180x360_aave/clim/10yr'
85+ param .test_name = " v2.LR.historical_0101"
86+ param .test_name = " F2010.PD.NGD_v3atm.0096484.compy"
87+ param .test_data_path = " /Users/zhang40/Documents/ACME_simulations/"
88+ param .test_data_path = " /compyfs/mahf708/E3SMv3_dev/F2010.PD.NGD_v3atm.0096484.compy/post/atm/180x360_aave/clim/10yr"
8489test_data = utils .dataset .Dataset (param , test = True )
8590
86- #rearth = 6.37122e6 #km
87- #UNITS_CONV = 86400.0*365.0*1e-9 # kg/s to Tg/yr
88- REARTH = 6.37122e6 # km
89- UNITS_CONV = 86400.0 * 365.0 * 1e-9 # kg/s to Tg/yr
91+ # rearth = 6.37122e6 #km
92+ # UNITS_CONV = 86400.0*365.0*1e-9 # kg/s to Tg/yr
93+ REARTH = 6.37122e6 # km
94+ UNITS_CONV = 86400.0 * 365.0 * 1e-9 # kg/s to Tg/yr
9095# TODO:
9196# Convert so4 unit to TgS
92- #mwso4 = 115.0
93- #mws = 32.066
94- #UNITS_CONV_S = UNITS_CONV/mwso4*mws # kg/s to TgS/yr
97+ # mwso4 = 115.0
98+ # mws = 32.066
99+ # UNITS_CONV_S = UNITS_CONV/mwso4*mws # kg/s to TgS/yr
95100
96101
97- species = ["bc" , "dst" , "mom" , "ncl" ,"pom" ,"so4" ,"soa" ]
98- SPECIES_NAMES = {"bc" : "Black Carbon" ,
102+ species = ["bc" , "dst" , "mom" , "ncl" , "pom" , "so4" , "soa" ]
103+ SPECIES_NAMES = {
104+ "bc" : "Black Carbon" ,
99105 "dst" : "Dust" ,
100106 "mom" : "Marine Organic Matter" ,
101107 "ncl" : "Sea Salt" ,
102108 "pom" : "Primary Organic Matter" ,
103109 "so4" : "Sulfate" ,
104- "soa" : "Secondary Organic Aerosol" }
110+ "soa" : "Secondary Organic Aerosol" ,
111+ }
105112MISSING_VALUE = 999.999
106113metrics_dict = {}
107114metrics_dict_ref = {}
108115
109116seasons = ["ANN" ]
110117
111118ref_data_path = os .path .join (
112- e3sm_diags .INSTALL_PATH ,
113- "control_runs" ,
114- "aerosol_global_metrics_benchmarks.json" ,
115- )
119+ e3sm_diags .INSTALL_PATH ,
120+ "control_runs" ,
121+ "aerosol_global_metrics_benchmarks.json" ,
122+ )
116123
117- with open (ref_data_path , 'r' ) as myfile :
118- ref_file = myfile .read ()
124+ with open (ref_data_path , "r" ) as myfile :
125+ ref_file = myfile .read ()
119126
120127metrics_ref = json .loads (ref_file )
121128
122129for season in seasons :
123130 for aerosol in species :
124- print (f' Aerosol species: { aerosol } ' )
131+ print (f" Aerosol species: { aerosol } " )
125132 metrics_dict [aerosol ] = generate_metrics_dic (test_data , aerosol , season )
126133 metrics_dict_ref [aerosol ] = metrics_ref [aerosol ]
127- #metrics_dict_ref[aerosol] = {
134+ # metrics_dict_ref[aerosol] = {
128135 # "Surface Emission (Tg/yr)": f'{MISSING_VALUE:.3f}',
129136 # "Sink (Tg/yr)": f'{MISSING_VALUE:.3f}',
130137 # "Dry Deposition (Tg/yr)": f'{MISSING_VALUE:.3f}',
131138 # "Wet Deposition (Tg/yr)": f'{MISSING_VALUE:.3f}',
132139 # "Burden (Tg)": f'{MISSING_VALUE:.3f}',
133140 # "Lifetime (Days)": f'{MISSING_VALUE:.3f}',
134141 # }
135-
136- with open (f' aerosol_table_{ season } .csv' , "w" ) as table_csv :
142+
143+ with open (f" aerosol_table_{ season } .csv" , "w" ) as table_csv :
137144 writer = csv .writer (
138145 table_csv ,
139146 delimiter = "," ,
140147 quotechar = "'" ,
141148 quoting = csv .QUOTE_MINIMAL ,
142- lineterminator = ' \n ' ,
149+ lineterminator = " \n " ,
143150 )
144- #writer.writerow([" ", "test","ref",])
151+ # writer.writerow([" ", "test","ref",])
145152 for key , values in metrics_dict .items ():
146153 writer .writerow ([SPECIES_NAMES [key ]])
147- print (' key' , key , values )
154+ print (" key" , key , values )
148155 for value in values :
149156 print (value )
150157 line = []
151158 line .append (value )
152159 line .append (values [value ])
153160 line .append (metrics_dict_ref [key ][value ])
154- print (line , ' line' )
161+ print (line , " line" )
155162 writer .writerows ([line ])
156163 writer .writerows (["" ])
157-
158-
159-
160-
0 commit comments