| Expand |
|---|
| title | Click here to see preparation instructions ... |
|---|
|
| Excerpt Include |
|---|
| Using FLEXPART with Metview |
|---|
| Using FLEXPART with Metview |
|---|
|
|
...
The macro starts with defining the release like this:
rel_volcano = flexpart_release(
name : "VOLCANO",
starting_date : 0,
starting_time : 15,
ending_date : 2,
ending_time : 12,
area : [63.63,-19.6,63.63,-19.6],
top_level : 9000,
bottom_level : 1651,
particle_count : 10000,
masses : 1000000
)
|
This says that the release will happen over a 45 h period between heights 1651 and 10000 m at the location of the volcano and we will release 1000 tons of material in total.
...
The actual simulation is carried out by calling flexpart_run():
r = flexpart_run(
output_path : "result_fwd_conc",
input_path : "../data",
starting_date : 20120517,
starting_time : 12,
ending_date : 20120519,
ending_time : 12,
output_field_type: "concentration",
output_flux : "on",
output_trajectory : "on",
output_area : [40,-25,66,10],
output_grid : [0.25,0.25],
output_levels : [500,1000,2000,3000,4000,5000,7500,10000,15000],
release_species : 8,
receptors : "on",
receptor_names : ["rec1","rec2"],
receptor_latitudes : [60,56.9],
receptor_longitudes : [6.43,-3.5],
releases : rel_volcano
)
print(r)
|
Here we defined both the input and output path and specified the simulation period, the output grid and levels as well. We also told FLEXPART to generate gridded concentration fields and plume trajectories on output..
...
First, we define the level (8000 m) and the parameter ("mdc") we want to plot. Then we call mvl_flexpart_read_hl() to filter the data into a fieldset.
dIn="result_fwd_conc/"
inFile=dIn & "conc_s001.grib"
lev=8000
par="mdc"
#Read fields on the given height level
g=mvl_flexpart_read_hl(inFile,par,lev,-1,-1)
|
Next, we define the contouring definition. The units we used here are ng m**-3 because for parameter "mdc" the native units (kg m**-3) are automatically scaled by the plotting library (see details about the this scaling for various FLEXPART GRIB fields here.
cont_list=[1,10,50,100,150,200,250,500,750,1000,2000,6000]
conc_shade = mcont(
legend : "on",
contour : "off",
contour_level_selection_type : "level_list",
contour_level_list : cont_list,
contour_label : "off",
contour_shade : "on",
contour_shade_method : "area_fill",
contour_shade_max_level_colour : "red",
contour_shade_min_level_colour : "RGB(0.14,0.37,0.86)",
contour_shade_colour_direction : "clockwise",
contour_method: "linear"
)
|
Next, we build the title with mvl_flexpart_title(). Please note that we need to explicitly specify the plotting units!
title=mvl_flexpart_title(g,0.3,"ng m**-3")
|
Finally we define the view with the map:
Defining the map view Expand source
and generate the plot:
plot(view,g,conc_shade,title)
|
Having run the we get a plot like this (after navigating to step 39h):
...