#Metview Macro

# **************************** 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 ************************************

# Program parameters

home      = getenv("HOME")
path      = home & "/metview/macro_tutorial/macro_tut1/"
file_name = path & "TUV_Data"
par       = "t"  #['u','v']
vf_date   = 2012-03-01
n_of_days = 5
the_area  = [30,-25,50,65] # S, W, N, E



acoast = mcoast(
        map_coastline_resolution        : "low",
        map_coastline_thickness         : 3,
        map_coastline_land_shade        : "on",
        map_coastline_land_shade_colour : "grey",
        map_coastline_sea_shade         : "on",
        map_coastline_sea_shade_colour  : "RGB(0.9,0.95,1)",
        map_grid_longitude_increment    : 10
        )

ps_atlantic = geoview(
        map_projection      : "polar_stereographic",
        map_area_definition : "corners",
        area                : the_area,
        coastlines          : acoast
        )

mollweide_view = geoview(
        map_projection : "mollweide",
        coastlines     : acoast
        )

#view = ps_atlantic       # choose the projection/view
view = mollweide_view   # choose the projection/view

# Define the display window

page = plot_page(
        view  : view
        )

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


# Check which visual definition to use: contour or wind arrows?

visdef = select_visdef (par)



# Derive the difference field

fa_diff = fc_an_diff(file_name, par, vf_date, n_of_days)


# Plot the result

plot (dw, fa_diff, visdef)



# Print some statistics on the data

print ('minvalue: ', minvalue(fa_diff[1]))
print ('maxvalue: ', maxvalue(fa_diff[1]))
print ('average:  ', average (fa_diff[1]))



# ===================================================================

# A function to compute a difference field
# Author, date, restrictions, argument types, etc,...

function fc_an_diff (fname, par, vf_date, n_of_days)

    infile = read(fname)

    fc = read(
        type    : "fc",
        param   : par,
        date    : vf_date - n_of_days,
        step    : n_of_days*24,
        data    : infile
        )

    an = read(
        type  : "an",
        param : par,
        date  : vf_date,
        step  : 0,
        data  : infile
        )

    return fc-an

end fc_an_diff


# ----------------------------------------------------------------------------
# Function      : select_visdef
#
# Author (date) : Iain Russell (21/02/2006)
#
# Description   : Selects a visual definition appropriate to the
#                 meteorological parameter supplied.
#
# Parameters    : The name, or list of names of parameters to plot
#
# Return value  : A visual definition
# ----------------------------------------------------------------------------

function select_visdef (par)


    # define our visual definitions

    coloured_wind = mwind(
            legend                                : "on",
            wind_arrow_legend_text                : "m/s",
            wind_advanced_method                  : "on",
            wind_advanced_colour_max_level_colour : "red",
            wind_advanced_colour_min_level_colour : "blue",
            wind_advanced_colour_direction        : "clockwise",
            wind_thinning_factor                  : 1
            )
    
    pos = mcont(
            contour_line_thickness       : 2,
            contour_line_colour          : "red",
            contour_highlight            : "off",
            contour_level_selection_type : "level_list",
            contour_max_level            : 10,
            contour_min_level            : 0.5,
            contour_level_list           : [0.5,1,2,4,10]
            )
    
    neg = mcont(
            contour_line_thickness       : 2,
            contour_highlight            : "off",
            contour_level_selection_type : "level_list",
            contour_max_level            : -0.5,
            contour_min_level            : -10,
            contour_level_list           : [-10,-4,-2,-1,-0.5]
            )
    

    # select the appropriate one

    if (par = ['u', 'v']) then
        selected_visdef = coloured_wind
    else if par = 't' then
        selected_visdef = [pos, neg]
    else
        print ('Could not find visual definition for parameter ', par)
        selected_visdef = ()
    end if


    # return the selected one
    
    return selected_visdef

end select_visdef



