PrepEmisSources/example2.py at main · saneku/PrepEmisSources · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
from src import *
import numpy as np

# Example 2:
# Emisson profiles of ash and SO2 [Mt/sec] are obtained by inversion (Ukhov et al. 2023).
# Emissions are for two scenarios: with and without radiative feedback.
# Note, masses of ash and SO2 are not the same in both scenarios.
# Water vapor emission is set to 100 Mt. Umbrella profile is used for water vapor emissions.
# Inverted profiles are interpolated to 10 minute intervals.
# Inverted ash and SO2 profiles [Mt/sec] are interpolated into the vertical grid from the
# provided wrfinput file.

# Ukhov, A., Stenchikov, G., Schnell, J., Ahmadov, R., Rizza, U., Grell, G., and Hoteit, I.:
# Enhancing volcanic eruption simulations with the WRF-Chem v4.8, Geosci. Model Dev., 18, 9805–9825,
# https://doi.org/10.5194/gmd-18-9805-2025, 2025.

if __name__ == "__main__":
    # Location of the Pinatubo volcano
    LAT, LON = 15.1429, 120.3496
    YEAR, MONTH, DAY = 1991, 6, 15
    DURATION = 14*3600     # seconds

    netcdf_handler = WRFNetCDFWriter(source_dir="./")
    y,x = netcdf_handler.findClosestGridCell(LAT,LON)
    staggerred_h=netcdf_handler.getColumn_H(x,y)

    #Water vapor Emissions
    watervapor_profiles = [(VerticalProfile_Umbrella, [staggerred_h, YEAR, MONTH, DAY, 1.666, DURATION, 17000, 1000, 0.75,1])]
    #Define the water vapor scenario
    watervapor_scenario = EmissionScenario(Emission_WaterVapor(mass_mt=100.0,lat=LAT, lon=LON))
    #Add the water vapor profile to the water vapor scenario
    for p, args in watervapor_profiles:
        watervapor_scenario.add_profile(p(*args))

    scenarios_with_disabled_rad_feedback = [
                EmissionScenario_Inverted_Pinatubo(Emission_Ash(mass_mt=62.67,lat=LAT, lon=LON,bin_n=10,mean_r=2.4,stddev=1.8),
                                                    './scenarios/Pinatubo_Ukhov_2023/ash_2d_emission_profiles_rad_off'),

                EmissionScenario_Inverted_Pinatubo(Emission_SO2(mass_mt=16.73,lat=LAT, lon=LON),
                                                     './scenarios/Pinatubo_Ukhov_2023/so2_2d_emission_profiles_rad_off'),
                watervapor_scenario
                #EmissionScenario_Inverted_Pinatubo(Emission_WaterVapor(mass_mt=100.0,lat=LAT, lon=LON),
                #                                     './scenarios/Pinatubo_Ukhov_2023/so2_2d_emission_profiles_rad_off')
                #Emission_Sulfate(mass_mt=0.1,lat=15,lon=165),
                ]

    scenarios_with_enabled_rad_feedback = [
                EmissionScenario_Inverted_Pinatubo(Emission_Ash(mass_mt=66.53,lat=LAT, lon=LON,bin_n=10,mean_r=2.4,stddev=1.8),
                                                    './scenarios/Pinatubo_Ukhov_2023/ash_2d_emission_profiles_rad_on'),

                EmissionScenario_Inverted_Pinatubo(Emission_SO2(mass_mt=15.54,lat=LAT, lon=LON),
                                                     './scenarios/Pinatubo_Ukhov_2023/so2_2d_emission_profiles_rad_on'),

                watervapor_scenario
                #EmissionScenario_Inverted_Pinatubo(Emission_WaterVapor(mass_mt=0.01,lat=LAT, lon=LON),
                #                                     './scenarios/Pinatubo_Ukhov_2023/so2_2d_emission_profiles_rad_on')
                #Emission_Sulfate(mass_mt=0.1,lat=15,lon=165),
                ]

    emission_writer = EmissionWriter_NonUniformInTimeHeightProfiles(scenarios_with_enabled_rad_feedback, netcdf_handler, output_interval_m=10)

    emission_writer.plot_scenarios()
    # Plot the scenarios
    #scenarios[0].plot()
    #scenarios[1].plot()
    #scenarios[2].plot()
    #scenarios[0].plot(linestyle='--', color='grey', marker='')


    emission_writer.write()

    #for p in scenarios[0].profiles:
    #        p.plot()
    #        print (f"Mass {sum(p.values):.2f}")

    emission_writer.plot_scenarios()
    # Plot the scenarios
    #scenarios[0].plot()
    #scenarios[1].plot()
    #scenarios[2].plot()