# **************************** LICENSE START ***********************************
#
# Copyright 2012 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 ************************************

# we include the file that contains the 'function 'gradientb'
# in order to make that available to us

include "gradientb"


# set the area we wish to retrieve data from
#          N,   W,  S,  E

area_xx = [70, -45, 10, 85]


# Retrieve the specific humidity
q = retrieve(
		date	:	-1,
		param	:	"q",
		level	:	700,
		grid	:	[1.5,1.5]
		)

# Get the u and v components of the wind
u = retrieve(
		date	:	-1,
		param	:	"u",
		level	:	700,
		area	:	area_xx,
		grid	:	[1.5,1.5]
		)
v = retrieve(
		date	:	-1,
		param	:	"v",
		level	:	700,
		area	:	area_xx,
		grid	:	[1.5,1.5]
		)


# Compute the gradient of Q
q = gradientb(q)

# Extract the area we are calculating on
q = read ( area : area_xx, data : q)



# Compute the advection of Q
a = q[1]*u + q[2]*v
a = -a * (10 ^ 8) # units will be 10e-8 (kg/kg)/sec


# Plot positive advection in blue, negative in red
contour_common = (
		contour_level_selection_type : "interval",
		contour_interval             : 3,
		contour_label                : "on",
		contour_label_height         : 0.25,	
		contour_highlight            : "off",
		contour_hilo                 : "on",
		contour_hilo_type            : "number",
		contour_hilo_format          : "F5.1",
		contour_hilo_height          : 0.3
		)

cont_n = mcont(
		contour_common,
		contour_max_level       : -0.0001,
		contour_line_colour     : "red",
		contour_label_colour    : "red",
		contour_lo_colour       : "red"
		)

cont_p = mcont(
		contour_common,
		contour_min_level       : 0.0001,
		contour_line_colour     : "blue",
		contour_label_colour    : "blue",
		contour_hi_colour       : "blue"
		)

# A plot window
acoast = mcoast(
		map_coastline_resolution        : "low",
		map_grid_longitude_increment    : 10,
		map_coastline_land_shade        : "on",
		map_coastline_land_shade_colour : "cream"
		)

ps_atlantic = geoview(
		map_projection  : "polar_stereographic",
		area            : [30,-25,50,65],
		coastlines      : acoast
		)

page = plot_page(
		view            : ps_atlantic
		)

dw = plot_superpage(
		custom_width    : 29.7,
		custom_height   : 21,
		pages           : page
		)



# Now plot the result

plot(dw,a,cont_p,cont_n)

