#Metview Macro

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2019 ECMWF. This software is distributed under the terms
#  of the Apache License version 2.0. In applying this license, ECMWF does not
#  waive the privileges and immunities granted to it by virtue of its status as
#  an Intergovernmental Organization or submit itself to any jurisdiction.
# 
#  ***************************** LICENSE END ************************************
# 

# read grib - contains fields on 850 hPa
f = read("advection_850.grib")

# extract fields  
q = read(data: f, param: "q")
u = read(data: f, param: "u")
v = read(data: f, param: "v")
z = read(data: f, param: "z")

# compute specific humidity gradient 
grad = gradient(q[1])

# compute advection
adv = u[1]*grad[1] + v[1]*grad[2]

# scale results for better visualisation
adv = adv * 1E6

# define contour shading for advection
adv_cont = mcont(
    legend                      : "on",
    contour                     : "off",
    contour_level_count         : 21,
    contour_label               : "off",
    contour_max_level           : 1,
    contour_min_level           : -1,
    contour_shade               : "on",
    contour_shade_colour_method : "palette",
    contour_shade_method        : "area_fill",
    contour_shade_palette_name  : "eccharts_black_red_21"
    )

# define contouring for geopotential
z_cont = mcont(
    contour_line_colour          : "black",
    contour_highlight_colour     : "black",
    contour_level_selection_type : "interval",
    contour_interval             : 5
    )
    
# define map view
view = geoview(
    map_area_definition : "corners",
    area                : [32,-30,62,10]
    )
    
# define title    
title = mtext(text_lines: ["Geopotential 850 hPa","Specific humidity advection 850 hPa"])    
 
# define the output plot file
setoutput(pdf_output(output_name : 'advection'))   
    
# generate plot    
plot(view, adv, adv_cont, z[1], z_cont, title)

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2019 ECMWF. This software is distributed under the terms
#  of the Apache License version 2.0. In applying this license, ECMWF does not
#  waive the privileges and immunities granted to it by virtue of its status as
#  an Intergovernmental Organization or submit itself to any jurisdiction.
# 
#  ***************************** LICENSE END ************************************
# 

import metview as mv

# read grib - contains fields on 850 hPa
f = mv.read("advection_850.grib")

# extract fields  
q = mv.read(data = f, param = "q")
u = mv.read(data = f, param = "u")
v = mv.read(data = f, param = "v")
z = mv.read(data = f, param = "z")

# compute specific humidity gradient 
grad = mv.gradient(q[0])

# compute advection
adv = u[0]*grad[0] + v[0]*grad[1]

# scale results for better visualisation
adv = adv * 1E6

# define contour shading for advection
adv_cont = mv.mcont(
    legend                      = "on",
    contour                     = "off",
    contour_level_count         = 21,
    contour_label               = "off",
    contour_max_level           = 1,
    contour_min_level           = -1,
    contour_shade               = "on",
    contour_shade_colour_method = "palette",
    contour_shade_method        = "area_fill",
    contour_shade_palette_name  = "eccharts_black_red_21"
    )

# define contouring for geopotential
z_cont = mv.mcont(
    contour_line_colour          = "black",
    contour_highlight_colour     = "black",
    contour_level_selection_type = "interval",
    contour_interval             = 5
    )
    
# define map view
view = mv.geoview(
    map_area_definition = "corners",
    area                = [32,-30,62,10]
    )
    
# define title    
title = mv.mtext(text_lines = ["Geopotential 850 hPa","Specific humidity advection 850 hPa"])    
  
# define the output plot file
mv.setoutput(mv.pdf_output(output_name = 'advection'))  
    
# generate plot    
mv.plot(view, adv, adv_cont, z[0], z_cont, title)