|
1 | 1 | """ |
2 | | -Prints how effective (in terms of the depleted moles / ug of dry air) was: |
3 | | - oxidation in total (looking at the depletes S_IV) |
| 2 | +Prints how effective (in terms of the total depleted moles in the domain) was: |
| 3 | + oxidation in total (looking at the depleted S_IV and gained S_VI) |
4 | 4 | oxidation by H2O2 (looking at the depleted H2O2) |
5 | 5 | oxidation by O3 (looking at the depleted O3) |
| 6 | +Prints the initial and final dry particulate mass (total value and in % of the final mass) |
6 | 7 | """ |
7 | 8 | import numpy as np |
8 | 9 | import h5py as h5 |
9 | 10 | import sys |
| 11 | +import math |
10 | 12 |
|
11 | 13 | # libcloudph++ bindings to python (to have access to library constants) |
12 | 14 | sys.path.insert(0, "../../../../../build/bindings/python/") |
13 | 15 | from libcloudphxx import common as cm |
14 | 16 |
|
15 | | -#for case in ('case_base', 'case_base_rk', 'case3', 'case4', 'case4_no_O3', 'case5'): |
16 | | -for case in ('case4', 'case5', 'case6'): |
17 | | - print " " |
18 | | - print case |
19 | | - |
| 17 | +# path to build directory with data |
| 18 | +dir_path = '../../../build/tests/chem_sandbox/' |
| 19 | + |
| 20 | +for case in ['case_base', 'case3', 'case4', 'case5', 'case6']: |
20 | 21 |
|
21 | 22 | # open hdf5 files with data |
22 | | - h5f_ini = h5.File('data/' + case + '/out_hall_pinsky_stratocumulus/timestep0000000000.h5', 'r') |
23 | | - h5f_end = h5.File('data/' + case + '/out_hall_pinsky_stratocumulus/timestep0000011800.h5', 'r') |
| 23 | + h5f_ini = h5.File(dir_path + 'out_' + case + '/timestep0000010000.h5', 'r') # model state afer spinup |
| 24 | + h5f_end = h5.File(dir_path + 'out_' + case + '/timestep0000011800.h5', 'r') # model state at the end of simulation |
| 25 | + h5f_cst = h5.File(dir_path + 'out_' + case + '/const.h5', 'r') # constant dry air density profile |
24 | 26 |
|
| 27 | + # dry air density |
| 28 | + rho_d = h5f_cst["G"][:] # Apparently it's important to multiply |
| 29 | + # volume of grid cells # by both in order to properly compare the initial |
| 30 | + dv = np.ones(shape=(76,76)) # and final mass in the domain. |
| 31 | + dv[0,:] = 0.5 # (The sedimentation process displaces super droplets. |
| 32 | + dv[-1,:] = 0.5 # The size distribution moments reported from the libcloud |
| 33 | + dv[:,0] = 0.5 # are divided by grid-cell volume and density. |
| 34 | + dv[:,-1] = 0.5 # Without taking into account dv and rhod one gets errors ~ 1-10% of the |
| 35 | + dv[0,0] = 0.25 # total created S_VI) |
| 36 | + dv[-1,-1] = 0.25 |
| 37 | + dv[0,-1] = 0.25 |
| 38 | + dv[-1,0] = 0.25 |
| 39 | + dv *= 20. |
| 40 | + # total domain volume |
| 41 | + total_vol = dv.sum() |
| 42 | + |
25 | 43 | # helper dict for chem names and molar mass |
26 | | - #name gas molar mass aqueous molar mass label in hdf5 ini spn end |
| 44 | + #name gas molar mass aqueous molar mass label in hdf5 spn end |
27 | 45 | help_dict = { |
28 | | - 'H2O2' : [cm.M_H2O2, cm.M_H2O2, 'H2O2', 0, 0, 0], |
29 | | - 'O3' : [cm.M_O3, cm.M_O3, 'O3', 0, 0, 0], |
30 | | - 'SO2' : [cm.M_SO2, cm.M_SO2_H2O, 'S_IV', 0, 0, 0], |
31 | | - 'H2SO4': [0, cm.M_H2SO4, 'S_VI', 0, 0, 0] |
| 46 | + 'H2O2' : [cm.M_H2O2, cm.M_H2O2, 'H2O2', 0, 0], |
| 47 | + 'O3' : [cm.M_O3, cm.M_O3, 'O3', 0, 0], |
| 48 | + 'SO2' : [cm.M_SO2, cm.M_SO2_H2O, 'S_IV', 0, 0], |
| 49 | + 'H2SO4': [0, cm.M_H2SO4, 'S_VI', 0, 0] |
32 | 50 | } |
33 | 51 |
|
34 | | - # calulate the initial and final number of moles per micro gram of dry air |
35 | | - # and store them in dict |
| 52 | + # calulate the initial and final number of moles and store them in dict |
36 | 53 | # (additionally one could also check state after spinup) |
37 | 54 | for key, val in help_dict.iteritems(): |
38 | 55 |
|
39 | | - name1 = key + "g" # gas phase chem species |
40 | | - size_all = 76. * 76 |
| 56 | + name1 = key + "g" # gas phase chem species |
41 | 57 |
|
42 | 58 | if key == 'H2SO4': |
43 | 59 | name2 = 'chem_S_VI_aq' |
44 | 60 |
|
45 | | - # total concentration [moles/ug of dry air] |
46 | | - ini = (h5f_ini[name2][:] / val[1]).sum() / size_all * 1e9 |
47 | | - end = (h5f_end[name2][:] / val[1]).sum() / size_all * 1e9 |
| 61 | + # total moles |
| 62 | + ini = (h5f_ini[name2][:] * rho_d * dv).sum() / val[1] |
| 63 | + end = (h5f_end[name2][:] * rho_d * dv).sum() / val[1] |
48 | 64 |
|
49 | 65 | else: |
50 | | - name2 = "chem_" + val[2] + "_aq" # aq phase chem species |
| 66 | + name2 = "chem_" + val[2] + "_aq" # aq phase chem species |
51 | 67 |
|
52 | | - # total concentration [moles/ug of dry air] |
53 | | - ini = (h5f_ini[name1][:].sum() / val[0] + h5f_ini[name2][:].sum() / val[1]) / size_all * 1e9 |
54 | | - print " " |
55 | | - print "ini ", val[2], ini |
56 | | - end = (h5f_end[name1][:].sum() / val[0] + h5f_end[name2][:].sum() / val[1]) / size_all * 1e9 |
57 | | - print "end", val[2], end |
58 | | - print " " |
59 | | - |
60 | | - # calculate average over all |
61 | | - #ini = (total_ini_conc[:]).sum() / size_all * 1e9 |
62 | | - #end = (total_end_conc[:]).sum() / size_all * 1e9 |
| 68 | + # total moles |
| 69 | + ini = (h5f_ini[name1][:] * dv * rho_d).sum() / val[0] + (h5f_ini[name2][:] * dv * rho_d).sum() / val[1] |
| 70 | + end = (h5f_end[name1][:] * dv * rho_d).sum() / val[0] + (h5f_end[name2][:] * dv * rho_d).sum() / val[1] |
63 | 71 |
|
64 | 72 | val[3] = ini |
65 | | - val[5] = end |
66 | | - |
67 | | - print "-------------- % difference % --------------------------------" |
68 | | - |
69 | | - print "iniO3 - endO3 / init O3 ", (help_dict['O3'][3] - help_dict['O3'][5]) / help_dict['O3'][3] * 100 |
70 | | - print "iniH2O2 - endH2O2 / init H2O2 ", (help_dict['H2O2'][3] - help_dict['H2O2'][5]) / help_dict['H2O2'][3] * 100 |
71 | | - print "iniS4 - endS4 / init S4 ", (help_dict['SO2'][3] - help_dict['SO2'][5]) / help_dict['SO2'][3] * 100 |
| 73 | + val[4] = end |
72 | 74 |
|
73 | | - #print "dO3 / dSO2 ", (help_dict['O3'][3] - help_dict['O3'][5]) / (help_dict['SO2'][3] - help_dict['SO2'][5]) |
74 | | - #print "dH2O2 / dSO2 ", (help_dict['H2O2'][3] - help_dict['H2O2'][5]) / (help_dict['SO2'][3] - help_dict['SO2'][5]) |
| 75 | + # changes in moles due to oxidation |
| 76 | + dn_S6 = (help_dict['H2SO4'][4] - help_dict['H2SO4'][3]) |
| 77 | + dn_S4 = (help_dict['SO2'][3] - help_dict['SO2'][4]) |
| 78 | + dn_O3 = (help_dict['O3'][3] - help_dict['O3'][4]) |
| 79 | + dn_H2O2 = (help_dict['H2O2'][3] - help_dict['H2O2'][4]) |
| 80 | + |
| 81 | + # change in dry particulate matter |
| 82 | + ini_dry_mass = 4./3 * math.pi * (h5f_ini["rd_rng000_mom3"][:] * dv * rho_d).sum() * 1800 / total_vol * 1e9 # ug/m3 dry air |
| 83 | + fin_dry_mass = 4./3 * math.pi * (h5f_end["rd_rng000_mom3"][:] * dv * rho_d).sum() * 1800 / total_vol * 1e9 # ug/m3 dry air |
| 84 | + ini_dry_mass_chem = help_dict['H2SO4'][3] * (cm.M_NH3 + cm.M_H2SO4) / total_vol * 1e9 # ug/m3 dry air |
75 | 85 |
|
76 | 86 | print " " |
77 | | - print "ini S6", help_dict['H2SO4'][3] * val[1] , "ug/kg of dry air" |
78 | | - print "fin S6", help_dict['H2SO4'][5] * val[1] , "ug/kg of dry air" |
79 | | - print "delta S6", help_dict['H2SO4'][5] * val[1] - help_dict['H2SO4'][3] * val[1], "ug/kg of dry air" |
80 | | - print "finS6 - iniS6 / final S6 ", (help_dict['H2SO4'][5] - help_dict['H2SO4'][3]) / help_dict['H2SO4'][5] * 100 |
81 | | - |
82 | | - |
| 87 | + print "---------------------- " + case + " ------------------------------" |
| 88 | + print "init. SO2 is depleted by ", (help_dict['SO2'][3] - help_dict['SO2'][4]) / help_dict['SO2'][3] * 100, "%" |
| 89 | + print "oxidation by O3 ", dn_O3 / dn_S6 * 100, "%" |
| 90 | + print "oxidation by H2O2 ", dn_H2O2 / dn_S6 * 100, "%" |
| 91 | + print " " |
| 92 | + print "ini particulate mass (r dry) ", ini_dry_mass, "ug/m3 dry air" |
| 93 | + print "ini particulate mass (chem) ", ini_dry_mass_chem, "ug/m3 dry air" |
| 94 | + print "fin particulate mass (r_dry) ", fin_dry_mass, "ug/m3 dry air" |
| 95 | + print "delta particulate mass ", fin_dry_mass - ini_dry_mass, "ug/m3 dry air" |
| 96 | + print " " |
| 97 | + print "init mass [% final mass] ", ini_dry_mass / fin_dry_mass * 100, "%" |
| 98 | + print "created mass [% final mass] ", (fin_dry_mass - ini_dry_mass) / fin_dry_mass * 100, "%" |
0 commit comments