Simulating a Point Source
In this tutorial, the steps to simulate a volcanic eruption via a point source are given.
Contents
Setting up ICON-ART
The first thing to do is having a working installation of ICON-ART. To check if your icon version has been built correctly, you can check if
your_icon_folder/bin/icon.exe
exists.
Setting up Directories
Now a directory structure has to be set up. Usually the following directories are used:
- Working Directory: Here the files that are needed for an individual run are saved. This usually includes the runscript and the relevant .xml files.
- Icon Code directory: This is where the Icon code is stored.
- External Data directory: Here external files which are needed for a run are stored. For example, to parametrize the optical properties of clouds a files like ECHAM6_CldOptProps.nc is used. These files of course can be switched out for others, however in most applications the same ones are used. A list of the used files is given in the runscript, which then creates a link of these files in the output directory.
- Output directory: This is where the new simulation data will be stored. Since most of the time large amounts of data are produced, this is stored in the work or scratch partitions on most HPC systems. The namelists produced by the runfile are also stored here.
Setting up the Runscript
The function of the runscript is to set up the directory structure, link the relevant files and create the namelists, which are grouped into several namelist files.
Setting up the .xml's
Here the .xml
pntSRC.xml
is set up. It contains all the necessary information to describe the emission of the here defined aerosols into the atmosphere. The following information is contained:
- where is the Pointsource
- when is it emitting
- what substances are emitted
- how much of each substance is emitted
- what size are the emitted substances (median and standard deviation)
This can be adapted as needed.
In this example the Raikoke eruption on 21 of June 2019 is simulated, as can be seen in the following .xml:
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!DOCTYPE tracers SYSTEM "sources_selTrnsp.dtd">
3
4 <sources>
5 <pntSrc id="Raikoke-SO2">
6 <lon type="real">153.24</lon>
7 <lat type="real">49.29</lat>
8 <substance type="char">TRSO2</substance>
9 <source_strength type="real">46300.0</source_strength>
10 <height type="real">14000</height>
11 <height_bot type="real">600</height_bot>
12 <unit type="char">kg s-1</unit>
13 <startTime type="char">2019-06-21T13:00:00</startTime>
14 <endTime type="char">2019-06-21T22:00:00</endTime>
15 </pntSrc>
16 <pntSrc id="Raikoke-ashacc">
17 <lon type="real">153.24</lon>
18 <lat type="real">49.29</lat>
19 <substance type="char">ash_insol_acc</substance>
20 <dg3_emiss type="real">0.8E-6</dg3_emiss>
21 <sigma_emiss type="real">1.4</sigma_emiss>
22 <source_strength type="real">19700.0</source_strength>
23 <height type="real">14000</height>
24 <height_bot type="real">600</height_bot>
25 <unit type="char">kg s-1</unit>
26 <startTime type="char">2019-06-21T13:00:00</startTime>
27 <endTime type="char">2019-06-21T22:00:00</endTime>
28 </pntSrc>
29 <pntSrc id="Raikoke-ashcoa">
30 <lon type="real">153.24</lon>
31 <lat type="real">49.29</lat>
32 <substance type="char">ash_insol_coa</substance>
33 <dg3_emiss type="real">2.98E-6</dg3_emiss>
34 <sigma_emiss type="real">1.4</sigma_emiss>
35 <source_strength type="real">19700.0</source_strength>
36 <height type="real">14000</height>
37 <height_bot type="real">600</height_bot>
38 <unit type="char">kg s-1</unit>
39 <startTime type="char">2019-06-21T13:00:00</startTime>
40 <endTime type="char">2019-06-21T22:00:00</endTime>
41 </pntSrc>
42 <pntSrc id="Raikoke-ashgiant">
43 <lon type="real">153.24</lon>
44 <lat type="real">49.29</lat>
45 <substance type="char">ash_giant</substance>
46 <dg3_emiss type="real">11.35E-6</dg3_emiss>
47 <sigma_emiss type="real">1.4</sigma_emiss>
48 <source_strength type="real">19700.0</source_strength>
49 <height type="real">14000</height>
50 <height_bot type="real">600</height_bot>
51 <unit type="char">kg s-1</unit>
52 <startTime type="char">2019-06-21T13:00:00</startTime>
53 <endTime type="char">2019-06-21T22:00:00</endTime>
54 </pntSrc>
55 </sources>
Running the Simulation
The setup of the model for the simulation happens with a runscript. It contains all the different control parameters, such as time control of the simulation, what equations are going to be used, whether ART is active, and a lot more. The runscript itself is just a text file, which is run via bash to start the simulation. Usually a template file is used and adapted to the specific needs of the user running the simulation. Here we will go through the runscript for our pointsource simulation.
The first step is to set the relevant directories, so that the XML data we provide as well as the path to the ICON model can be found and accessed easily. In this simulation, external data such as a grid containing external parameter data is required. The directory containing the data is set as DATADIR. The output directory is where the model will write the output NetCDF files into, as well as linking a lot of relevant files used in the simulation.
1 # working directory
2 export XMLDIR=/your/path/to/Training2023/Volcano
3
4 # datadir2 --> do not change
5 export DATADIR=/your/path/to/DATA_Volcano
6
7 # output directory
8 export OUTDIR=/your/path/to/Volcano_externally
9
10 # Code directory
11 export ICONDIR=/your/path/to/icon-kit
12 export ARTDIR=${ICONDIR}/externals/art
In the next step, the previously defined output directory is created, if not already there, and the icon executable is linked into it for easier access and overview.
1 if [ ! -d $OUTDIR ]; then
2 mkdir -p $OUTDIR
3 fi
4
5 cd $OUTDIR
6
7 # copy icon binary to OUTDIR
8 cp ${ICONDIR}/bin/icon icon.exe
Now, the external parameters and chemistry files are linked into the output directory, so that ICON and the user can find it in the output directory
1 ln -sf ${DATADIR}/icon_grid_0024_R02B06_G.nc ${OUTDIR}/iconR2B06-DOM01.nc
2 ln -sf ${DATADIR}/icon_grid_0024_R02B06_G-grfinfo.nc ${OUTDIR}/iconR2B06-DOM01-grfinfo.nc
3 ln -sf ${DATADIR}/icon_extpar_0024_R02B06_G_20200917_tiles.nc ${OUTDIR}/extpar_iconR2B06-DOM01.nc
4 ln -sf ${DATADIR}/icon_init_0024_R02B06_2019062112.nc ${OUTDIR}/igfff00000000
5
6 ln -sf ${ICONDIR}/data/rrtmg_lw.nc ${OUTDIR}/rrtmg_lw.nc
7 ln -sf ${ICONDIR}/data/ECHAM6_CldOptProps.nc ${OUTDIR}/ECHAM6_CldOptProps.nc
8
9 # Dictionary for the mapping: DWD GRIB2 names <-> ICON internal names
10 ln -sf ${ICONDIR}/run/ana_varnames_map_file.txt ${OUTDIR}/map_file.ana
11
12 ## chemistry
13 ln -sf ${ARTDIR}/runctrl_examples/init_ctrl/Simnoy2002.dat ${OUTDIR}/Simnoy2002.dat
14 ln -sf ${ARTDIR}/runctrl_examples/init_ctrl/Linoz2004Br.dat ${OUTDIR}/Linoz2004Br.dat
15
16 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_scat-aer.dat ${OUTDIR}/FJX_scat-aer.dat
17 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_j2j.dat ${OUTDIR}/FJX_j2j.dat
18 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_j2j_extended.dat ${OUTDIR}/FJX_j2j_extended.dat
19 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_scat-cld.dat ${OUTDIR}/FJX_scat-cld.dat
20 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_scat-ssa.dat ${OUTDIR}/FJX_scat-ssa.dat
21 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_scat-UMa.dat ${OUTDIR}/FJX_scat-UMa.dat
22 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_spec_extended.dat ${OUTDIR}/FJX_spec_extended.dat
23 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_spec_extended_lyman.dat ${OUTDIR}/FJX_spec_extended_lyman.dat
24 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/FJX_spec.dat ${OUTDIR}/FJX_spec.dat
25 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/atmos_std.dat ${OUTDIR}/atmos_std.dat
26 ln -sf ${ARTDIR}/runctrl_examples/photo_ctrl/atmos_h2och4.dat ${OUTDIR}/atmos_h2och4.dat
In the next step, namelist files are created. They contain the individual namelists which control the settings of the different modules of icon. They are written into the output directory and read by icon during a run. First the master namelist is created, which most importantly contains the start and stop date, as well as the possibility to restart a run. The next namelist file created is called "NAMELIST_Raikoke-June-2019". Here a lot of the model setup is taking place. It is subdivided into different namelists, which all do their part to set up a different part of the model. Some important namelists for this simulation are:
- grid_nml (line 35): Here the grid on which the simulation is taking place is given
- run_nml (line 49): Contains the switch to turn on ICON-ART, as well as set model height levels and timestep length.
- output_nml (line 225): Defines which variables at which intervals are written in the output directory
- art_nml (line 238): Contains the switches relevant for ICON-ART, see Namelist for more details. Most importantly, the xml file for the point source emission is linked here
1 cd $OUTDIR
2 cat > icon_master.namelist << EOF
3 &master_nml
4 lRestart = .false.
5 /
6 &master_time_control_nml
7 experimentStartDate = "2019-06-21T12:00:00"
8 forecastLeadTime = "P1D"
9 checkpointTimeIntval = "P10D"
10 /
11 &master_model_nml
12 model_type=1
13 model_name="ATMO"
14 model_namelist_filename="NAMELIST_Raikoke-June-2019"
15 model_min_rank=1
16 model_max_rank=65536
17 model_inc_rank=1
18 /
19 &time_nml
20 ini_datetime_string = "2019-06-21T12:00:00"
21 dt_restart = 864000 ! 10 days
22 /
23 EOF
24
25 cat > NAMELIST_Raikoke-June-2019 << EOF
26 ¶llel_nml
27 nproma = 8 ! optimal setting 8 for CRAY; use 16 or 24 for IBM
28 p_test_run = .false.
29 l_test_openmp = .false.
30 l_log_checks = .false.
31 num_io_procs = 1 ! up to one PE per output stream is possible
32 itype_comm = 1
33 iorder_sendrecv = 3 ! best value for CRAY (slightly faster than option 1)
34 /
35 &grid_nml
36 dynamics_grid_filename = 'iconR2B06-DOM01.nc'
37 dynamics_parent_grid_id = 0
38 lredgrid_phys = .false.
39 /
40 &initicon_nml
41 init_mode = 7 ! 2: IFS; 7: Initialized DWD Analysis
42 lread_ana = .FALSE. ! no analysis data will be read
43 dwdfg_filename = "igfff00000000" ! initial data filename
44 filetype = 4
45 ana_varnames_map_file = "map_file.ana" ! dictionary mapping internal names onto GRIB2 shortNames
46 ltile_coldstart = .TRUE. ! coldstart for surface tiles
47 ltile_init = .FALSE. ! set it to .TRUE. if FG data originate from run without tiles
48 /
49 &run_nml
50 num_lev = 90
51 lvert_nest = .true. ! use vertical nesting if a nest is active
52 !nsteps = ${nsteps} ! 50 ! 1200 ! 7200 !
53 dtime = 180 ! timestep in seconds
54 ldynamics = .TRUE. ! dynamics
55 ltransport = .true.
56 iforcing = 3 ! NWP forcing
57 ltestcase = .false. ! false: run with real data
58 msg_level = 7 ! print maximum wind speeds every 5 time steps
59 ltimer = .true. ! set .TRUE. for timer output
60 timers_level = 15 ! can be increased up to 10 for detailed timer output
61 output = "nml"
62 lart = .true.
63 /
64 &nwp_phy_nml
65 inwp_gscp = 1
66 inwp_convection = 1
67 inwp_radiation = 4
68 inwp_cldcover = 1
69 inwp_turb = 1
70 inwp_satad = 1
71 inwp_sso = 1
72 inwp_gwd = 1
73 inwp_surface = 1
74 icapdcycl = 3 ! apply CAPE modification to improve diurnalcycle over tropical land (optimizes NWP scores)
75 latm_above_top = .false., .true. ! the second entry refers to the nested domain (if present)
76 efdt_min_raylfric = 3600 ! 7200.
77 itype_z0 = 2
78 icpl_aero_conv = 0 ! check meaning -> default 0 - off
79 icpl_aero_gscp = 0 ! check meaning -> default 0 - off
80 ! resolution-dependent settings - please choose the appropriate one
81 dt_rad = 2160.
82 dt_conv = 720.
83 dt_sso = 1440.
84 dt_gwd = 1440.
85 mu_rain = 0.5
86 rain_n0_factor = 0.1
87 lshallowconv_only = .true.
88 lgrayzone_deepconv = .false.
89 /
90 &nwp_tuning_nml
91 itune_albedo = 1 ! somewhat reduced albedo (w.r.t. MODIS data) over Sahara in order to reduce cold bias
92 tune_gkdrag = 0.09
93 tune_gkwake = 1.8
94 tune_gfrcrit = 0.333
95 tune_dust_abs = 1.
96 tune_zvz0i = 1.1
97 tune_box_liq_asy = 4.0
98 tune_gust_factor = 7.0
99 tune_rcucov = 0.075
100 tune_rhebc_land = 0.825
101 tune_zvz0i = 0.85
102 icpl_turb_clc = 2
103 max_calibfac_clcl = 2.0
104 tune_box_liq = 0.04
105 /
106 &turbdiff_nml
107 tkhmin = 0.75 ! new default since rev. 16527
108 tkmmin = 0.75 ! "
109 pat_len = 750. ! 750 in exp.nh_oper
110 c_diff = 0.2
111 rat_sea = 0.8 ! ** new value since for v2.0.15; previously 8.0 ** ! ** new since r20191: 8.5 for R3B7, 8.0 for R2B6 in order to get similar temperature biases in the tropics **
112 ltkesso = .true.
113 frcsmot = 0.2 ! these 2 switches together apply vertical smoothing of the TKE source terms
114 imode_frcsmot = 2 ! in the tropics (only), which reduces the moist bias in the tropical lower troposphere:
115 itype_sher = 3
116 ltkeshs = .true.
117 a_hshr = 2.0
118 alpha1 = 0.125
119 icldm_turb = 1 ! found in exp.nh_oper ** new recommendation for v2.0.15 in conjunction with evaporation fix for grid-scale rain **
120 rlam_heat = 10.0
121 /
122 &lnd_nml
123 ntiles = 3 !!! 1 for assimilation cycle and forecast
124 nlev_snow = 3 !!! 1 for assimilation cycle and forecast
125 lmulti_snow = .false. !!! .false. for assimilation cycle and forecast
126 itype_heatcond = 3
127 idiag_snowfrac = 20
128 lsnowtile = .true. !! later on .true. if GRIB encoding issues are solved
129 lseaice = .true.
130 llake = .true.
131 lprog_albsi = .true.
132 itype_lndtbl = 4 ! minimizes moist/cold bias in lower tropical troposphere
133 itype_evsl = 4
134 itype_root = 2
135 itype_trvg = 3
136 cwimax_ml = 5.e-4
137 c_soil = 1.25
138 c_soil_urb = 0.5
139 sstice_mode = 2
140 /
141 &radiation_nml
142 irad_o3 = 7
143 irad_aero = 6
144 albedo_type = 2 ! Modis albedo
145 vmr_co2 = 390.e-06 ! values representative for 2012
146 vmr_ch4 = 1800.e-09
147 vmr_n2o = 322.0e-09
148 vmr_o2 = 0.20946
149 vmr_cfc11 = 240.e-12
150 vmr_cfc12 = 532.e-12
151 direct_albedo_water = 3
152 albedo_whitecap = 1
153 ecrad_data_path = "${ICONDIR}/externals/ecrad/data"
154 /
155 &nonhydrostatic_nml
156 iadv_rhotheta = 2
157 ivctype = 2
158 itime_scheme = 4
159 exner_expol = 0.333
160 vwind_offctr = 0.2
161 damp_height = 44000.
162 rayleigh_coeff = 1
163 lhdiff_rcf = .true.
164 divdamp_order = 24 ! for data assimilation runs, '2' provides extra-strong filtering of gravity waves
165 divdamp_type = 32 !!! optional: 2 for assimilation cycle if very strong gravity-wave filtering is desired
166 divdamp_fac = 0.004
167 l_open_ubc = .false.
168 igradp_method = 3
169 l_zdiffu_t = .true.
170 thslp_zdiffu = 0.02
171 thhgtd_zdiffu = 125.
172 htop_moist_proc= 22500.
173 hbot_qvsubstep = 22500. ! use 19000. with R3B7, must be at least as large as htop_moist_proc
174 htop_aero_proc = 25000. !height limit for aerosol tracer processes
175 /
176 &sleve_nml
177 min_lay_thckn = 20.
178 max_lay_thckn = 400. ! maximum layer thickness below htop_thcknlimit
179 htop_thcknlimit = 14000. ! this implies that the upcoming COSMO-EU nest will have 60 levels
180 top_height = 75000.
181 stretch_fac = 0.9
182 decay_scale_1 = 4000.
183 decay_scale_2 = 2500.
184 decay_exp = 1.2
185 flat_height = 16000.
186 /
187 &dynamics_nml
188 iequations = 3
189 idiv_method = 1
190 divavg_cntrwgt = 0.50
191 lcoriolis = .TRUE.
192 /
193 &transport_nml
194 ivadv_tracer = 3,3,3,3,3,3,3,3,3,3,3
195 itype_hlimit = 3,4,4,4,4,3,3,3,3,3,3
196 ihadv_tracer = 52,2,2,2,2,22,22,22,22,22,22
197 iadv_tke = 0
198 /
199 &diffusion_nml
200 hdiff_order = 5
201 itype_vn_diffu = 1
202 itype_t_diffu = 2
203 hdiff_efdt_ratio = 24.0
204 hdiff_smag_fac = 0.025
205 lhdiff_vn = .TRUE.
206 lhdiff_temp = .TRUE.
207 /
208 &interpol_nml
209 nudge_zone_width = 8
210 lsq_high_ord = 3
211 l_intp_c2l = .true.
212 l_mono_c2l = .true.
213 /
214 &extpar_nml
215 itopo = 1
216 n_iter_smooth_topo = 1
217 heightdiff_threshold = 3000.
218 hgtdiff_max_smooth_topo = 750., ! found in exp.nh_oper ** should be changed to 750.,750 with next Extpar update! **
219 itype_lwemiss = 2
220 /
221 &io_nml
222 itype_pres_msl = 5 ! 5: new DWD method; 4: IFS method with bug fix for self-consistency between SLP and geopotential
223 itype_rh = 1 ! RH w.r.t. water
224 /
225 &output_nml
226 filetype = 4 ! output format: 2=GRIB2, 4=NETCDFv2
227 dom = -1 ! write all domains
228 mode = 1 ! 1 = forecast mode with TAXIS_RELATIVE, works only with output in multiples of 1h; 2 = climate mode, default, TAXIS_ABSOLUTE
229 output_bounds = 0., 691200., 600. ! start, end, increment
230 steps_per_file = 1
231 include_last = .TRUE.
232 output_filename = 'Raikoke-June-2019-forecast_mode-remap' ! file name base
233 ml_varlist = 'TRSO2_chemtr','TRH2SO4_chemtr','group:ART_AEROSOL','rho','pres','temp', 'z_mc','z_ifc'
234 remap = 1
235 reg_lon_def = 120.,0.1,280.
236 reg_lat_def = 85.,-0.1, 30.
237 /
238 &art_nml
239 lart_diag_out = .TRUE.
240 lart_aerosol = .TRUE.
241 lart_pntSrc = .TRUE.
242 lart_chem = .TRUE.
243 lart_chemtracer = .TRUE.
244 iart_ari = 0
245
246 cart_aerosol_xml = '${XMLDIR}/Ex01_aerotracer.xml'
247 cart_modes_xml = '${XMLDIR}/Ex01_modes.xml'
248 cart_pntSrc_xml = '${XMLDIR}/Ex01_pntSrc.xml'
249 cart_chemtracer_xml = '${XMLDIR}/Ex01_chemtracer.xml'
250 /
251 EOF
With the namelists set, we can now prepare the simulation for running. Since a lot of setup is required for running, this is also prepared within the runscript and then later executed. The following part is highly dependant on the HPC system configuration. The following setup works on the DKRZ Levante HPC system (July 2023):
1 cat > ${OUTDIR}/job_ICON << ENDFILE
2 #!/bin/bash
3 #SBATCH --job-name=EXP01_ART # Specify job name
4 #SBATCH --partition=compute # Specify partition name
5 #SBATCH --nodes=4 # Specify number of nodes
6 #SBATCH --ntasks-per-node=128 # Specify number of (MPI) tasks on each node
7 #SBATCH --time=01:00:00 # Set a limit on the total run time
8 #SBATCH --mail-type=FAIL # Notify user by email in case of job failure
9 #SBATCH --mail-user=example@mail.com # Set your e-mail address
10 #SBATCH --account=your_project_account # Charge resources on this project account
11
12
13 module load gcc
14 module load intel-oneapi-compilers/2022.0.1-gcc-11.2.0
15 module load intel-oneapi-mkl/2022.0.1-gcc-11.2.0
16 module load openmpi/4.1.2-intel-2021.5.0
17 module load netcdf-c/4.8.1-openmpi-4.1.2-intel-2021.5.0
18 module load netcdf-fortran/4.5.3-openmpi-4.1.2-intel-2021.5.0
19 module load parallel-netcdf/1.12.2-openmpi-4.1.2-intel-2021.5.0
20 module load hdf5/1.12.1-openmpi-4.1.2-intel-2021.5.0
21 module load eccodes/2.21.0-intel-2021.5.0
22
23 unset SLURM_EXPORT_ENV
24 unset SLURM_MEM_PER_NODE
25 unset SBATCH_EXPORT
26
27 export LD_LIBRARY_PATH="/usr/lib:/usr/lib64:/sw/spack-levante/netcdf-c-4.8.1-2k3cmu/lib:/sw/spack-levante/netcdf-fortran-4.5.3-k6xq5g/lib:/sw/spack-levante/hdf5-1.12.1-tvymb5/lib:/sw/spack-levante/eccodes-2.21.0-3ehkbb/lib64:/sw/spack-levante/intel-oneapi-mkl-2022.0.1-ttdktf/mkl/2022.0.1/lib/intel64/"
28
29 export OMPI_MCA_pml="ucx"
30 export OMPI_MCA_btl=self
31 export OMPI_MCA_osc="pt2pt"
32 export UCX_IB_ADDR_TYPE=ib_global
33
34 export OMPI_MCA_coll="^ml,hcoll"
35 export OMPI_MCA_coll_hcoll_enable="0"
36 export HCOLL_ENABLE_MCAST_ALL="0"
37 export HCOLL_MAIN_IB=mlx5_0:1
38 export UCX_NET_DEVICES=mlx5_0:1
39 export UCX_TLS=mm,knem,cma,dc_mlx5,dc_x,self
40 export UCX_UNIFIED_MODE=y
41 export HDF5_USE_FILE_LOCKING=FALSE
42 export OMPI_MCA_io="romio321"
43 export UCX_HANDLE_ERRORS=bt
44
45
46 ulimit -s 102400
47 ulimit -c 0
48
49 srun -l --cpu_bind=cores --distribution=block:cyclic --propagate=STACK,CORE ./icon.exe
50
51
52 ENDFILE
This concludes the runfile. The Job_ICON-file that has been produced in the last step can then be executed. For that, the runfile hast to be run using bash:
$ chmod +x myrunfile
$ ./myrunfile
This will then create all the Namelist files and the runscript job_ICON. To submit the icon run, go to the output directory, make the runscript executable, and run it:
$ cd /path/to/your/output
$ chmod +x job_ICON
$ sbatch job_ICON
Now the Job should be submitted. The submission status can be seen with
$ squeue --user your_username
If everything went well, in the output directory there will be output files generated with names like Raikoke-June-2019-forecast_mode-remap_DOM01_ML_0001.nc
Inspecting the Output
In this example, the output is created in the form of several NetCDF4 (.nc) files. To avoid having one giant file, it is subdivided into several files containing a fixed amount of timesteps. The quickest way to inspect a NetCDF file is with CDO:
$ cdo sinfo Raikoke-June-2019-forecast_mode-remap_DOM01_ML_0001.nc
This way it is possible to obtain information about the variables, grid parameters and timesteps of the simulation. To actually see the values of a specific variable, tools like ncdump can be used, however in many cases this is not practical, since the amount of single values is way too large to look at individual data points.
Making a Plot of the Simulation Data
In order to gain insight into the simulation, as well as an overview, the most common way is to create some plots of different variables of interest. In our case, we want to know what happens to the ash emitted from the Raikoke eruption. A quick and easy way to create plots is using python and a reference script to produce a plot similar to the Pointsource animation is given here:
1 # first we import the following packages:
2 import numpy as np # general utility for trigonomic functions and arrays
3 import matplotlib.pyplot as plt # Library to produce Plots
4 import xarray as xr # reading in of netCDF data, optimized for large datasets
5 from glob import glob
6 import cartopy.crs as ccrs # plotting things on a map
7
8 #Next step: read in the Netcdf Files
9
10 #datadir_externally = 'path/to/your/output/Training2023/Volcano_externally/'
11 files_externally = glob(datadir_externally+'Raikoke-June-2019-forecast_mode-remap_DOM01_ML*')
12 files_externally.sort()
13
14 #Next step: Calculate relevant data to be displayed
15 # get model level height
16 def get_dz(ds):
17 dz = -1 * ds.z_ifc.diff('height_2')
18 dz = dz.rename({'height_2':'height'})
19 dz = dz.assign_coords(height=(dz.height - 1))
20 return dz
21
22 # Calculate column integrated tracer mass (tracer load)
23 def tracer_load(tr,rho,dz):
24 m_tr = tr * rho * dz
25 m_tr = m_tr.sum('height')
26 m_tr = m_tr.squeeze()
27 return m_tr
28
29 #
30 MM_SO2 = 64/28.96 # molar weight of SO2 / dry air
31
32
33 # set the timestep you want to plot
34 # in this simulation there are 144 timesteps and 144 files each containing one timestep
35 # To produce a plot for each timestep, you can also loop the rest of this file over i
36 i = 30
37
38 #load the NetCDF file containing the timestep
39 ds = xr.open_dataset(files_externally[i],autoclose=True)
40
41 #calculate Heightlevels
42 dz = get_dz(ds)
43
44 # Read in SO2 Column
45 SO2 = tracer_load(ds.TRSO2_chemtr.values*MM_SO2,ds.rho,dz)
46
47 # set colorbar levels
48 color_so2 = np.linspace(1e-5, 0.02, 11,endpoint=True)
49
50
51 # Plot SO2 column
52 # set the size of the figure and the background color to white
53 fig = plt.figure(figsize=(8,6),facecolor='w')
54
55 # Set the axes using the specified map projection
56 ax=plt.axes(projection=ccrs.NearsidePerspective(central_longitude=180, central_latitude=49.29,satellite_height=35785831/5, globe=None))
57 ax.stock_img() #add a stock image as background
58 ax.coastlines() # add coastlines
59 ax.set_title(str(ds.time.values)[2:12]+" "+str(ds.time.values)[13:18]) #set title to the current time displayed
60 ax.gridlines(draw_labels=False, dms=False, x_inline=True, y_inline=True) #add lat lon gridlines
61
62 # Make a filled contour plot
63 plot1=ax.contourf(ds.lon, ds.lat, SO2, levels=color_so2, cmap='plasma',
64 transform = ccrs.PlateCarree())
65 kwargs = {'format': '%.3f'}
66 cbar = plt.colorbar(plot1,location='right',ticks=color_so2,**kwargs)
67 cbar.set_label('SO2 Column Mass (kg m-2)')
68
69 #save the plot at a path of your choice
70
71 #plt.savefig("path/to/your/Plots/folder/SO2_{:03d}.png".format(i))
72 plt.close(fig)
To run the script, you need a working python version and the packages listed in the beginning of the script. Dont forget to make the file executable in order to run it using chmod +x plotscript .
The resulting Plot will look something like this: 
