# -*- coding: iso-8859-15 -*-
#import pylab as plt
# From http://stackoverflow.com/questions/13336823/matplotlib-python-error
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
from matplotlib.pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import os
from netCDF4 import Dataset as NetCDFFile
import nc_var_tools as ncvar
import generic_tools as gen

errormsg = 'ERROR -- error -- ERROR -- error'
warnmsg = 'WARNING -- waring -- WARNING -- warning'

fillValue = 1.e20
fillValueF = 1.e20

colorsauto = ['#FF0000', '#00FF00', '#0000FF', '#FF00FF', '#00FFFF', '#FFAA00',     \
  '#AA0000', '#00AA00', '#0000AA', '#AA00AA', '#00AAAA', '#AA3200']
pointkindsauto = ['.', ',', 'x', '+', '*', '|', '_', 'o', '<', '>', 'v', '^',       \
  's', 'D', 'p' ,'h' ,'H']
linekindsauto = ['-', '--', '-.', ':']
linewidthsauto = [1.]
pointsizesauto = [7.]

####### Funtions
# searchInlist:
# datetimeStr_datetime:
# dateStr_date:
# numVector_String:
# timeref_datetime:
# slice_variable: Function to return a slice of a given variable according to values to its dimension
# interpolate_locs:
# datetimeStr_conversion:
# percendone:
# netCDFdatetime_realdatetime:
# file_nlines:
# variables_values:
# check_colorBar:
# units_lunits:
# ASCII_LaTeX:
# pretty_int:
# DegGradSec_deg:
# intT2dt:
# lonlat2D: Function to return lon, lat 2D matrices from any lon,lat matrix
# lonlat_values:
# date_CFtime:
# pot_values:
# CFtimes_plot:
# color_lines:
# output_kind:
# check_arguments:
# Str_Bool:
# plot_points:
# plot_2Dfield:
# plot_2Dfield_easy:
# plot_topo_geogrid:
# plot_topo_geogrid_boxes:
# plot_2D_shadow:
# plot_2D_shadow_time: Plotting a 2D field with one of the axes being time
# plot_Neighbourghood_evol:Plotting neighbourghood evolution# plot_Trajectories
# plot_2D_shadow_contour:
# plot_2D_shadow_contour_time:
# dxdy_lonlat: Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
# plot_2D_shadow_line:
# plot_lines: Function to plot a collection of lines
# plot_ZQradii: Function to plot following radial averages only at exact grid poins
# transform: Function to transform the values and the axes

# From nc_var_tools.py
def reduce_spaces(string):
    """ Function to give words of a line of text removing any extra space
    """
    values = string.replace('\n','').split(' ')
    vals = []
    for val in values:
         if len(val) > 0:
             vals.append(val)

    return vals

def searchInlist(listname, nameFind):
    """ Function to search a value within a list
    listname = list
    nameFind = value to find
    >>> searInlist(['1', '2', '3', '5'], '5')
    True
    """
    for x in listname:
      if x == nameFind:
        return True
        break
    return False

def datetimeStr_datetime(StringDT):
    """ Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
    >>> datetimeStr_datetime('1976-02-17_00:00:00')
    1976-02-17 00:00:00
    """
    import datetime as dt

    fname = 'datetimeStr_datetime'

    dateD = np.zeros((3), dtype=int)
    timeT = np.zeros((3), dtype=int)

#    quit()

    dateD[0] = int(StringDT[0:4])
    dateD[1] = int(StringDT[5:7])
    dateD[2] = int(StringDT[8:10])

    trefT = StringDT.find(':')
    if not trefT == -1:
#        print '  ' + fname + ': refdate with time!'
        timeT[0] = int(StringDT[11:13])
        timeT[1] = int(StringDT[14:16])
        timeT[2] = int(StringDT[17:19])

    if int(dateD[0]) == 0:
        print warnmsg
        print '    ' + fname + ': 0 reference year!! changing to 1'
        dateD[0] = 1 
 
    newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])

    return newdatetime

def dateStr_date(StringDate):
  """ Function to transform a string date ([YYYY]-[MM]-[DD] format) to a date object
  >>> dateStr_date('1976-02-17')
  1976-02-17
  """
  import datetime as dt

  dateD = StringDate.split('-')
  if int(dateD[0]) == 0:
    print warnmsg
    print '    dateStr_date: 0 reference year!! changing to 1'
    dateD[0] = 1
  newdate = dt.date(int(dateD[0]), int(dateD[1]), int(dateD[2]))
  return newdate

def numVector_String(vec,char):
    """ Function to transform a vector of numbers to a single string [char] separated
    numVector_String(vec,char)
      vec= vector with the numerical values
      char= single character to split the values
    >>> print numVector_String(np.arange(10),' ')
    0 1 2 3 4 5 6 7 8 9
    """
    fname = 'numVector_String'

    if vec == 'h':
        print fname + '_____________________________________________________________'
        print numVector_String.__doc__
        quit()

    Nvals = len(vec)

    string=''
    for i in range(Nvals):
        if i == 0:
            string = str(vec[i])
        else:
            string = string + char + str(vec[i])

    return string

def timeref_datetime(refd, timeval, tu):
    """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to datetime object
    refd: time of reference (as datetime object)
    timeval: time value (as [tu] from [tref])
    tu: time units
    >>> timeref = date(1949,12,1,0,0,0)
    >>> timeref_datetime(timeref, 229784.36, hours)
    1976-02-17 08:21:36
    """
    import datetime as dt
    import numpy as np

## Not in timedelta
#    if tu == 'years':
#        realdate = refdate + dt.timedelta(years=float(timeval))
#    elif tu == 'months':
#        realdate = refdate + dt.timedelta(months=float(timeval))
    if tu == 'weeks':
        realdate = refd + dt.timedelta(weeks=float(timeval))
    elif tu == 'days':
        realdate = refd + dt.timedelta(days=float(timeval))
    elif tu == 'hours':
        realdate = refd + dt.timedelta(hours=float(timeval))
    elif tu == 'minutes':
        realdate = refd + dt.timedelta(minutes=float(timeval))
    elif tu == 'seconds':
        realdate = refd + dt.timedelta(seconds=float(timeval))
    elif tu == 'milliseconds':
        realdate = refd + dt.timedelta(milliseconds=float(timeval))
    else:
          print errormsg
          print '    timeref_datetime: time units "' + tu + '" not ready!!!!'
          quit(-1)

    return realdate

def slice_variable(varobj, dimslice):
    """ Function to return a slice of a given variable according to values to its 
      dimensions
    slice_variable(varobj, dims)
      varobj= object wit the variable
      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension 
        [value]: 
          * [integer]: which value of the dimension
          * -1: all along the dimension
          * [beg]:[end] slice from [beg] to [end]
          * -9: last value of the dimension

    """
    fname = 'slice_variable'

    if varobj == 'h':
        print fname + '_____________________________________________________________'
        print slice_variable.__doc__
        quit()

    vardims = varobj.dimensions
    Ndimvar = len(vardims)

    Ndimcut = len(dimslice.split('|'))
    if Ndimcut == 0:
        Ndimcut = 1
        dimcut = list(dimslice)

    dimsl = dimslice.split('|')

    varvalsdim = []
    dimnslice = []
    monodim = []
    for idd in range(Ndimvar):
        found = False
        for idc in range(Ndimcut):
            dimcutn = dimsl[idc].split(':')[0]
            dimcutv = dimsl[idc].split(':')[1]
            if vardims[idd] == dimcutn:
                posfrac = dimcutv.find('@')
                if posfrac != -1:
                    inifrac = int(dimcutv.split('@')[0])
                    endfrac = int(dimcutv.split('@')[1])
                    varvalsdim.append(slice(inifrac,endfrac))
                    dimnslice.append(vardims[idd])
                else:
                    if int(dimcutv) == -1:
                        varvalsdim.append(slice(0,varobj.shape[idd]))
                        dimnslice.append(vardims[idd])
                    elif int(dimcutv) == -9:
                        varvalsdim.append(varobj.shape[idd]-1)
                        monodim.append(vardims[idd])
                    else:
                        if int(dimcutv) > varobj.shape[idd]-1:
                            print errormsg
                            print '  ' + fname + ': given value:', int(dimcutv),     \
                              "is larger than length of dimension '" + dimcutn +     \
                              "':", varobj.shape[idd],'(C-like!) !!'
                            quit(-1)
                        varvalsdim.append(int(dimcutv))
                        monodim.append(vardims[idd])
                found = True
                break
        if not found and not searchInlist(dimnslice,vardims[idd]) and                \
          not searchInlist(monodim,vardims[idd]):
            varvalsdim.append(slice(0,varobj.shape[idd]))
            dimnslice.append(vardims[idd])

    varvalues = varobj[tuple(varvalsdim)]
    varvalues = np.squeeze(varobj[tuple(varvalsdim)])

    return varvalues, dimnslice

def get_str_nc(varo, Lchar):
    """ Function to get string values in a netCDF variable as a chains of 1char values
    varo= netCDF variable object
    Lchar = length of the string in the netCDF file
    """
    fname = 'get_str_nc'

    Nvalues = varo.shape[0]

    strings = []

    for ival in range(Nvalues):
        stringv = varo[ival,:]
        string = str('')
        for ich in range(Lchar):
            charval = str(stringv[ich:ich+1]).replace('[','').replace(']','').replace("'","")
            if charval == '--' or len(charval) == 0:
                break
            else:
                if ich == 0:
                    string = charval
                else:
                    string = string + charval

        strings.append(string.strip())

    return strings

def interpolate_locs(locs,coords,kinterp):
    """ Function to provide interpolate locations on a given axis
    interpolate_locs(locs,axis,kinterp)
      locs= locations to interpolate
      coords= axis values with the reference of coordinates
      kinterp: kind of interpolation
        'lin': linear
    >>> coordinates = np.arange((10), dtype=np.float)
    >>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
    >>> interpolate_locs(values,coordinates,'lin')
    [ -1.2   2.4   5.6   7.8  13. ]
    >>> coordinates[0] = 0.5
    >>> coordinates[2] = 2.5
    >>> interpolate_locs(values,coordinates,'lin')
    [ -3.4          1.93333333   5.6          7.8         13.        ]
    """

    fname = 'interpolate_locs'

    if locs == 'h':
        print fname + '_____________________________________________________________'
        print interpolate_locs.__doc__
        quit()

    Nlocs = locs.shape[0]
    Ncoords = coords.shape[0]

    dcoords = coords[Ncoords-1] - coords[0]

    intlocs = np.zeros((Nlocs), dtype=np.float)
    minc = np.min(coords)
    maxc = np.max(coords)

    for iloc in range(Nlocs):
        for icor in range(Ncoords-1):
            if locs[iloc] < minc and dcoords > 0.:
                a = 0.
                b = 1. / (coords[1] - coords[0])
                c = coords[0]
            elif locs[iloc] > maxc and dcoords > 0.:
                a = (Ncoords-1)*1.
                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
                c = coords[Ncoords-2]
            elif locs[iloc] < minc and dcoords < 0.:
                a = (Ncoords-1)*1.
                b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
                c = coords[Ncoords-2]
            elif locs[iloc] > maxc and dcoords < 0.:
                a = 0.
                b = 1. / (coords[1] - coords[0])
                c = coords[0]
            elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.:
                a = icor*1.
                b = 1. / (coords[icor+1] - coords[icor])
                c = coords[icor]
                print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b
            elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.:
                a = icor*1.
                b = 1. / (coords[icor+1] - coords[icor])
                c = coords[icor]

        if kinterp == 'lin':
            intlocs[iloc] = a + (locs[iloc] - c)*b
        else:
            print errormsg
            print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
            quit(-1)

    return intlocs

def datetimeStr_conversion(StringDT,typeSi,typeSo):
    """ Function to transform a string date to an another date object
    StringDT= string with the date and time
    typeSi= type of datetime string input
    typeSo= type of datetime string output
      [typeSi/o]
        'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
        'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
        'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
        'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
        'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
        'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
        'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H], 
          [H], ':', [M], [M], ':', [S], [S]
    >>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
    [1976    2   17    8   32    5]
    >>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
    1980/03/05 18-00-00
    """
    import datetime as dt

    fname = 'datetimeStr_conversion'

    if StringDT[0:1] == 'h':
        print fname + '_____________________________________________________________'
        print datetimeStr_conversion.__doc__
        quit()

    if typeSi == 'cfTime':
        timeval = np.float(StringDT.split(',')[0])
        tunits = StringDT.split(',')[1].split(' ')[0]
        Srefdate = StringDT.split(',')[1].split(' ')[2]

# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
##
        yrref=Srefdate[0:4]
        monref=Srefdate[5:7]
        dayref=Srefdate[8:10]

        trefT = Srefdate.find(':')
        if not trefT == -1:
#            print '  ' + fname + ': refdate with time!'
            horref=Srefdate[11:13]
            minref=Srefdate[14:16]
            secref=Srefdate[17:19]
            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
              '_' + horref + ':' + minref + ':' + secref)
        else:
            refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref +    \
              '_00:00:00')

        if tunits == 'weeks':
            newdate = refdate + dt.timedelta(weeks=float(timeval))
        elif tunits == 'days':
            newdate = refdate + dt.timedelta(days=float(timeval))
        elif tunits == 'hours':
            newdate = refdate + dt.timedelta(hours=float(timeval))
        elif tunits == 'minutes':
            newdate = refdate + dt.timedelta(minutes=float(timeval))
        elif tunits == 'seconds':
            newdate = refdate + dt.timedelta(seconds=float(timeval))
        elif tunits == 'milliseconds':
            newdate = refdate + dt.timedelta(milliseconds=float(timeval))
        else:
              print errormsg
              print '    timeref_datetime: time units "' + tunits + '" not ready!!!!'
              quit(-1)

        yr = newdate.year
        mo = newdate.month
        da = newdate.day
        ho = newdate.hour
        mi = newdate.minute
        se = newdate.second
    elif typeSi == 'matYmdHMS':
        yr = StringDT[0]
        mo = StringDT[1]
        da = StringDT[2]
        ho = StringDT[3]
        mi = StringDT[4]
        se = StringDT[5]
    elif typeSi == 'YmdHMS':
        yr = int(StringDT[0:4])
        mo = int(StringDT[4:6])
        da = int(StringDT[6:8])
        ho = int(StringDT[8:10])
        mi = int(StringDT[10:12])
        se = int(StringDT[12:14])
    elif typeSi == 'Y-m-d_H:M:S':
        dateDT = StringDT.split('_')
        dateD = dateDT[0].split('-')
        timeT = dateDT[1].split(':')
        yr = int(dateD[0])
        mo = int(dateD[1])
        da = int(dateD[2])
        ho = int(timeT[0])
        mi = int(timeT[1])
        se = int(timeT[2])
    elif typeSi == 'Y-m-d H:M:S':
        dateDT = StringDT.split(' ')
        dateD = dateDT[0].split('-')
        timeT = dateDT[1].split(':')
        yr = int(dateD[0])
        mo = int(dateD[1])
        da = int(dateD[2])
        ho = int(timeT[0])
        mi = int(timeT[1])
        se = int(timeT[2])
    elif typeSi == 'Y/m/d H-M-S':
        dateDT = StringDT.split(' ')
        dateD = dateDT[0].split('/')
        timeT = dateDT[1].split('-')
        yr = int(dateD[0])
        mo = int(dateD[1])
        da = int(dateD[2])
        ho = int(timeT[0])
        mi = int(timeT[1])
        se = int(timeT[2])
    elif typeSi == 'WRFdatetime':
        yr = int(StringDT[0])*1000 + int(StringDT[1])*100 + int(StringDT[2])*10 +    \
          int(StringDT[3])
        mo = int(StringDT[5])*10 + int(StringDT[6])
        da = int(StringDT[8])*10 + int(StringDT[9])
        ho = int(StringDT[11])*10 + int(StringDT[12])
        mi = int(StringDT[14])*10 + int(StringDT[15])
        se = int(StringDT[17])*10 + int(StringDT[18])
    else:
        print errormsg
        print '  ' + fname + ': type of String input date "' + typeSi +              \
          '" not ready !!!!'
        quit(-1)

    if typeSo == 'matYmdHMS':
        dateYmdHMS = np.zeros((6), dtype=int)
        dateYmdHMS[0] =  yr
        dateYmdHMS[1] =  mo
        dateYmdHMS[2] =  da
        dateYmdHMS[3] =  ho
        dateYmdHMS[4] =  mi
        dateYmdHMS[5] =  se
    elif typeSo == 'YmdHMS':
        dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) +        \
          str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
    elif typeSo == 'Y-m-d_H:M:S':
        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
          str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
          str(se).zfill(2)
    elif typeSo == 'Y-m-d H:M:S':
        dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' +               \
          str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
          str(se).zfill(2)
    elif typeSo == 'Y/m/d H-M-S':
        dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' +               \
          str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
          str(se).zfill(2)
    elif typeSo == 'WRFdatetime':
        dateYmdHMS = []
        yM = yr/1000
        yC = (yr-yM*1000)/100
        yD = (yr-yM*1000-yC*100)/10
        yU = yr-yM*1000-yC*100-yD*10

        mD = mo/10
        mU = mo-mD*10
        
        dD = da/10
        dU = da-dD*10

        hD = ho/10
        hU = ho-hD*10

        miD = mi/10
        miU = mi-miD*10

        sD = se/10
        sU = se-sD*10

        dateYmdHMS.append(str(yM))
        dateYmdHMS.append(str(yC))
        dateYmdHMS.append(str(yD))
        dateYmdHMS.append(str(yU))
        dateYmdHMS.append('-')
        dateYmdHMS.append(str(mD))
        dateYmdHMS.append(str(mU))
        dateYmdHMS.append('-')
        dateYmdHMS.append(str(dD))
        dateYmdHMS.append(str(dU))
        dateYmdHMS.append('_')
        dateYmdHMS.append(str(hD))
        dateYmdHMS.append(str(hU))
        dateYmdHMS.append(':')
        dateYmdHMS.append(str(miD))
        dateYmdHMS.append(str(miU))
        dateYmdHMS.append(':')
        dateYmdHMS.append(str(sD))
        dateYmdHMS.append(str(sU))
    else:
        print errormsg
        print '  ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
        quit(-1)

    return dateYmdHMS

def percendone(nvals,tot,percen,msg):
    """ Function to provide the percentage of an action across the matrix
    nvals=number of values
    tot=total number of values
    percen=percentage frequency for which the message is wanted
    msg= message
    """
    from sys import stdout

    num = int(tot * percen/100)
    if (nvals%num == 0): 
        print '\r        ' + msg + '{0:8.3g}'.format(nvals*100./tot) + ' %',
        stdout.flush()

    return ''

def netCDFdatetime_realdatetime(units, tcalendar, times):
    """ Function to transfrom from netCDF CF-compilant times to real time
    """
    import datetime as dt

    txtunits = units.split(' ')
    tunits = txtunits[0]
    Srefdate = txtunits[len(txtunits) - 1]

# Calendar type
##
    is360 = False
    if tcalendar is not None:
      print '  netCDFdatetime_realdatetime: There is a calendar attribute'
      if tcalendar == '365_day' or tcalendar == 'noleap':
          print '    netCDFdatetime_realdatetime: No leap years!'
          isleapcal = False
      elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian':
          isleapcal = True
      elif tcalendar == '360_day':
          is360 = True
          isleapcal = False
      else:
          print errormsg
          print '    netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!'
          quit(-1)

# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
##
    timeval = Srefdate.find(':')

    if not timeval == -1:
        print '  netCDFdatetime_realdatetime: refdate with time!'
        refdate = datetimeStr_datetime(Srefdate)
    else:
        refdate = dateStr_date(Srefdate + '_00:00:00')

    dimt = len(times)
#    datetype = type(dt.datetime(1972,02,01))
#    realdates = np.array(dimt, datetype)
#    print realdates

## Not in timedelta
#  if tunits == 'years':
#    for it in range(dimt):
#      realdate = refdate + dt.timedelta(years=float(times[it]))
#      realdates[it] = int(realdate.year)
#  elif tunits == 'months':
#    for it in range(dimt):
#      realdate = refdate + dt.timedelta(months=float(times[it]))
#      realdates[it] = int(realdate.year)
#    realdates = []
    realdates = np.zeros((dimt, 6), dtype=int)
    if tunits == 'weeks':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(weeks=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'days':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(days=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'hours':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(hours=float(times[it]))
#            if not isleapcal:
#                Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year))
#                realdate = realdate - dt.timedelta(days=Nleapdays)
#            if is360:
#                Nyears360 = int(realdate.year) - int(refdate.year) + 1
#                realdate = realdate -dt.timedelta(days=Nyears360*5)
#            realdates[it] = realdate
#        realdates = refdate + dt.timedelta(hours=float(times))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'minutes':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(minutes=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'seconds':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(seconds=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'milliseconds':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    elif tunits == 'microseconds':
        for it in range(dimt):
            realdate = refdate + dt.timedelta(microseconds=float(times[it]))
            realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
    else:
        print errormsg
        print '  netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!'
        quit(-1)

    return realdates

def file_nlines(filen):
    """ Function to provide the number of lines of a file
    filen= name of the file
    >>> file_nlines('trajectory.dat')
    49
    """
    fname = 'file_nlines'

    if not os.path.isfile(filen):
        print errormsg
        print '  ' + fname + ' file: "' + filen + '" does not exist !!'
        quit(-1)

    fo = open(filen,'r')

    nlines=0
    for line in fo: nlines = nlines + 1

    fo.close()

    return nlines

def realdatetime1_CFcompilant(time, Srefdate, tunits):
    """ Function to transform a matrix with a real time value ([year, month, day, 
      hour, minute, second]) to a netCDF one
        time= matrix with time
        Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
        tunits= units of time respect to Srefdate
    >>> realdatetime1_CFcompilant([1976, 2, 17, 8, 20, 0], '19491201000000', 'hours')
    229784.33333333
    """ 

    import datetime as dt
    yrref=int(Srefdate[0:4])
    monref=int(Srefdate[4:6])
    dayref=int(Srefdate[6:8])
    horref=int(Srefdate[8:10])
    minref=int(Srefdate[10:12])
    secref=int(Srefdate[12:14])
 
    refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)

    if tunits == 'weeks':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5])-refdate
        cfdates = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
    elif tunits == 'days':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
        cfdates = cfdate.days + cfdate.seconds/(3600.*24.)
    elif tunits == 'hours':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
        cfdates = cfdate.days*24. + cfdate.seconds/3600.
    elif tunits == 'minutes':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
        cfdates = cfdate.days*24.*60. + cfdate.seconds/60.
    elif tunits == 'seconds':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
        cfdates = cfdate.days*24.*3600. + cfdate.seconds
    elif tunits == 'milliseconds':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
        cfdates = cfdate.days*1000.*24.*3600. + cfdate.seconds*1000.
    elif tunits == 'microseconds':
        cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],times[5]) - refdate
        cfdates = cfdate.days*1000000.*24.*3600. + cfdate.seconds*1000000.
    else:
        print errormsg
        print '  ' + fname + ': time units "' + tunits + '" is not ready!!!'
        quit(-1)

    return cfdates

def basicvardef(varobj, vstname, vlname, vunits):
    """ Function to give the basic attributes to a variable
    varobj= netCDF variable object
    vstname= standard name of the variable
    vlname= long name of the variable
    vunits= units of the variable
    """
    attr = varobj.setncattr('standard_name', vstname)
    attr = varobj.setncattr('long_name', vlname)
    attr = varobj.setncattr('units', vunits)

    return 

def variables_values(varName):
    """ Function to provide values to plot the different variables values from ASCII file
      'variables_values.dat'
    variables_values(varName)
      [varName]= name of the variable
        return: [var name], [std name], [minimum], [maximum], 
          [long name]('|' for spaces), [units], [color palette] (following: 
          http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
     [varn]: original name of the variable
       NOTE: It might be better doing it with an external ASII file. But then we 
         got an extra dependency...
    >>> variables_values('WRFght')
    ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
    """
    import subprocess as sub

    fname='variables_values'

    if varName == 'h':
        print fname + '_____________________________________________________________'
        print variables_values.__doc__
        quit()

# This does not work....
#    folderins = sub.Popen(["pwd"], stdout=sub.PIPE)
#    folder = list(folderins.communicate())[0].replace('\n','')
# From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python
    folder = os.path.dirname(os.path.realpath(__file__))

    infile = folder + '/variables_values.dat'

    if not os.path.isfile(infile):
        print errormsg
        print '  ' + fname + ": File '" + infile + "' does not exist !!"
        quit(-1)

# Variable name might come with a statistical surname...
    stats=['min','max','mean','stdv', 'sum']

# Variables with a statistical section on their name...
    NOstatsvars = ['zmaxth', 'zmax_th', 'lmax_th', 'lmaxth']

    ifst = False
    if not searchInlist(NOstatsvars, varName.lower()):
        for st in stats:
            if varName.find(st) > -1:
                print '    '+ fname + ": varibale '" + varName + "' with a " +       \
                  "statistical surname: '",st,"' !!"
                Lst = len(st)
                LvarName = len(varName)
                varn = varName[0:LvarName - Lst]
                ifst = True
                break
    if not ifst:
        varn = varName

    ncf = open(infile, 'r')

    for line in ncf:
        if line[0:1] != '#':
            values = line.replace('\n','').split(',')
            if len(values) != 8:
                print errormsg
                print "problem in varibale: '", values[0],                           \
                  "' it should have 8 values and it has:", len(values)
                quit(-1)

            if varn[0:6] == 'varDIM': 
# Variable from a dimension (all with 'varDIM' prefix)
                Lvarn = len(varn)
                varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                 \
                  "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1',  \
                   'rainbow']
            else:
                varvals = [values[1].replace(' ',''), values[2].replace(' ',''),     \
                  np.float(values[3]), np.float(values[4]),values[5].replace(' ',''),\
                  values[6].replace(' ',''), values[7].replace(' ','')]
            if values[0] == varn:
                ncf.close()
                return varvals
                break

    print errormsg
    print '  ' + fname + ": variable '" + varn + "' not defined !!!"
    ncf.close()
    quit(-1)

    return 

def variables_values_old(varName):
    """ Function to provide values to plot the different variables
    variables_values(varName)
      [varName]= name of the variable
        return: [var name], [std name], [minimum], [maximum], 
          [long name]('|' for spaces), [units], [color palette] (following: 
          http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
     [varn]: original name of the variable
       NOTE: It might be better doing it with an external ASII file. But then we 
         got an extra dependency...
    >>> variables_values('WRFght')
    ['z', 'geopotential_height', 0.0, 80000.0, 'geopotential|height', 'm2s-2', 'rainbow']
    """
    fname='variables_values'

    if varName == 'h':
        print fname + '_____________________________________________________________'
        print variables_values.__doc__
        quit()

# Variable name might come with a statistical surname...
    stats=['min','max','mean','stdv', 'sum']

    ifst = False
    for st in stats:
        if varName.find(st) > -1:
            print '    '+ fname + ": varibale '" + varName + "' with a statistical "+\
              " surname: '",st,"' !!"
            Lst = len(st)
            LvarName = len(varName)
            varn = varName[0:LvarName - Lst]
            ifst = True
            break
    if not ifst:
        varn = varName

    if varn[0:6] == 'varDIM': 
# Variable from a dimension (all with 'varDIM' prefix)
        Lvarn = len(varn)
        varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1.,                         \
          "variable|from|size|of|dimension|'" + varn[6:Lvarn+1] + "'", '1', 'rainbox']
    elif varn == 'a_tht' or varn == 'LA_THT':
        varvals = ['ath', 'total_thermal_plume_cover', 0., 1.,                       \
        'total|column|thermal|plume|cover', '1', 'YlGnBu']
    elif varn == 'acprc' or varn == 'RAINC':
        varvals = ['acprc', 'accumulated_cmulus_precipitation', 0., 3.e4,            \
          'accumulated|cmulus|precipitation', 'mm', 'Blues']
    elif varn == 'acprnc' or varn == 'RAINNC':
        varvals = ['acprnc', 'accumulated_non-cmulus_precipitation', 0., 3.e4,       \
          'accumulated|non-cmulus|precipitation', 'mm', 'Blues']
    elif varn == 'bils' or varn == 'LBILS':
        varvals = ['bils', 'surface_total_heat_flux', -100., 100.,                   \
          'surface|total|heat|flux', 'Wm-2', 'seismic']
    elif varn == 'landcat' or varn == 'category':
        varvals = ['landcat', 'land_categories', 0., 22., 'land|categories', '1',    \
          'rainbow']
    elif varn == 'c' or varn == 'QCLOUD' or varn == 'oliq' or varn == 'OLIQ':
        varvals = ['c', 'condensed_water_mixing_ratio', 0., 3.e-4,                   \
          'condensed|water|mixing|ratio', 'kgkg-1', 'BuPu']
    elif varn == 'ci' or varn == 'iwcon' or varn == 'LIWCON':
        varvals = ['ci', 'cloud_iced_water_mixing_ratio', 0., 0.0003,                \
         'cloud|iced|water|mixing|ratio', 'kgkg-1', 'Purples']
    elif varn == 'cl' or varn == 'lwcon' or varn == 'LLWCON':
        varvals = ['cl', 'cloud_liquidwater_mixing_ratio', 0., 0.0003,               \
         'cloud|liquid|water|mixing|ratio', 'kgkg-1', 'Blues']
    elif varn == 'cld' or varn == 'CLDFRA' or varn == 'rneb' or varn == 'lrneb' or   \
      varn == 'LRNEB':
        varvals = ['cld', 'cloud_area_fraction', 0., 1., 'cloud|fraction', '1',      \
          'gist_gray']
    elif varn == 'cldc' or varn == 'rnebcon' or varn == 'lrnebcon' or                \
      varn == 'LRNEBCON':
        varvals = ['cldc', 'convective_cloud_area_fraction', 0., 1.,                 \
          'convective|cloud|fraction', '1', 'gist_gray']
    elif varn == 'cldl' or varn == 'rnebls' or varn == 'lrnebls' or varn == 'LRNEBLS':
        varvals = ['cldl', 'large_scale_cloud_area_fraction', 0., 1.,                \
          'large|scale|cloud|fraction', '1', 'gist_gray']
    elif varn == 'clt' or varn == 'CLT' or varn == 'cldt' or                         \
      varn == 'Total cloudiness':
        varvals = ['clt', 'cloud_area_fraction', 0., 1., 'total|cloud|cover', '1',   \
          'gist_gray']
    elif varn == 'cll' or varn == 'cldl' or varn == 'LCLDL' or                       \
      varn == 'Low-level cloudiness':
        varvals = ['cll', 'low_level_cloud_area_fraction', 0., 1.,                   \
          'low|level|(p|>|680|hPa)|cloud|fraction', '1', 'gist_gray']
    elif varn == 'clm' or varn == 'cldm' or varn == 'LCLDM' or                       \
      varn == 'Mid-level cloudiness':
        varvals = ['clm', 'mid_level_cloud_area_fraction', 0., 1.,                   \
          'medium|level|(440|<|p|<|680|hPa)|cloud|fraction', '1', 'gist_gray']
    elif varn == 'clh' or varn == 'cldh' or varn == 'LCLDH' or                       \
      varn == 'High-level cloudiness':
        varvals = ['clh', 'high_level_cloud_area_fraction', 0., 1.,                  \
          'high|level|(p|<|440|hPa)|cloud|fraction', '1', 'gist_gray']
    elif varn == 'clmf' or varn == 'fbase' or varn == 'LFBASE':
        varvals = ['clmf', 'cloud_base_max_flux', -0.3, 0.3, 'cloud|base|max|flux',  \
          'kgm-2s-1', 'seismic']
    elif varn == 'clp' or varn == 'pbase' or varn == 'LPBASE':
        varvals = ['clp', 'cloud_base_pressure', -0.3, 0.3, 'cloud|base|pressure',   \
          'Pa', 'Reds']
    elif varn == 'cpt' or varn == 'ptconv' or varn == 'LPTCONV':
        varvals = ['cpt', 'convective_point', 0., 1., 'convective|point', '1',       \
          'seismic']
    elif varn == 'dqajs' or varn == 'LDQAJS':
        varvals = ['dqajs', 'dry_adjustment_water_vapor_tendency', -0.0003, 0.0003,  \
        'dry|adjustment|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqcon' or varn == 'LDQCON':
        varvals = ['dqcon', 'convective_water_vapor_tendency', -3e-8, 3.e-8,         \
        'convective|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqdyn' or varn == 'LDQDYN':
        varvals = ['dqdyn', 'dynamics_water_vapor_tendency', -3.e-7, 3.e-7,          \
        'dynamics|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqeva' or varn == 'LDQEVA':
        varvals = ['dqeva', 'evaporation_water_vapor_tendency', -3.e-6, 3.e-6,       \
        'evaporation|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqlscst' or varn == 'LDQLSCST':
        varvals = ['dqlscst', 'stratocumulus_water_vapor_tendency', -3.e-7, 3.e-7,   \
        'stratocumulus|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqlscth' or varn == 'LDQLSCTH': 
        varvals = ['dqlscth', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,        \
        'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqlsc' or varn == 'LDQLSC':
        varvals = ['dqlsc', 'condensation_water_vapor_tendency', -3.e-6, 3.e-6,      \
        'condensation|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqphy' or varn == 'LDQPHY':
        varvals = ['dqphy', 'physics_water_vapor_tendency', -3.e-7, 3.e-7,           \
        'physics|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqthe' or varn == 'LDQTHE':
        varvals = ['dqthe', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7,          \
        'thermal|plumes|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqvdf' or varn == 'LDQVDF':
        varvals = ['dqvdf', 'vertical_difussion_water_vapor_tendency', -3.e-8, 3.e-8,\
        'vertical|difussion|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dqwak' or varn == 'LDQWAK':
        varvals = ['dqwak', 'wake_water_vapor_tendency', -3.e-7, 3.e-7,              \
        'wake|water|vapor|tendency', 'kg/kg/s', 'seismic']
    elif varn == 'dta' or varn == 'tnt' or varn == 'LTNT':
        varvals = ['dta', 'tendency_air_temperature', -3.e-3, 3.e-3,                 \
        'tendency|of|air|temperature', 'K/s', 'seismic']
    elif varn == 'dtac' or varn == 'tntc' or varn == 'LTNTC':
        varvals = ['dtac', 'moist_convection_tendency_air_temperature', -3.e-3,      \
        3.e-3, 'moist|convection|tendency|of|air|temperature', 'K/s', 'seismic']
    elif varn == 'dtar' or varn == 'tntr' or varn == 'LTNTR':
        varvals = ['dtar', 'radiative_heating_tendency_air_temperature', -3.e-3,     \
          3.e-3, 'radiative|heating|tendency|of|air|temperature', 'K/s', 'seismic']
    elif varn == 'dtascpbl' or varn == 'tntscpbl' or varn == 'LTNTSCPBL':
        varvals = ['dtascpbl',                                                       \
          'stratiform_cloud_precipitation_BL_mixing_tendency_air_temperature',       \
          -3.e-6, 3.e-6,                                                             \
          'stratiform|cloud|precipitation|Boundary|Layer|mixing|tendency|air|'       +
          'temperature', 'K/s', 'seismic']
    elif varn == 'dtajs' or varn == 'LDTAJS':
        varvals = ['dtajs', 'dry_adjustment_thermal_tendency', -3.e-5, 3.e-5,        \
        'dry|adjustment|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtcon' or varn == 'LDTCON':
        varvals = ['dtcon', 'convective_thermal_tendency', -3.e-5, 3.e-5,            \
        'convective|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtdyn' or varn == 'LDTDYN':
        varvals = ['dtdyn', 'dynamics_thermal_tendency', -3.e-4, 3.e-4,              \
        'dynamics|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dteva' or varn == 'LDTEVA':
        varvals = ['dteva', 'evaporation_thermal_tendency', -3.e-3, 3.e-3,           \
        'evaporation|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtlscst' or varn == 'LDTLSCST':
        varvals = ['dtlscst', 'stratocumulus_thermal_tendency', -3.e-4, 3.e-4,       \
        'stratocumulus|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtlscth' or varn == 'LDTLSCTH':
        varvals = ['dtlscth', 'thermals_thermal_tendency', -3.e-4, 3.e-4,            \
        'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtlsc' or varn == 'LDTLSC':
        varvals = ['dtlsc', 'condensation_thermal_tendency', -3.e-3, 3.e-3,          \
        'condensation|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtlwr' or varn == 'LDTLWR':
        varvals = ['dtlwr', 'long_wave_thermal_tendency', -3.e-3, 3.e-3, \
        'long|wave|radiation|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtphy' or varn == 'LDTPHY':
        varvals = ['dtphy', 'physics_thermal_tendency', -3.e-4, 3.e-4,               \
        'physics|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtsw0' or varn == 'LDTSW0':
        varvals = ['dtsw0', 'cloudy_sky_short_wave_thermal_tendency', -3.e-4, 3.e-4, \
        'cloudy|sky|short|wave|radiation|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtthe' or varn == 'LDTTHE':
        varvals = ['dtthe', 'thermals_thermal_tendency', -3.e-4, 3.e-4,              \
        'thermal|plumes|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtvdf' or varn == 'LDTVDF':
        varvals = ['dtvdf', 'vertical_difussion_thermal_tendency', -3.e-5, 3.e-5,    \
        'vertical|difussion|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'dtwak' or varn == 'LDTWAK':
        varvals = ['dtwak', 'wake_thermal_tendency', -3.e-4, 3.e-4,                  \
        'wake|thermal|tendency', 'K/s', 'seismic']
    elif varn == 'ducon' or varn == 'LDUCON':
        varvals = ['ducon', 'convective_eastward_wind_tendency', -3.e-3, 3.e-3,      \
        'convective|eastward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'dudyn' or varn == 'LDUDYN':
        varvals = ['dudyn', 'dynamics_eastward_wind_tendency', -3.e-3, 3.e-3,        \
        'dynamics|eastward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'duvdf' or varn == 'LDUVDF':
        varvals = ['duvdf', 'vertical_difussion_eastward_wind_tendency', -3.e-3,     \
         3.e-3, 'vertical|difussion|eastward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'dvcon' or varn == 'LDVCON':
        varvals = ['dvcon', 'convective_difussion_northward_wind_tendency', -3.e-3,  \
         3.e-3, 'convective|northward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'dvdyn' or varn == 'LDVDYN':
        varvals = ['dvdyn', 'dynamics_northward_wind_tendency', -3.e-3,              \
         3.e-3, 'dynamics|difussion|northward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'dvvdf' or varn == 'LDVVDF':
        varvals = ['dvvdf', 'vertical_difussion_northward_wind_tendency', -3.e-3,    \
         3.e-3, 'vertical|difussion|northward|wind|tendency', 'ms-2', 'seismic']
    elif varn == 'etau' or varn == 'ZNU':
        varvals = ['etau', 'etau', 0., 1, 'eta values on half (mass) levels', '-',   \
        'reds']
    elif varn == 'evspsbl' or varn == 'LEVAP' or varn == 'evap' or varn == 'SFCEVPde':
        varvals = ['evspsbl', 'water_evaporation_flux', 0., 1.5e-4,                  \
          'water|evaporation|flux', 'kgm-2s-1', 'Blues']
    elif varn == 'evspsbl' or varn == 'SFCEVPde':
        varvals = ['evspsblac', 'water_evaporation_flux_ac', 0., 1.5e-4,             \
          'accumulated|water|evaporation|flux', 'kgm-2', 'Blues']
    elif varn == 'g' or varn == 'QGRAUPEL':
        varvals = ['g', 'grauepl_mixing_ratio', 0., 0.0003, 'graupel|mixing|ratio',  \
          'kgkg-1', 'Purples']
    elif varn == 'h2o' or varn == 'LH2O':
        varvals = ['h2o', 'water_mass_fraction', 0., 3.e-2,                          \
          'mass|fraction|of|water', '1', 'Blues']
    elif varn == 'h' or varn == 'QHAIL':
        varvals = ['h', 'hail_mixing_ratio', 0., 0.0003, 'hail|mixing|ratio',        \
          'kgkg-1', 'Purples']
    elif varn == 'hfls' or varn == 'LH' or varn == 'LFLAT' or varn == 'flat':
        varvals = ['hfls', 'surface_upward_latent_heat_flux', -400., 400.,           \
          'upward|latnt|heat|flux|at|the|surface', 'Wm-2', 'seismic']
    elif varn == 'hfss' or varn == 'LSENS' or varn == 'sens' or varn == 'HFX':
        varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150.,         \
          'upward|sensible|heat|flux|at|the|surface', 'Wm-2', 'seismic']
    elif varn == 'hfso' or varn == 'GRDFLX':
        varvals = ['hfso', 'downward_heat_flux_in_soil', -150., 150.,                \
          'Downward|soil|heat|flux', 'Wm-2', 'seismic']
    elif varn == 'hus' or varn == 'WRFrh' or varn == 'LMDZrh' or varn == 'rhum' or   \
      varn == 'LRHUM':
        varvals = ['hus', 'specific_humidity', 0., 1., 'specific|humidty', '1',      \
          'BuPu']
    elif varn == 'huss' or varn == 'WRFrhs' or varn == 'LMDZrhs' or varn == 'rh2m' or\
      varn == 'LRH2M':
        varvals = ['huss', 'specific_humidity', 0., 1., 'specific|humidty|at|2m',    \
          '1', 'BuPu']
    elif varn == 'i' or varn == 'QICE':
        varvals = ['i', 'iced_water_mixing_ratio', 0., 0.0003,                       \
         'iced|water|mixing|ratio', 'kgkg-1', 'Purples']
    elif varn == 'lat' or varn == 'XLAT' or varn == 'XLAT_M' or varn == 'latitude':
        varvals = ['lat', 'latitude', -90., 90., 'latitude', 'degrees North',        \
          'seismic']
    elif varn == 'lcl' or varn == 's_lcl' or varn == 'ls_lcl' or varn == 'LS_LCL':
        varvals = ['lcl', 'condensation_level', 0., 2500., 'level|of|condensation',  \
          'm', 'Greens']
    elif varn == 'lambdath' or varn == 'lambda_th' or varn == 'LLAMBDA_TH':
        varvals = ['lambdath', 'thermal_plume_vertical_velocity', -30., 30.,         \
          'thermal|plume|vertical|velocity', 'm/s', 'seismic']
    elif varn == 'lmaxth' or varn == 'LLMAXTH':
        varvals = ['lmaxth', 'upper_level_thermals', 0., 100., 'upper|level|thermals'\
          , '1', 'Greens']
    elif varn == 'lon' or varn == 'XLONG' or varn == 'XLONG_M':
        varvals = ['lon', 'longitude', -180., 180., 'longitude', 'degrees East',     \
          'seismic']
    elif varn == 'longitude':
        varvals = ['lon', 'longitude', 0., 360., 'longitude', 'degrees East',        \
          'seismic']
    elif varn == 'orog' or varn == 'HGT' or varn == 'HGT_M':
        varvals = ['orog', 'orography',  0., 3000., 'surface|altitude', 'm','terrain']
    elif varn == 'pfc' or varn == 'plfc' or varn == 'LPLFC':
        varvals = ['pfc', 'pressure_free_convection', 100., 1100.,                   \
          'pressure|free|convection', 'hPa', 'BuPu']
    elif varn == 'plcl' or varn == 'LPLCL':
        varvals = ['plcl', 'pressure_lifting_condensation_level', 700., 1100.,       \
          'pressure|lifting|condensation|level', 'hPa', 'BuPu']
    elif varn == 'pr' or varn == 'RAINTOT' or varn == 'precip' or                    \
      varn == 'LPRECIP' or varn == 'Precip Totale liq+sol':
        varvals = ['pr', 'precipitation_flux', 0., 1.e-4, 'precipitation|flux',      \
          'kgm-2s-1', 'BuPu']
    elif varn == 'prprof' or varn == 'vprecip' or varn == 'LVPRECIP':
        varvals = ['prprof', 'precipitation_profile', 0., 1.e-3,                     \
          'precipitation|profile', 'kg/m2/s', 'BuPu']
    elif varn == 'prprofci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
        varvals = ['prprofci', 'precipitation_profile_convective_i', 0., 1.e-3,      \
          'precipitation|profile|convective|i', 'kg/m2/s', 'BuPu']
    elif varn == 'prprofcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
        varvals = ['prprofcl', 'precipitation_profile_convective_l', 0., 1.e-3,      \
          'precipitation|profile|convective|l', 'kg/m2/s', 'BuPu']
    elif varn == 'prprofli' or varn == 'pr_lsc_i' or varn == 'LPR_LSC_I':
        varvals = ['prprofli', 'precipitation_profile_large_scale_i', 0., 1.e-3,     \
          'precipitation|profile|large|scale|i', 'kg/m2/s', 'BuPu']
    elif varn == 'prprofll' or varn == 'pr_lsc_l' or varn == 'LPR_LSC_L':
        varvals = ['prprofll', 'precipitation_profile_large_scale_l', 0., 1.e-3,     \
          'precipitation|profile|large|scale|l', 'kg/m2/s', 'BuPu']
    elif varn == 'pracc' or varn == 'ACRAINTOT':
        varvals = ['pracc', 'precipitation_amount', 0., 100.,                        \
          'accumulated|precipitation', 'kgm-2', 'BuPu']
    elif varn == 'prc' or varn == 'LPLUC' or varn == 'pluc' or varn == 'WRFprc' or   \
      varn == 'RAINCde':
        varvals = ['prc', 'convective_precipitation_flux', 0., 2.e-4,                \
          'convective|precipitation|flux', 'kgm-2s-1', 'Blues']
    elif varn == 'prci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
        varvals = ['prci', 'convective_ice_precipitation_flux', 0., 0.003,           \
          'convective|ice|precipitation|flux', 'kgm-2s-1', 'Purples']
    elif varn == 'prcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
        varvals = ['prcl', 'convective_liquid_precipitation_flux', 0., 0.003,        \
          'convective|liquid|precipitation|flux', 'kgm-2s-1', 'Blues']
    elif varn == 'pres' or varn == 'presnivs' or varn == 'pressure' or               \
      varn == 'lpres' or varn == 'LPRES':
        varvals = ['pres', 'air_pressure', 0., 103000., 'air|pressure', 'Pa',        \
          'Blues']
    elif varn == 'prls' or varn == 'WRFprls' or varn == 'LPLUL' or varn == 'plul' or \
       varn == 'RAINNCde':
        varvals = ['prls', 'large_scale_precipitation_flux', 0., 2.e-4,              \
          'large|scale|precipitation|flux', 'kgm-2s-1', 'Blues']
    elif varn == 'prsn' or varn == 'SNOW' or varn == 'snow' or varn == 'LSNOW':
        varvals = ['prsn', 'snowfall', 0., 1.e-4, 'snowfall|flux', 'kgm-2s-1', 'BuPu']
    elif varn == 'prw' or varn == 'WRFprh':
        varvals = ['prw', 'atmosphere_water_vapor_content', 0., 10.,                 \
          'water|vapor"path', 'kgm-2', 'Blues']
    elif varn == 'ps' or varn == 'psfc' or varn =='PSFC' or varn == 'psol' or        \
      varn == 'Surface Pressure':
        varvals=['ps', 'surface_air_pressure', 85000., 105400., 'surface|pressure',  \
          'hPa', 'cool']
    elif varn == 'psl' or varn == 'mslp' or varn =='WRFmslp':
        varvals=['psl', 'air_pressure_at_sea_level', 85000., 104000.,                \
          'mean|sea|level|pressure', 'Pa', 'Greens']
    elif varn == 'qth' or varn == 'q_th' or varn == 'LQ_TH':
        varvals = ['qth', 'thermal_plume_total_water_content', 0., 25.,              \
          'total|water|cotent|in|thermal|plume', 'mm', 'YlOrRd']
    elif varn == 'r' or varn == 'QVAPOR' or varn == 'ovap' or varn == 'LOVAP':
        varvals = ['r', 'water_mixing_ratio', 0., 0.03, 'water|mixing|ratio',        \
          'kgkg-1', 'BuPu']
    elif varn == 'r2' or varn == 'Q2':
        varvals = ['r2', 'water_mixing_ratio_at_2m', 0., 0.03, 'water|mixing|' +     \
          'ratio|at|2|m','kgkg-1', 'BuPu']
    elif varn == 'rsds' or varn == 'SWdnSFC' or varn == 'SWdn at surface' or         \
      varn == 'SWDOWN':
        varvals=['rsds', 'surface_downwelling_shortwave_flux_in_air',  0., 1200.,    \
          'downward|SW|surface|radiation', 'Wm-2' ,'Reds']
    elif varn == 'rsdsacc':
        varvals=['rsdsacc', 'accumulated_surface_downwelling_shortwave_flux_in_air', \
          0., 1200., 'accumulated|downward|SW|surface|radiation', 'Wm-2' ,'Reds']
    elif varn == 'rvor' or varn == 'WRFrvor':
        varvals = ['rvor', 'air_relative_vorticity', -2.5E-3, 2.5E-3,                \
          'air|relative|vorticity', 's-1', 'seismic']
    elif varn == 'rvors' or varn == 'WRFrvors':
        varvals = ['rvors', 'surface_air_relative_vorticity', -2.5E-3, 2.5E-3,       \
          'surface|air|relative|vorticity', 's-1', 'seismic']
    elif varn == 's' or varn == 'QSNOW':
        varvals = ['s', 'snow_mixing_ratio', 0., 0.0003, 'snow|mixing|ratio',        \
          'kgkg-1', 'Purples']
    elif varn == 'stherm' or varn == 'LS_THERM':
        varvals = ['stherm', 'thermals_excess', 0., 0.8, 'thermals|excess', 'K',     \
          'Reds']
    elif varn == 'ta' or varn == 'WRFt' or varn == 'temp' or varn == 'LTEMP' or      \
      varn == 'Air temperature':
        varvals = ['ta', 'air_temperature', 195., 320., 'air|temperature', 'K',      \
          'YlOrRd']
    elif varn == 'tah' or varn == 'theta' or varn == 'LTHETA':
        varvals = ['tah', 'potential_air_temperature', 195., 320.,                   \
          'potential|air|temperature', 'K', 'YlOrRd']
    elif varn == 'tas' or varn == 'T2' or varn == 't2m' or varn == 'T2M' or          \
      varn == 'Temperature 2m':
        varvals = ['tas', 'air_temperature', 240., 310., 'air|temperature|at|2m', '  \
          K', 'YlOrRd']
    elif varn == 'tds' or varn == 'TH2':
        varvals = ['tds', 'air_dew_point_temperature', 240., 310.,                   \
          'air|dew|point|temperature|at|2m', 'K', 'YlGnBu']
    elif varn == 'tke' or varn == 'TKE' or varn == 'tke' or varn == 'LTKE':
        varvals = ['tke', 'turbulent_kinetic_energy', 0., 0.003,                     \
          'turbulent|kinetic|energy', 'm2/s2', 'Reds']
    elif varn == 'time'or varn == 'time_counter':
        varvals = ['time', 'time', 0., 1000., 'time',                                \
          'hours|since|1949/12/01|00:00:00', 'Reds']
    elif varn == 'tmla' or varn == 's_pblt' or varn == 'LS_PBLT':
        varvals = ['tmla', 'atmosphere_top_boundary_layer_temperature', 250., 330.,  \
          'atmosphere|top|boundary|layer|temperature', 'K', 'Reds']
    elif varn == 'ua' or varn == 'vitu' or varn == 'U' or varn == 'Zonal wind' or    \
      varn == 'LVITU':
        varvals = ['ua', 'eastward_wind', -30., 30., 'eastward|wind', 'ms-1',        \
          'seismic']
    elif varn == 'uas' or varn == 'u10m' or varn == 'U10' or varn =='Vent zonal 10m':
        varvals = ['uas', 'eastward_wind', -30., 30., 'eastward|2m|wind',    \
          'ms-1', 'seismic']
    elif varn == 'va' or varn == 'vitv' or varn == 'V' or varn == 'Meridional wind'  \
      or varn == 'LVITV':
        varvals = ['va', 'northward_wind', -30., 30., 'northward|wind', 'ms-1',      \
          'seismic']
    elif varn == 'vas' or varn == 'v10m' or varn == 'V10' or                         \
      varn =='Vent meridien 10m':
        varvals = ['vas', 'northward_wind', -30., 30., 'northward|2m|wind', 'ms-1',  \
          'seismic']
    elif varn == 'wakedeltaq' or varn == 'wake_deltaq' or varn == 'lwake_deltaq' or  \
      varn == 'LWAKE_DELTAQ':
        varvals = ['wakedeltaq', 'wake_delta_vapor', -0.003, 0.003,                  \
          'wake|delta|mixing|ratio', '-', 'seismic']
    elif varn == 'wakedeltat' or varn == 'wake_deltat' or varn == 'lwake_deltat' or  \
      varn == 'LWAKE_DELTAT':
        varvals = ['wakedeltat', 'wake_delta_temp', -0.003, 0.003,                   \
          'wake|delta|temperature', '-', 'seismic']
    elif varn == 'wakeh' or varn == 'wake_h' or varn == 'LWAKE_H':
        varvals = ['wakeh', 'wake_height', 0., 1000., 'height|of|the|wakes', 'm',    \
          'YlOrRd']
    elif varn == 'wakeomg' or varn == 'wake_omg' or varn == 'lwake_omg' or           \
      varn == 'LWAKE_OMG':
        varvals = ['wakeomg', 'wake_omega', 0., 3., 'wake|omega', \
          '-', 'BuGn']
    elif varn == 'wakes' or varn == 'wake_s' or varn == 'LWAKE_S':
        varvals = ['wakes', 'wake_area_fraction', 0., 0.5, 'wake|spatial|fraction',  \
          '1', 'BuGn']
    elif varn == 'wa' or varn == 'W' or varn == 'Vertical wind':
        varvals = ['wa', 'upward_wind', -10., 10., 'upward|wind', 'ms-1',            \
          'seismic']
    elif varn == 'wap' or varn == 'vitw' or varn == 'LVITW':
        varvals = ['wap', 'upward_wind', -3.e-10, 3.e-10, 'upward|wind', 'mPa-1',    \
          'seismic']
    elif varn == 'wss' or varn == 'SPDUV':
        varvals = ['wss', 'air_velocity',  0., 30., 'surface|horizontal|wind|speed', \
          'ms-1', 'Reds']
# Water budget
# Water budget de-accumulated
    elif varn == 'ccond' or varn == 'CCOND' or varn == 'ACCCONDde':
        varvals = ['ccond', 'cw_cond',  0., 30.,                                     \
          'cloud|water|condensation', 'mm', 'Reds']
    elif varn == 'wbr' or varn == 'ACQVAPORde':
        varvals = ['wbr', 'wbr',  0., 30., 'Water|Budget|water|wapor', 'mm', 'Blues']
    elif varn == 'diabh' or varn == 'DIABH' or varn == 'ACDIABHde':
        varvals = ['diabh', 'diabh',  0., 30., 'diabatic|heating', 'K', 'Reds']
    elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
        varvals = ['wbpw', 'water_budget_pw',  0., 30., 'Water|Budget|water|content',\
           'mms-1', 'Reds']
    elif varn == 'wbf' or varn == 'WBACF' or varn == 'WBACFde':
        varvals = ['wbf', 'water_budget_hfcqv',  0., 30.,                       \
          'Water|Budget|horizontal|convergence|of|water|vapour|(+,|' +   \
          'conv.;|-,|div.)', 'mms-1', 'Reds']
    elif varn == 'wbfc' or varn == 'WBFC' or varn == 'WBACFCde':
        varvals = ['wbfc', 'water_budget_fc',  0., 30.,                         \
          'Water|Budget|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
          'div.)', 'mms-1', 'Reds']
    elif varn == 'wbfp' or varn == 'WBFP' or varn == 'WBACFPde':
        varvals = ['wbfp', 'water_budget_cfp',  0., 30.,                        \
          'Water|Budget|horizontal|convergence|of|precipitation|(+,|' +  \
          'conv.;|-,|div.)', 'mms-1', 'Reds']
    elif varn == 'wbz' or varn == 'WBZ' or varn == 'WBACZde':
        varvals = ['wbz', 'water_budget_z',  0., 30.,                           \
          'Water|Budget|vertical|convergence|of|water|vapour|(+,|conv.' +\
          ';|-,|div.)', 'mms-1', 'Reds']
    elif varn == 'wbc' or varn == 'WBC' or varn == 'WBACCde':
        varvals = ['wbc', 'water_budget_c',  0., 30.,                           \
          'Water|Budget|Cloud|water|species','mms-1', 'Reds']
    elif varn == 'wbqvd' or varn == 'WBQVD' or varn == 'WBACQVDde':
        varvals = ['wbqvd', 'water_budget_qvd',  0., 30.,                       \
          'Water|Budget|water|vapour|divergence', 'mms-1', 'Reds']
    elif varn == 'wbqvblten' or varn == 'WBQVBLTEN' or varn == 'WBACQVBLTENde':
        varvals = ['wbqvblten', 'water_budget_qv_blten',  0., 30.,              \
          'Water|Budget|QV|tendency|due|to|pbl|parameterization',        \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbqvcuten' or varn == 'WBQVCUTEN' or varn == 'WBACQVCUTENde':
        varvals = ['wbqvcuten', 'water_budget_qv_cuten',  0., 30.,              \
          'Water|Budget|QV|tendency|due|to|cu|parameterization',         \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbqvshten' or varn == 'WBQVSHTEN' or varn == 'WBACQVSHTENde':
        varvals = ['wbqvshten', 'water_budget_qv_shten',  0., 30.,              \
          'Water|Budget|QV|tendency|due|to|shallow|cu|parameterization', \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbpr' or varn == 'WBP' or varn == 'WBACPde':
        varvals = ['wbpr', 'water_budget_pr',  0., 30.,                         \
          'Water|Budget|recipitation', 'mms-1', 'Reds']
    elif varn == 'wbpw' or varn == 'WBPW' or varn == 'WBACPWde':
        varvals = ['wbpw', 'water_budget_pw',  0., 30.,                         \
          'Water|Budget|water|content', 'mms-1', 'Reds']
    elif varn == 'wbcondt' or varn == 'WBCONDT' or varn == 'WBACCONDTde':
        varvals = ['wbcondt', 'water_budget_condt',  0., 30.,                   \
          'Water|Budget|condensation|and|deposition', 'mms-1', 'Reds']
    elif varn == 'wbqcm' or varn == 'WBQCM' or varn == 'WBACQCMde':
        varvals = ['wbqcm', 'water_budget_qcm',  0., 30.,                       \
          'Water|Budget|hydrometeor|change|and|convergence', 'mms-1', 'Reds']
    elif varn == 'wbsi' or varn == 'WBSI' or varn == 'WBACSIde':
        varvals = ['wbsi', 'water_budget_si',  0., 30.,                         \
          'Water|Budget|hydrometeor|sink', 'mms-1', 'Reds']
    elif varn == 'wbso' or varn == 'WBSO' or varn == 'WBACSOde':
        varvals = ['wbso', 'water_budget_so',  0., 30.,                         \
          'Water|Budget|hydrometeor|source', 'mms-1', 'Reds']
# Water Budget accumulated
    elif varn == 'ccondac' or varn == 'ACCCOND':
        varvals = ['ccondac', 'cw_cond_ac',  0., 30.,                                \
          'accumulated|cloud|water|condensation', 'mm', 'Reds']
    elif varn == 'rac' or varn == 'ACQVAPOR':
        varvals = ['rac', 'ac_r',  0., 30., 'accumualted|water|wapor', 'mm', 'Blues']
    elif varn == 'diabhac' or varn == 'ACDIABH':
        varvals = ['diabhac', 'diabh_ac',  0., 30., 'accumualted|diabatic|heating',  \
          'K', 'Reds']
    elif varn == 'wbpwac' or varn == 'WBACPW':
        varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
          'Water|Budget|accumulated|water|content', 'mm', 'Reds']
    elif varn == 'wbfac' or varn == 'WBACF':
        varvals = ['wbfac', 'water_budget_hfcqv_ac',  0., 30.,                       \
          'Water|Budget|accumulated|horizontal|convergence|of|water|vapour|(+,|' +   \
          'conv.;|-,|div.)', 'mm', 'Reds']
    elif varn == 'wbfcac' or varn == 'WBACFC':
        varvals = ['wbfcac', 'water_budget_fc_ac',  0., 30.,                         \
          'Water|Budget|accumulated|horizontal|convergence|of|cloud|(+,|conv.;|-,|' +\
          'div.)', 'mm', 'Reds']
    elif varn == 'wbfpac' or varn == 'WBACFP':
        varvals = ['wbfpac', 'water_budget_cfp_ac',  0., 30.,                        \
          'Water|Budget|accumulated|horizontal|convergence|of|precipitation|(+,|' +  \
          'conv.;|-,|div.)', 'mm', 'Reds']
    elif varn == 'wbzac' or varn == 'WBACZ':
        varvals = ['wbzac', 'water_budget_z_ac',  0., 30.,                           \
          'Water|Budget|accumulated|vertical|convergence|of|water|vapour|(+,|conv.' +\
          ';|-,|div.)', 'mm', 'Reds']
    elif varn == 'wbcac' or varn == 'WBACC':
        varvals = ['wbcac', 'water_budget_c_ac',  0., 30.,                           \
          'Water|Budget|accumulated|Cloud|water|species','mm', 'Reds']
    elif varn == 'wbqvdac' or varn == 'WBACQVD':
        varvals = ['wbqvdac', 'water_budget_qvd_ac',  0., 30.,                       \
          'Water|Budget|accumulated|water|vapour|divergence', 'mm', 'Reds']
    elif varn == 'wbqvbltenac' or varn == 'WBACQVBLTEN':
        varvals = ['wbqvbltenac', 'water_budget_qv_blten_ac',  0., 30.,              \
          'Water|Budget|accumulated|QV|tendency|due|to|pbl|parameterization',        \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbqvcutenac' or varn == 'WBACQVCUTEN':
        varvals = ['wbqvcutenac', 'water_budget_qv_cuten_ac',  0., 30.,              \
          'Water|Budget|accumulated|QV|tendency|due|to|cu|parameterization',         \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbqvshtenac' or varn == 'WBACQVSHTEN':
        varvals = ['wbqvshtenac', 'water_budget_qv_shten_ac',  0., 30.,              \
          'Water|Budget|accumulated|QV|tendency|due|to|shallow|cu|parameterization', \
          'kg kg-1 s-1', 'Reds']
    elif varn == 'wbprac' or varn == 'WBACP':
        varvals = ['wbprac', 'water_budget_pr_ac',  0., 30.,                         \
          'Water|Budget|accumulated|precipitation', 'mm', 'Reds']
    elif varn == 'wbpwac' or varn == 'WBACPW':
        varvals = ['wbpwac', 'water_budget_pw_ac',  0., 30.,                         \
          'Water|Budget|accumulated|water|content', 'mm', 'Reds']
    elif varn == 'wbcondtac' or varn == 'WBACCONDT':
        varvals = ['wbcondtac', 'water_budget_condt_ac',  0., 30.,                   \
          'Water|Budget|accumulated|condensation|and|deposition', 'mm', 'Reds']
    elif varn == 'wbqcmac' or varn == 'WBACQCM':
        varvals = ['wbqcmac', 'water_budget_qcm_ac',  0., 30.,                       \
          'Water|Budget|accumulated|hydrometeor|change|and|convergence', 'mm', 'Reds']
    elif varn == 'wbsiac' or varn == 'WBACSI':
        varvals = ['wbsiac', 'water_budget_si_ac',  0., 30.,                         \
          'Water|Budget|accumulated|hydrometeor|sink', 'mm', 'Reds']
    elif varn == 'wbsoac' or varn == 'WBACSO':
        varvals = ['wbsoac', 'water_budget_so_ac',  0., 30.,                         \
          'Water|Budget|accumulated|hydrometeor|source', 'mm', 'Reds']

    elif varn == 'xtime' or varn == 'XTIME':
        varvals = ['xtime', 'time',  0., 1.e5, 'time',                               \
          'minutes|since|simulation|start', 'Reds']
    elif varn == 'x' or varn == 'X':
        varvals = ['x', 'x',  0., 100., 'x', '-', 'Reds']
    elif varn == 'y' or varn == 'Y':
        varvals = ['y', 'y',  0., 100., 'y', '-', 'Blues']
    elif varn == 'z' or varn == 'Z':
        varvals = ['z', 'z',  0., 100., 'z', '-', 'Greens']
    elif varn == 'zg' or varn == 'WRFght' or varn == 'Geopotential height' or        \
      varn == 'geop' or varn == 'LGEOP':
        varvals = ['zg', 'geopotential_height', 0., 80000., 'geopotential|height',   \
          'm2s-2', 'rainbow']
    elif varn == 'zmaxth' or varn == 'zmax_th'  or varn == 'LZMAX_TH':
        varvals = ['zmaxth', 'thermal_plume_height', 0., 4000.,                     \
          'maximum|thermals|plume|height', 'm', 'YlOrRd']
    elif varn == 'zmla' or varn == 's_pblh' or varn == 'LS_PBLH':
        varvals = ['zmla', 'atmosphere_boundary_layer_thickness', 0., 2500.,         \
          'atmosphere|boundary|layer|thickness', 'm', 'Blues']
    else:
        print errormsg
        print '  ' + fname + ": variable '" + varn + "' not defined !!!"
        quit(-1)

    return varvals

def lonlat2D(lon,lat):
    """ Function to return lon, lat 2D matrices from any lon,lat matrix
      lon= matrix with longitude values
      lat= matrix with latitude values
    """
    fname = 'lonlat2D'

    if len(lon.shape) != len(lat.shape):
        print errormsg
        print '  ' + fname + ': longitude values with shape:', lon.shape,            \
          'is different that latitude values with shape:', lat.shape, '(dif. size) !!'
        quit(-1)

    if len(lon.shape) == 3:
        lonvv = lon[0,:,:]
        latvv = lat[0,:,:]
    elif len(lon.shape) == 2:
        lonvv = lon[:]
        latvv = lat[:]
    elif len(lon.shape) == 1:
        lonlatv = np.meshgrid(lon[:],lat[:])
        lonvv = lonlatv[0]
        latvv = lonlatv[1]

    return lonvv, latvv


def transform(vals, trans, dxv='0', dyv='0', dxt='0', dyt='0', dxl='0',         \
  dyl='0', dxtit='0', dytit='0'):
    """ Function to transform the values and the axes
      vals= values to transform
      trans= '|' separated list of operations of transformation
        'transpose': Transpose matrix of values (x-->y, y-->x)
        'flip@[x/y]': Flip the given axis
      d[x/y]v= original values for the [x/y]-axis
      d[x/y]t= original ticks for the [x/y]-axis
      d[x/y]l= original tick-labels for the [x/y]-axis
      d[x/y]tit= original titels for the [x/y]-axis

    >>> mat = np.arange(15).reshape(3,5)
    >>> xvals = np.zeros((3,5), dtype=np.float)
    >>> yvals = np.zeros((3,5), dtype=np.float)
    >>> for i in range(5):
    >>>     xvals[:,i] = -1.25 + i*0.5
    >>> for j in range(3):
    >>>     yvals[j,:] = 40.25 + j*0.5
    >>> xticks = np.array([-1.25, -0.25, 0.25])
    >>> yticks = np.array([40.25, 41.25])
    >>> xticklab = ['1.25 W', '0.25 W', '0,25 E']
    >>> yticklab = ['40.25 N', '41.25 N']
    >>> xtit = 'Longitude'
    >>> ytit = 'Latitude'
    >>> nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT = transform(mat, 'transpose', \
      xvals, yvals, xticks, yticks, xticklab, yticklab, xtit, ytit)
    [[ 0  5 10]
     [ 1  6 11]
     [ 2  7 12]
     [ 3  8 13]
     [ 4  9 14]]
    [[ 40.25  40.75  41.25]
     [ 40.25  40.75  41.25]
     [ 40.25  40.75  41.25]
     [ 40.25  40.75  41.25]
     [ 40.25  40.75  41.25]]
    [[-1.25 -1.25 -1.25]
     [-0.75 -0.75 -0.75]
     [-0.25 -0.25 -0.25]
     [ 0.25  0.25  0.25]
     [ 0.75  0.75  0.75]] 
    [ 40.25  41.25] 
    [-1.25 -0.25  0.25] 
    ['41.25 N', '40.25 N'] 
    ['0,25 E', '0.25 W', '1.25 W'] 
    Latitude 
    Longitude
    """
    fname = 'transform'

    transforms = trans.split('|')
    Ntransforms = len(transforms)

    newvals = vals.copy()
    wdvx = False
    if dxv != '0':
        newdxv = dxv.copy()
        wdvx = True
    else:
        newdxv = None
    wdvy = False
    if dyv != '0':
        newdyv = dyv.copy()
        wdvy = True
    else:
        newdyv = None

    wdtx = False
    if dxt != '0':
        newdxt = dxt.copy()
        wdtx = True
    else:
        newdxt = None
    wdty = False
    if dyt != '0':
        newdyt = dyt.copy()
        wdty = True
    else:
        newdyt = None

    wdlx = False
    if dxl != '0':
        newdxl = list(dxl)
        wdlx = True
    else:
        newdxl = None
    wdly = False
    if dyl != '0':
        newdyl = list(dyl)
        wdly = True
    else:
        newdyl = None

    wdTx = False
    if dxt != '0':
        newdxtit = str(dxtit)
        wdTx = True
    else:
        newdxtit = None
    wdTy = False
    if dyt != '0':
        newdytit = str(dytit)
        wdTy = True
    else:
        newdytit = None

    dx = newvals.shape[0]
    dy = newvals.shape[1]

    for transform in transforms:
        if transform == 'transpose':
            newvals = np.transpose(newvals)
            if wdvx and wdvy:
                copy = newdxv.copy()
                newdxv = np.transpose(newdyv)
                newdyv = np.transpose(copy)
#            else:
#                print errormsg
#                print '  ' + fname + ": it can not be '" + transform + "' without " +\
#                  ' both dimxvals:', dxv, ' and dimyvals:', dyv, ' !!'
#                quit(-1)

            if wdtx and wdtx:
                copy = newdxt.copy()
                newdxt = np.transpose(newdyt)
                newdyt = np.transpose(copy)
#            else:
#                print errormsg
#                print '  ' + fname + ": it can not be '" + transform + "' without " +\
#                  ' both dimxtickvals:', dxt, ' and dimytickvals:', dyt, ' !!'
#                quit(-1)

            if wdlx and wdly:
                copy = list(newdxl)
                newdxl = list(newdyl[::-1])
                newdyl = list(copy[::-1])
#            else:
#                print errormsg
#                print '  ' + fname + ": it can not be '" + transform + "' without " +\
#                  ' both dimxlabels:', dxl, ' and dimylabels:', dyl, ' !!'
#                quit(-1)

            if wdTx and wdTy:
                copy = str(newdxtit)
                newdxtit = str(newdytit)
                newdytit = str(copy)
#            else:
#                print errormsg
#                print '  ' + fname + ": it can not be '" + transform + "' without " +\
#                  ' both dimx title:', dxtit, ' and dimytitle:', dytit, ' !!'
#                quit(-1)

        elif transform[0:4] == 'flip':
            if transform.find('@') == -1:
                print errormsg
                print '  ' + fname + ": no ax given for '" + transform + "' !!"
                quit(-1)
                
            flip = transform.split('@')[1]
            if flip == 'x':
                newvals = newvals[...,::-1]
                if dv:
                    newdxv = newdxv[...,::-1]
                if dt:
                    newdxt = newdxt[::-1]
                if dl:
                    newdxl = newdxl[::-1]
            elif flip == 'y':
                newvals = newvals[...,::-1,:]
                if dv:
                    newdyv = newdyv[...,::-1,:]
                if dt:
                    newdyt = newdyt[::-1]
                if dl:
                    newdyl = newdyl[::-1]
            elif flip == 'z':
                newvals = newvals[...,::-1,:,:]
        else:
            print errormsg
            print '  '+ fname + ": transformation '" + transform + "' not ready!!"
            quit(-1)

    return newvals, newdxv, newdyv, newdxt, newdyt, newdxl, newdyl, newdxtit, newdytit

#mat = np.arange(15).reshape(3,5)
#xvals = np.zeros((3,5), dtype=np.float)
#yvals = np.zeros((3,5), dtype=np.float)
#for i in range(5):
#    xvals[:,i] = -1.25 + i*0.5
#for j in range(3):
#    yvals[j,:] = 40.25 + j*0.5
#xticks = np.array([-1.25, -0.25, 0.25])
#yticks = np.array([40.25, 41.25])
#xticklab = ['1.25 W', '0.25 W', '0,25 E']
#yticklab = ['40.25 N', '41.25 N']
#xtit = 'Longitude'
#ytit = 'Latitude'
#nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT = transform(mat, 'transpose', \
#      xvals, yvals, xticks, yticks, xticklab, yticklab, xtit, ytit)

#print nv, ndxv, ndyv, ndxt, ndyt, dxl, ndyl, ndxT, ndyT
#quit()

#######    #######    #######    #######    #######    #######    #######    #######    #######    #######

def check_colorBar(cbarn):
    """ Check if the given colorbar exists in matplotlib
    """
    fname = 'check_colorBar'

# Possible color bars
    colorbars = ['binary', 'Blues', 'BuGn', 'BuPu', 'gist_yarg', 'GnBu', 'Greens',   \
      'Greys', 'Oranges', 'OrRd', 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',       \
      'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'afmhot', 'autumn', 'bone',      \
      'cool', 'copper', 'gist_gray', 'gist_heat', 'gray', 'hot', 'pink', 'spring',   \
      'summer', 'winter', 'BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr', 'RdBu', \
      'RdGy', 'RdYlBu', 'RdYlGn', 'seismic', 'Accent', 'Dark2', 'hsv', 'Paired',     \
      'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3', 'spectral', 'gist_earth',        \
      'gist_ncar', 'gist_rainbow', 'gist_stern', 'jet', 'brg', 'CMRmap', 'cubehelix',\
      'gnuplot', 'gnuplot2', 'ocean', 'rainbow', 'terrain', 'flag', 'prism']

    if not searchInlist(colorbars,cbarn):
        print warnmsg
        print '  ' + fname + ' color bar: "' + cbarn + '" does not exist !!'
        print '  a standard one will be use instead !!'

    return

def units_lunits(u):
    """ Fucntion to provide LaTeX equivalences from a given units
      u= units to transform
    >>> units_lunits('kgkg-1')
    '$kgkg^{-1}$'    
    """
    fname = 'units_lunits'

    if u == 'h':
        print fname + '_____________________________________________________________'
        print units_lunits.__doc__
        quit()

# Units which does not change
    same = ['1', 'category', 'day', 'deg', 'degree', 'degrees', 'degrees East',      \
      'degrees Nord', 'degrees North', 'g', 'gpm', 'hour', 'hPa', 'K', 'Km', 'kg',   \
      'km', 'm', 'minute', 'mm', 'month', 'Pa', 's', 'second', 'um', 'year', '-']

    if searchInlist(same,u):
        lu = '$' + u + '$'
    elif len(u.split(' ')) > 1 and u.split(' ')[1] == 'since':
        uparts = u.split(' ')
        ip=0
        for up in uparts:
            if ip == 0:
               lu = '$' + up
            else:
               lu = lu + '\ ' + up
            ip=ip+1
        lu = lu + '$'
    else:
        if u == '': lu='-'
        elif u == 'C': lu='$^{\circ}C$'
        elif u == 'days': lu='$day$'
        elif u == 'degrees_East': lu='$degrees\ East$'
        elif u == 'degrees_east': lu='$degrees\ East$'
        elif u == 'degree_east': lu='$degrees\ East$'
        elif u == 'degrees longitude': lu='$degrees\ East$'
        elif u == 'degrees latitude': lu='$degrees\ North$'
        elif u == 'degrees_North': lu='$degrees\ North$'
        elif u == 'degrees_north': lu='$degrees\ North$'
        elif u == 'degree_north': lu='$degrees\ North$'
        elif u == 'deg C': lu='$^{\circ}C$'
        elif u == 'degC': lu='$^{\circ}C$'
        elif u == 'deg K': lu='$K$'
        elif u == 'degK': lu='$K$'
        elif u == 'g/g': lu='$gg^{-1}$'
        elif u == 'hours': lu='$hour$'
        elif u == 'J/kg': lu='$Jkg^{-1}$'
        elif u == 'Jkg-1': lu='$Jkg^{-1}$'
        elif u == 'K/m': lu='$Km^{-1}$'
        elif u == 'Km-1': lu='$Km^{-1}$'
        elif u == 'Km-3': lu='$Km^{-3}$'
        elif u == 'K/s': lu='$Ks^{-1}$'
        elif u == 'Ks-1': lu='$Ks^{-1}$'
        elif u == 'K s-1': lu='$Ks^{-1}$'
        elif u == 'kg/kg': lu='$kgkg^{-1}$'
        elif u == 'kgkg-1': lu='$kgkg^{-1}$'
        elif u == 'kg kg-1': lu='$kgkg^{-1}$'
        elif u == '(kg/kg)/s': lu='$kgkg^{-1}s^{-1}$'
        elif u == 'kgkg-1s-1': lu='$kgkg^{-1}s^{-1}$'
        elif u == 'kg kg-1 s-1': lu='$kgkg^{-1}s^{-1}$'
        elif u == 'kg/m2': lu='$kgm^{-2}$'
        elif u == 'kgm-2': lu='$kgm^{-2}$'
        elif u == 'kg m-2': lu='$kgm^{-2}$'
        elif u == 'Kg m-2': lu='$kgm^{-2}$'
        elif u == 'kg/m2/s': lu='$kgm^{-2}s^{-1}$'
        elif u == 'kg/m^2/s': lu='$kgm^{-2}s^{-1}$'
        elif u == 'kg/(m2*s)': lu='$kgm^{-2}s^{-1}$'
        elif u == 'kg/(s*m2)': lu='$kgm^{-2}s^{-1}$'
        elif u == 'kgm-2s-1': lu='$kgm^{-2}s^{-1}$'
        elif u == 'kg m-2 s-1': lu='$kgm^{-2}s^{-1}$'
        elif u == '1/m': lu='$m^{-1}$'
        elif u == 'm-1': lu='$m^{-1}$'
        elif u == 'm^2': lu='$m^{2}$'
        elif u == 'm2/s': lu='$m2s^{-1}$'
        elif u == 'm2s-1': lu='$m2s^{-1}$'
        elif u == 'm2/s2': lu='$m2s^{-2}$'
        elif u == 'm3/s': lu='$m^{3}s^{-1}$'
        elif u == 'm^3/s': lu='$m^{3}s^{-1}$'
        elif u == 'm/s': lu='$ms^{-1}$'
        elif u == 'mmh-3': lu='$mmh^{-3}$'
        elif u == 'ms-1': lu='$ms^{-1}$'
        elif u == 'm s-1': lu='$ms^{-1}$'
        elif u == 'm/s2': lu='$ms^{-2}$'
        elif u == 'ms-2': lu='$ms^{-2}$'
        elif u == 'minutes': lu='$minute$'
        elif u == 'Pa/s': lu='$Pas^{-1}$'
        elif u == 'Pas-1': lu='$Pas^{-1}$'
        elif u == 'W m-2': lu='$Wm^{-2}$'
        elif u == 'Wm-2': lu='$Wm^{-2}$'
        elif u == 'W/m2': lu='$Wm^{-2}$'
        elif u == '1/s': lu='$s^{-1}$'
        elif u == 's-1': lu='$s^{-1}$'
        elif u == 'seconds': lu='$second$'
        elif u == '%': lu='\%'
        elif u == '?': lu='-'
        else:
            print errormsg
            print '  ' + fname + ': units "' + u + '" not ready!!!!'
            quit(-1)

    return lu

def ASCII_LaTeX(ln):
    """ Function to transform from an ASCII line to LaTeX codification
      >>> ASCII_LaTeX('Laboratoire de Météorologie Dynamique però Hovmöller')
      Laboratoire de M\'et\'eorologie Dynamique per\`o Hovm\"oller
    """
    fname='ASCII_LaTeX'

    if ln == 'h':
        print fname + '_____________________________________________________________'
        print ASCII_LaTeX.__doc__
        quit()

    newln = ln.replace('\\', '\\textbackslash')

    newln = newln.replace('á', "\\'a")
    newln = newln.replace('é', "\\'e")
    newln = newln.replace('í', "\\'i")
    newln = newln.replace('ó', "\\'o")
    newln = newln.replace('ú', "\\'u")

    newln = newln.replace('à', "\\`a")
    newln = newln.replace('è', "\\`e")
    newln = newln.replace('ì', "\\`i")
    newln = newln.replace('ò', "\\`o")
    newln = newln.replace('ù', "\\`u")

    newln = newln.replace('â', "\\^a")
    newln = newln.replace('ê', "\\^e")
    newln = newln.replace('î', "\\^i")
    newln = newln.replace('ô', "\\^o")
    newln = newln.replace('û', "\\^u")

    newln = newln.replace('ä', '\\"a')
    newln = newln.replace('ë', '\\"e')
    newln = newln.replace('ï', '\\"i')
    newln = newln.replace('ö', '\\"o')
    newln = newln.replace('ü', '\\"u')

    newln = newln.replace('ç', '\c{c}')
    newln = newln.replace('ñ', '\~{n}')

    newln = newln.replace('Á', "\\'A")
    newln = newln.replace('É', "\\'E")
    newln = newln.replace('Í', "\\'I")
    newln = newln.replace('Ó', "\\'O")
    newln = newln.replace('Ú', "\\'U")

    newln = newln.replace('À', "\\`A")
    newln = newln.replace('È', "\\`E")
    newln = newln.replace('Ì', "\\`I")
    newln = newln.replace('Ò', "\\`O")
    newln = newln.replace('Ù', "\\`U")

    newln = newln.replace('Â', "\\^A")
    newln = newln.replace('Ê', "\\^E")
    newln = newln.replace('Î', "\\^I")
    newln = newln.replace('Ô', "\\^O")
    newln = newln.replace('Û', "\\^U")

    newln = newln.replace('Ä', '\\"A')
    newln = newln.replace('Ë', '\\"E')
    newln = newln.replace('Ï', '\\"I')
    newln = newln.replace('Ö', '\\"O')
    newln = newln.replace('Ü', '\\"U')

    newln = newln.replace('Ç', '\\c{C}')
    newln = newln.replace('Ñ', '\\~{N}')

    newln = newln.replace('¡', '!`')
    newln = newln.replace('¿', '¿`')
    newln = newln.replace('%', '\\%')
    newln = newln.replace('#', '\\#')
    newln = newln.replace('&', '\\&')
    newln = newln.replace('$', '\\$')
    newln = newln.replace('_', '\\_')
    newln = newln.replace('·', '\\textperiodcentered')
    newln = newln.replace('<', '$<$')
    newln = newln.replace('>', '$>$')
    newln = newln.replace('', '*')
#    newln = newln.replace('º', '$^{\\circ}$')
    newln = newln.replace('ª', '$^{a}$')
    newln = newln.replace('º', '$^{o}$')
    newln = newln.replace('°', '$^{\\circ}$')
    newln = newln.replace('\n', '\\\\\n')
    newln = newln.replace('\t', '\\medskip')

    return newln

def pretty_int(minv,maxv,Nint):
    """ Function to plot nice intervals
      minv= minimum value
      maxv= maximum value
      Nint= number of intervals
    >>> pretty_int(23.50,67.21,5)
    [ 25.  30.  35.  40.  45.  50.  55.  60.  65.]
    >>> pretty_int(-23.50,67.21,15)
    [  0.  20.  40.  60.]
    pretty_int(14.75,25.25,5)
    [ 16.  18.  20.  22.  24.]
    """ 
    fname = 'pretty_int'
    nice_int = [1,2,5]

#    print 'minv: ',minv,'maxv:',maxv,'Nint:',Nint

    interval = np.abs(maxv - minv)

    potinterval = np.log10(interval)
    Ipotint = int(potinterval)
    intvalue = np.float(interval / np.float(Nint))

# new
    potinterval = np.log10(intvalue)
    Ipotint = int(potinterval)

#    print 'interval:', interval, 'intavlue:', intvalue, 'potinterval:', potinterval, \
#     'Ipotint:', Ipotint, 'intvalue:', intvalue

    mindist = 10.e15
    for inice in nice_int:
#        print inice,':',inice*10.**Ipotint,np.abs(inice*10.**Ipotint - intvalue),mindist
        if np.abs(inice*10.**Ipotint - intvalue) < mindist:
            mindist = np.abs(inice*10.**Ipotint - intvalue)
            closestint = inice

    Ibeg = int(minv / (closestint*10.**Ipotint))

    values = []
    val = closestint*(Ibeg)*10.**(Ipotint)

#    print 'closestint:',closestint,'Ibeg:',Ibeg,'val:',val

    while val < maxv:
        values.append(val)
        val = val + closestint*10.**Ipotint

    return np.array(values, dtype=np.float)

def DegGradSec_deg(grad,deg,sec):
    """ Function to transform from a coordinate in grad deg sec to degrees (decimal)
    >>> DegGradSec_deg(39.,49.,26.)
    39.8238888889
    """
    fname = 'DegGradSec_deg'

    if grad == 'h':
        print fname + '_____________________________________________________________'
        print DegGradSec_deg.__doc__
        quit()

    deg = grad + deg/60. + sec/3600.

    return deg

def intT2dt(intT,tu):
    """ Function to provide an 'timedelta' object from a given interval value
      intT= interval value
      tu= interval units, [tu]= 'd': day, 'w': week, 'h': hour, 'i': minute, 's': second,
        'l': milisecond

      >>> intT2dt(3.5,'s')
      0:00:03.500000

      >>> intT2dt(3.5,'w')
      24 days, 12:00:00
    """
    import datetime as dt 

    fname = 'intT2dt'

    if tu == 'w':
        dtv = dt.timedelta(weeks=np.float(intT))
    elif tu == 'd':
        dtv = dt.timedelta(days=np.float(intT))
    elif tu == 'h':
        dtv = dt.timedelta(hours=np.float(intT))
    elif tu == 'i':
        dtv = dt.timedelta(minutes=np.float(intT))
    elif tu == 's':
        dtv = dt.timedelta(seconds=np.float(intT))
    elif tu == 'l':
        dtv = dt.timedelta(milliseconds=np.float(intT))
    else:
        print errormsg
        print '  ' + fname + ': time units "' + tu + '" not ready!!!!'
        quit(-1)

    return dtv

def lonlat_values(mapfile,lonvn,latvn):
    """ Function to obtain the lon/lat matrices from a given netCDF file
    lonlat_values(mapfile,lonvn,latvn)
      [mapfile]= netCDF file name
      [lonvn]= variable name with the longitudes
      [latvn]= variable name with the latitudes
    """

    fname = 'lonlat_values'

    if mapfile == 'h':
        print fname + '_____________________________________________________________'
        print lonlat_values.__doc__
        quit()

    if not os.path.isfile(mapfile):
        print errormsg
        print '  ' + fname + ": map file '" + mapfile + "' does not exist !!"
        quit(-1)

    ncobj = NetCDFFile(mapfile, 'r')
    lonobj = ncobj.variables[lonvn]
    latobj = ncobj.variables[latvn]

    if len(lonobj.shape) == 3:
        lonv = lonobj[0,:,:]
        latv = latobj[0,:,:]
    elif len(lonobj.shape) == 2:
        lonv = lonobj[:,:]
        latv = latobj[:,:]
    elif len(lonobj.shape) == 1:
        lon0 = lonobj[:]
        lat0 = latobj[:]
        lonv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
        latv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
        for iy in range(len(lat0)):
            lonv[iy,:] = lon0
        for ix in range(len(lon0)):
            latv[:,ix] = lat0
    else:
        print errormsg
        print '  ' + fname + ': lon/lat variables shape:',lonobj.shape,'not ready!!'
        quit(-1)

    return lonv, latv

def date_CFtime(ind,refd,tunits):
    """ Function to transform from a given date object a CF-convention time
      ind= date object to transform
      refd= reference date
      tunits= units for time
        >>> date_CFtime(dt.datetime(1976,02,17,08,30,00), dt.datetime(1949,12,01,00,00,00), 'seconds')
        827224200.0
    """
    import datetime as dt

    fname = 'date_CFtime'

    dt = ind - refd

    if tunits == 'weeks':
        value = dt.days/7. + dt.seconds/(3600.*24.*7.)
    elif tunits == 'days':
        value = dt.days + dt.seconds/(3600.*24.)
    elif tunits == 'hours':
        value = dt.days*24. + dt.seconds/(3600.)
    elif tunits == 'minutes':
        value = dt.days*24.*60. + dt.seconds/(60.)
    elif tunits == 'seconds':
        value = dt.days*24.*3600. + dt.seconds
    elif tunits == 'milliseconds':
        value = dt.days*24.*3600.*1000. + dt.seconds*1000.
    else:
        print errormsg
        print '  ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
        quit(-1)

    return value

def pot_values(values, uvals):
    """ Function to modify a seies of values by their potency of 10
      pot_values(values, uvals)
      values= values to modify
      uvals= units of the values
      >>> vals = np.sin(np.arange(20)*np.pi/5.+0.01)*10.e-5
      >>> pot_values(vals,'ms-1')
      (array([  0.00000000e+00,   5.87785252e-01,   9.51056516e-01,
         9.51056516e-01,   5.87785252e-01,   1.22464680e-16,
        -5.87785252e-01,  -9.51056516e-01,  -9.51056516e-01,
        -5.87785252e-01,  -2.44929360e-16,   5.87785252e-01,
         9.51056516e-01,   9.51056516e-01,   5.87785252e-01,
         3.67394040e-16,  -5.87785252e-01,  -9.51056516e-01,
        -9.51056516e-01,  -5.87785252e-01]), -4, 'x10e-4 ms-1', 'x10e-4')
    """

    fname = 'pot_values'

    if np.min(values) != 0.:
        potmin = int( np.log10( np.abs(np.min(values)) ) )
    else:
        potmin = 0

    if np.max(values) != 0.:
        potmax = int( np.log10( np.abs(np.max(values)) ) )
    else:
        potmax = 0

    if potmin * potmax > 9:
        potval = -np.min([np.abs(potmin), np.abs(potmax)]) * np.abs(potmin) / potmin

        newvalues = values*10.**potval
        potvalue = - potval
        potS = 'x10e' + str(potvalue)
        newunits = potS + ' ' + uvals
    else:
        newvalues = values
        potvalue = None
        potS = ''
        newunits = uvals

    return newvalues, potvalue, newunits, potS

def CFtimes_plot(timev,units,kind,tfmt):
    """ Function to provide a list of string values from a CF time values in order 
      to use them in a plot, according to the series of characteristics. 
      String outputs will be suited to the 'human-like' output
        timev= time values (real values)
        units= units string according to CF conventions ([tunits] since 
          [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]])
        kind= kind of output
          'Nval': according to a given number of values as 'Nval',[Nval]
          'exct': according to an exact time unit as 'exct',[tunit]; 
            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 
              'l': milisecond
        tfmt= desired format
          >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'Nval,5',"%Y/%m/%d %H:%M:%S")
            0.0 1979/12/01 00:00:00
            24.75 1979/12/02 00:45:00
            49.5 1979/12/03 01:30:00
            74.25 1979/12/04 02:15:00
            99.0 1979/12/05 03:00:00
          >>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'exct,2,d',"%Y/%m/%d %H:%M:%S")
            0.0 1979/12/01 00:00:00
            48.0 1979/12/03 00:00:00
            96.0 1979/12/05 00:00:00
            144.0 1979/12/07 00:00:00
    """ 
    import datetime as dt 

# Seconds between 0001 and 1901 Jan - 01
    secs0001_1901=59958144000.

    fname = 'CFtimes_plot'

    if timev == 'h':
        print fname + '_____________________________________________________________'
        print CFtimes_plot.__doc__
        quit()

    secsYear = 365.*24.*3600.
    secsWeek = 7.*24.*3600.
    secsDay = 24.*3600.
    secsHour = 3600.
    secsMinute = 60.
    secsMilisecond = 1./1000.
    secsMicrosecond = 1./1000000.

# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
##
    trefT = units.find(':')
    txtunits = units.split(' ')
    Ntxtunits = len(txtunits)

    if Ntxtunits == 3:
        Srefdate = txtunits[Ntxtunits - 1]
    else:
        Srefdate = txtunits[Ntxtunits - 2]

    if not trefT == -1:
#        print '  ' + fname + ': refdate with time!'
        if Ntxtunits == 3:
            refdate = datetimeStr_datetime(Srefdate)
        else:
            refdate = datetimeStr_datetime(Srefdate + '_' + txtunits[Ntxtunits - 1])
    else:
        refdate = datetimeStr_datetime(Srefdate + '_00:00:00')

    trefunits=units.split(' ')[0]
    if trefunits == 'weeks':
        trefu = 'w'
    elif trefunits == 'days':
        trefu = 'd'
    elif trefunits == 'hours':
        trefu = 'h'
    elif trefunits == 'minutes':
        trefu = 'm'
    elif trefunits == 'seconds':
        trefu = 's'
    elif trefunits == 'milliseconds':
        trefu = 'l'
    else:
        print errormsg
        print '  ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
        quit(-1)

    okind=kind.split(',')[0]
    dtv = len(timev)
  
    if refdate.year == 1:
        print warnmsg
        print '  ' + fname + ': changing reference date: ',refdate,                  \
          'to 1901-01-01_00:00:00 !!!'
        refdate = datetimeStr_datetime('1901-01-01_00:00:00')
        if trefu == 'w': timev = timev - secs0001_1901/(7.*24.*3600.)
        if trefu == 'd': timev = timev - secs0001_1901/(24.*3600.)
        if trefu == 'h': timev = timev - secs0001_1901/(3600.)
        if trefu == 'm': timev = timev - secs0001_1901/(60.)
        if trefu == 's': timev = timev - secs0001_1901
        if trefu == 'l': timev = timev - secs0001_1901*1000.
  
#    if refdate.year < 1901:
#        print warnmsg
#        print '  ' + fname + ': changing reference date: ',refdate,                  \
#          'to 1901-01-01_00:00:00 !!!'
#        diffrefs = dt.datetime(1901,1,1,0,0,0) - refdate
#        if gen.searchInlist(dir(diffrefs), 'total_seconds'):
#            secsref_1901 = diffrefs.total_seconds()
#        else:
#            secsref_1901 = diffrefs.days()*3600*24. + diffrefs.seconds()*1.
#        refdate = datetimeStr_datetime('1901-01-01_00:00:00')
#        print 'time 0 :', timev[0]
#        if trefu == 'w': timev = timev - secsref_1901/(7.*24.*3600.)
#        if trefu == 'd': timev = timev - secsref_1901/(24.*3600.)
#        if trefu == 'h': timev = timev - secsref_1901/(3600.)
#        if trefu == 'm': timev = timev - secsref_1901/(60.)
#        if trefu == 's': timev = timev - secsref_1901
#        if trefu == 'l': timev = timev - secsref_1901*1000.
#        print 'time 0 after:', timev[0]

    firstT = timev[0]
    lastT = timev[dtv-1]

# First and last times as datetime objects
    firstTdt = timeref_datetime(refdate, firstT, trefunits)
    lastTdt = timeref_datetime(refdate, lastT, trefunits)

# First and last times as [year, mon, day, hour, minut, second] vectors
    firstTvec = np.zeros((6), dtype= np.float)
    lastTvec = np.zeros((6), dtype= np.float)
    chTvec = np.zeros((6), dtype= bool)

    firstTvec = np.array([firstTdt.year, firstTdt.month, firstTdt.day, firstTdt.hour,\
      firstTdt.minute, firstTdt.second])
    lastTvec = np.array([lastTdt.year, lastTdt.month, lastTdt.day, lastTdt.hour,     \
      lastTdt.minute, lastTdt.second])

    chdate= lastTvec - firstTvec
    chTvec = np.where (chdate != 0., True, False)
   
    TOTdt = lastTdt - firstTdt
    TOTdtsecs = TOTdt.days*secsDay + TOTdt.seconds + TOTdt.microseconds*secsMicrosecond

    timeout = []
    if okind == 'Nval':
        nvalues = int(kind.split(',')[1])
        intervT = (lastT - firstT)/(nvalues-1)
        dtintervT = intT2dt(intervT, trefu)

        for it in range(nvalues):
            timeout.append(firstTdt + dtintervT*it)
    elif okind == 'exct':
        Nunits = int(kind.split(',')[1])
        tu = kind.split(',')[2]

# Generic incremental dt [seconds] according to all the possibilities ['c', 'y', 'm',
#   'w', 'd', 'h', 'i', 's', 'l'], some of them approximated (because they are not 
#   already necessary!)
        basedt = np.zeros((9), dtype=np.float)
        basedt[0] = (365.*100. + 25.)*24.*3600.
        basedt[1] = secsYear
        basedt[2] = 31.*24.*3600.
        basedt[3] = secsWeek
        basedt[4] = secsDay
        basedt[5] = secsHour
        basedt[6] = secsMinute
        basedt[7] = 1.
        basedt[8] = secsMilisecond

# Increment according to the units of the CF dates
        if trefunits == 'weeks':
            basedt = basedt/(secsWeek)
        elif trefunits == 'days':
            basedt = basedt/(secsDay)
        elif trefunits == 'hours':
            basedt = basedt/(secsHour)
        elif trefunits == 'minutes':
            basedt = basedt/(secsMinute)
        elif trefunits == 'seconds':
            basedt = basedt
        elif trefunits == 'milliseconds':
            basedt = basedt*secsMilisecond

        if tu == 'c':
            ti = firstTvec[0]
            tf = lastTvec[0] 
            centi = firstTvec[0] / 100

            datev = firstTdt
            while datev < lastTdt:
                yr = datev.year + Nunits*100
                mon = datev.month
                datev = dt.datetime(yr, mon, 1, 0, 0, 0)
                timeout.append(datev)

        elif tu == 'y':
            ti = firstTvec[0]
            tf = lastTvec[0]
            year = firstTvec[0]

            TOTsteps = int(TOTdtsecs/(Nunits*31*secsDay)) + 1

            if np.sum(firstTvec[1:5]) > 0:
                firstTdt = dt.datetime(year, 1, 1, 0, 0, 0)
                datev = dt.datetime(year, 1, 1, 0, 0, 0)

            datev = firstTdt
            timeout.append(datev)
            yr = datev.year
            while yr + Nunits <= tf :
                yr = datev.year + Nunits
                mon = datev.month
                datev = dt.datetime(yr, mon, 1, 0, 0, 0)
                timeout.append(datev)

        elif tu == 'm':
            ti = firstTvec[1]
            tf = lastTvec[1]
            
            yr = firstTvec[0]
            mon = firstTvec[1]

            TOTsteps = int(TOTdtsecs/(Nunits*31*secsDay)) + 1

            datev = firstTdt
            while datev < lastTdt:
                mon = datev.month + Nunits
                if mon > 12:
                    yr = yr + 1
                    mon = 1
                datev = dt.datetime(yr, mon, 1, 0, 0, 0)
                timeout.append(datev)

        elif tu == 'w':
            datev=firstTdt
            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(days=7*Nunits*it)
                timeout.append(datev)
                it = it + 1
        elif tu == 'd':
#            datev=firstTdt
            yr = firstTvec[0]
            mon = firstTvec[1]
            day = firstTvec[2]
 
            Iunits = np.mod(day,Nunits)
            if np.sum(firstTvec[2:5]) > 0:
                firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
                datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
            elif Iunits != 0:
                nNunits = int(day/Nunits)
                firstTdt = dt.datetime(yr, mon, nNunits*Nunits, 0, 0, 0)
                datev = dt.datetime(yr, mon, nNunits*Nunits, 0, 0, 0)
            else:
                firstTdt = dt.datetime(yr, mon, day, 0, 0, 0)
                datev = dt.datetime(yr, mon, day, 0, 0, 0)

            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(days=Nunits*it)
                timeout.append(datev)
                it = it + 1

        elif tu == 'h':
            datev=firstTdt
            yr = firstTvec[0]
            mon = firstTvec[1]
            day = firstTvec[2]
            hour = firstTvec[3]

            Iunits = np.mod(hour,Nunits)
            if np.sum(firstTvec[4:5]) > 0 or Iunits != 0:
                tadvance = 2*Nunits
                if tadvance >= 24:
                    firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
                    datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
                else:
                    nNunits = int(hour/Nunits)
                    firstTdt = dt.datetime(yr, mon, day, nNunits*Nunits, 0, 0)
                    datev = dt.datetime(yr, mon, day, nNunits*Nunits, 0, 0)
            else:
                firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
                datev = dt.datetime(yr, mon, day, hour, 0, 0)

            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(seconds=Nunits*3600*it)
                timeout.append(datev)
                it = it + 1                 
        elif tu == 'i':
            datev=firstTdt
            yr = firstTvec[0]
            mon = firstTvec[1]
            day = firstTvec[2]
            hour = firstTvec[3]
            minu = firstTvec[4]

            Iunits = np.mod(minu,Nunits)
            if firstTvec[5] > 0 or Iunits != 0:
                tadvance = 2*Nunits
                if tadvance >= 60:
                    firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
                    datev = dt.datetime(yr, mon, day, hour, 0, 0)
                else:
                    nNunits = int(minu/Nunits)
                    firstTdt = dt.datetime(yr, mon, day, hour, nNunits*Nunits, 0)
                    datev = dt.datetime(yr, mon, day, hour, nNunits*Nunits, 0)
            else:
                firstTdt = dt.datetime(yr, mon, day, hour, minu, 0)
                datev = dt.datetime(yr, mon, day, hour, minu, 0)
            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(seconds=Nunits*60*it)
                timeout.append(datev)
                it = it + 1                 
        elif tu == 's':
            datev=firstTdt
            yr = firstTvec[0]
            mon = firstTvec[1]
            day = firstTvec[2]
            hour = firstTvec[3]
            min = firstTvec[4]
            secu = firstTvec[5]

            Iunits = np.mod(secu,Nunits)
            if firstTvec[5] > 0 or Iunits != 0:
                tadvance = 2*Nunits
                if tadvance >= 60:
                    firstTdt = dt.datetime(yr, mon, day, hour, min, 0)
                    datev = dt.datetime(yr, mon, day, hour, min, 0)
                else:
                    nNunits = int(minu/Nunits)
                    firstTdt = dt.datetime(yr, mon, day, hour, min, nNunits*Nunits)
                    datev = dt.datetime(yr, mon, day, hour, min, nNunits*Nunits)
            else:
                firstTdt = dt.datetime(yr, mon, day, hour, min, secu)
                datev = dt.datetime(yr, mon, day, hour, min, secu)
            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(seconds=Nunits*it)
                timeout.append(datev)
                it = it + 1                 
        elif tu == 'l':
            datev=firstTdt
            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(seconds=Nunits*it/1000.)
                timeout.append(datev)
                it = it + 1
        else:
            print errormsg
            print '  '  + fname + ': exact units "' + tu + '" not ready!!!!!'
            quit(-1)

    else:
        print errormsg
        print '  ' + fname + ': output kind "' + okind + '" not ready!!!!'
        quit(-1)

    dtout = len(timeout)

    timeoutS = []
    timeoutv = np.zeros((dtout), dtype=np.float)

    for it in range(dtout):
        timeoutS.append(timeout[it].strftime(tfmt))
        timeoutv[it] = date_CFtime(timeout[it], refdate, trefunits)
       
#        print it,':',timeoutv[it], timeoutS[it]

    if len(timeoutv) < 1 or len(timeoutS) < 1:
        print errormsg
        print '  ' + fname + ': no time values are generated!'
        print '    values passed:',timev
        print '    units:',units
        print '    desired kind:',kind
        print '    format:',tfmt
        print '    function values ___ __ _'
        print '      reference date:',refdate
        print '      first date:',firstTdt
        print '      last date:',lastTdt
        print '      icrement:',basedt,trefunits

        quit(-1)

    return timeoutv, timeoutS

def color_lines(Nlines):
    """ Function to provide a color list to plot lines
    color_lines(Nlines)
      Nlines= number of lines
    """

    fname = 'color_lines'

    colors = ['r', 'b', 'g', 'p', 'g']

    colorv = []

    colorv.append('k')
    for icol in range(Nlines):
        colorv.append(colors[icol])


    return colorv

def output_kind(kindf, namef, close):
    """ Function to generate the output of the figure
      kindf= kind of the output 
        null: show in screen
        [jpg/pdf/png/ps]: standard output types
      namef= name of the figure (without extension)
      close= if the graph has to be close or not [True/False]
    """
    fname = 'output_kind'

    if kindf == 'h':
        print fname + '_____________________________________________________________'
        print output_kind.__doc__
        quit()

    if kindf == 'null':
        print 'showing figure...'
        plt.show()
    elif kindf == 'gif':
        plt.savefig(namef + ".gif")
        if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
    elif kindf == 'jpg':
        plt.savefig(namef + ".jpg")
        if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
    elif kindf == 'pdf':
        plt.savefig(namef + ".pdf")
        if close: print "Successfully generation of figure '" + namef + ".pdf' !!!"
    elif kindf == 'png':
        plt.savefig(namef + ".png")
        if close: print "Successfully generation of figure '" + namef + ".png' !!!"
    elif kindf == 'ps':
        plt.savefig(namef + ".ps")
        if close: print "Successfully generation of figure '" + namef + ".ps' !!!"
    else:
        print errormsg
        print '  ' + fname + ' output format: "' + kindf + '" not ready !!'
        print errormsg
        quit(-1)

    if close:
        plt.close()

    return

def check_arguments(funcname,args,expectargs,char):
    """ Function to check the number of arguments if they are coincident
    check_arguments(funcname,Nargs,args,char)
      funcname= name of the function/program to check
      args= passed arguments
      expectargs= expected arguments
      char= character used to split the arguments
    """

    fname = 'check_arguments'

    Nvals = len(args.split(char))
    Nargs = len(expectargs.split(char))

    if Nvals != Nargs:
        print errormsg
        print '  ' + fname + ': wrong number of arguments:',Nvals," passed to  '",   \
          funcname, "' which requires:",Nargs,'!!'
        print '     # given expected _______'
        Nmax = np.max([Nvals, Nargs])
        for ia in range(Nmax):
            if ia > Nvals-1:
                aval = ' '
            else:
                aval = args.split(char)[ia]
            if ia > Nargs-1:
                aexp = ' '
            else:
                aexp = expectargs.split(char)[ia]

            print '    ', ia, aval, aexp
        quit(-1)

    return

def Str_Bool(val):
    """ Function to transform from a String value to a boolean one
    >>> Str_Bool('True')
    True
    >>> Str_Bool('0')
    False
    >>> Str_Bool('no')
    False
    """

    fname = 'Str_Bool'

    if val == 'True' or val == 'true' or val == '1' or val == 'yes': 
        boolv = True
    elif val == 'False' or val == 'false' or val == '0' or val== 'no':
        boolv = False
    else:
        print errormsg
        print '  ' + fname + ": value '" + val + "' not ready!!"
        quit(-1)

    return boolv

def coincident_CFtimes(tvalB, tunitA, tunitB):
    """ Function to make coincident times for two different sets of CFtimes
    tvalB= time values B
    tunitA= time units times A to which we want to make coincidence
    tunitB= time units times B
    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
      'hours since 1949-12-01 00:00:00')
    [     0.   3600.   7200.  10800.  14400.  18000.  21600.  25200.  28800.  32400.]
    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
      'hours since 1979-12-01 00:00:00')
    [  9.46684800e+08   9.46688400e+08   9.46692000e+08   9.46695600e+08
       9.46699200e+08   9.46702800e+08   9.46706400e+08   9.46710000e+08
       9.46713600e+08   9.46717200e+08]
    """
    import datetime as dt
    fname = 'coincident_CFtimes'

    tunitA = tunitA.replace('_',' ')
    tunitB = tunitB.replace('_',' ')

    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
    tuA = tunitA.split(' ')[0]
    tuB = tunitB.split(' ')[0]

#    if tuB == 'months':
#        tvalB, tunitB  = gen.CFmonthU_daysU(tvalB, tunitB)
#        tuB = tunitB.split(' ')[0]

    if tuA != tuB:
        if tuA == 'microseconds':
            if tuB == 'microseconds':
                tB = tvalB*1.
            elif tuB == 'seconds':
                tB = tvalB*10.e6
            elif tuB == 'minutes':
                tB = tvalB*60.*10.e6
            elif tuB == 'hours':
                tB = tvalB*3600.*10.e6
            elif tuB == 'days':
                tB = tvalB*3600.*24.*10.e6
            else:
                print errormsg
                print '  ' + fname + ": combination of time untis: '" + tuA +        \
                  "' & '" + tuB + "' not ready !!"
                quit(-1)
        elif tuA == 'seconds':
            if tuB == 'microseconds':
                tB = tvalB/10.e6
            elif tuB == 'seconds':
                tB = tvalB*1.
            elif tuB == 'minutes':
                tB = tvalB*60.
            elif tuB == 'hours':
                tB = tvalB*3600.
            elif tuB == 'days':
                tB = tvalB*3600.*24.
            else:
                print errormsg
                print '  ' + fname + ": combination of time untis: '" + tuA +        \
                  "' & '" + tuB + "' not ready !!"
                quit(-1)
        elif tuA == 'minutes':
            if tuB == 'microseconds':
                tB = tvalB/(60.*10.e6)
            elif tuB == 'seconds':
                tB = tvalB/60.
            elif tuB == 'minutes':
                tB = tvalB*1.
            elif tuB == 'hours':
                tB = tvalB*60.
            elif tuB == 'days':
                tB = tvalB*60.*24.
            else:
                print errormsg
                print '  ' + fname + ": combination of time untis: '" + tuA +        \
                  "' & '" + tuB + "' not ready !!"
                quit(-1)
        elif tuA == 'hours':
            if tuB == 'microseconds':
                tB = tvalB/(3600.*10.e6)
            elif tuB == 'seconds':
                tB = tvalB/3600.
            elif tuB == 'minutes':
                tB = tvalB/60.
            elif tuB == 'hours':
                tB = tvalB*1.
            elif tuB == 'days':
                tB = tvalB*24.
            else:
                print errormsg
                print '  ' + fname + ": combination of time untis: '" + tuA +        \
                  "' & '" + tuB + "' not ready !!"
                quit(-1)
        elif tuA == 'days':
            if tuB == 'microseconds':
                tB = tvalB/(24.*3600.*10.e6)
            elif tuB == 'seconds':
                tB = tvalB/(24.*3600.)
            elif tuB == 'minutes':
                tB = tvalB/(24.*60.)
            elif tuB == 'hours':
                tB = tvalB/24.
            elif tuB == 'days':
                tB = tvalB*1.
            else:
                print errormsg
                print '  ' + fname + ": combination of time untis: '" + tuA +        \
                  "' & '" + tuB + "' not ready !!"
                quit(-1)
        else:
            print errormsg
            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
            quit(-1)
    else:
        tB = tvalB*1.

    if trefA != trefB:
        trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
        trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')

        difft = trefTB - trefTA
        diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds
        print '  ' + fname + ': different reference refA:',trefTA,'refB',trefTB
        print '    difference:',difft,':',diffv,'microseconds'

        if tuA == 'microseconds':
            tB = tB + diffv
        elif tuA == 'seconds':
            tB = tB + diffv/10.e6
        elif tuA == 'minutes':
            tB = tB + diffv/(60.*10.e6)
        elif tuA == 'hours':
            tB = tB + diffv/(3600.*10.e6)
        elif tuA == 'days':
            tB = tB + diffv/(24.*3600.*10.e6)
        else:
            print errormsg
            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
            quit(-1)

    return tB

####### ###### ##### #### ### ## #

def plot_TimeSeries(valtimes, vunits, tunits, hfileout, vtit, ttit, tkind, tformat,  \
  tit, linesn, lloc, kfig,coll,ptl,ptf):
    """ Function to draw time-series
      valtimes= list of arrays to plot [vals1[1values, 1times], [...,valsM[Mvals,Mtimes]]) 
      vunits= units of the values
      tunits= units of the times
      hfileout= header of the output figure. Final name: [hfileout]_[vtit].[kfig]
      vtit= variable title to be used in the graph
      ttit= time title to be used in the graph
      tkind= kind of time values to appear in the x-axis
        'Nval': according to a given number of values as 'Nval',[Nval]
        'exct': according to an exact time unit as 'exct',[tunit]; 
          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second, 
            'l': milisecond
      tformat= desired format of times
      tit= title of the graph
      linesn= list of values fot the legend
      lloc= location of the legend (0, autmoatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'
      kfig= type of figure: jpg, png, pds, ps
      coll= ',' list of colors for the lines, None for automatic, single
          value all the same
      ptl= ',' list of type of points for the lines, None for automatic, single
          value all the same
      ptf= frequency of point plotting, 'all' for all time steps
    """
    fname = 'plot_TimeSeries'

    if valtimes == 'h':
        print fname + '_____________________________________________________________'
        print plot_TimeSeries.__doc__
        quit()

    Nlines = len(valtimes)
# Canging line kinds every 7 lines (end of standard colors)
    Npts = len(pointkindsauto)
    linekinds = []
    pointkinds = []
    if ptl is None:
        if ptf is None:
            for ptype in range(Npts):
                for ip in range(7):
                    linekinds.append(pointkindsauto[ptype] + '-')
        else:
            for ptype in range(Npts):
                for ip in range(7):
                    linekinds.append('-')
                    pointkinds.append(pointkindsauto[ptype])
    else:
        if ptf is None:
            if len(ptl) > 1:
               for pt in ptl:
                    linekinds.append(pt + '-')
            else:
                for il in range(Nlines):
                    linekinds.append(ptl+'-')
        else:
            if len(ptl) > 1:
               for pt in ptl:
                    linekinds.append('-')
                    pointkinds.append(pt)
            else:
                for il in range(Nlines):
                    linekinds.append('-')
                pointkinds = ptl

    Nvalues = []
    Ntimes = []

    for il in range(Nlines):
        array = valtimes[il]

        if Nlines == 1:
            print warnmsg
            print '  ' + fname + ': drawing only one line!'

            Nvalues.append(array.shape[1])
            Ntimes.append(array.shape[0])
            tmin = np.min(array[1])
            tmax = np.max(array[1])
            vmin = np.min(array[0])
            vmax = np.max(array[0])
        else:
            Nvalues.append(array.shape[1])
            Ntimes.append(array.shape[0])
            tmin = np.min(array[1,:])
            tmax = np.max(array[1,:])
            vmin = np.min(array[0,:])
            vmax = np.max(array[0,:])

        if il == 0:
            xmin = tmin
            xmax = tmax
            ymin = vmin
            ymax = vmax
        else:
            if tmin < xmin: xmin = tmin
            if tmax > xmax: xmax = tmax
            if vmin < ymin: ymin = vmin
            if vmax > ymax: ymax = vmax

    dx = np.max(Ntimes)
    dy = np.min(Nvalues)

    plt.rc('text', usetex=True)

    print vtit
    if vtit == 'ps':
        plt.ylim(98000.,ymax)
    else:
        plt.ylim(ymin,ymax)

    plt.xlim(xmin,xmax)
#    print 'x lim:',xmin,xmax
#    print 'y lim:',ymin,ymax

    N7lines=0
    for il in range(Nlines):
        array = valtimes[il]
        if vtit == 'ps':
            array[0,:] = np.where(array[0,:] < 98000., None, array[0,:])

        if ptf is None:
            if coll is None:
                if ptf is None:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il])
                else:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il])
                    plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il],          \
                      label=linesn[il])
            else:
                if ptf is None:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \
                      color=coll[il])
                else:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \
                      color=coll[il])
                    plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il],          \
                      label=linesn[il], color=coll[il])
        else:
            if coll is None:
                if ptf is None:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il])
                else:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il])
                    plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il],          \
                      label=linesn[il])
            else:
                if ptf is None:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \
                      color=coll[il])
                else:
                    plt.plot(array[1,:],array[0,:], linekinds[il], label=linesn[il], \
                      color=coll[il])
                    plt.plot(array[1,::ptf],array[0,::ptf], pointkinds[il],          \
                      label=linesn[il], color=coll[il])

    timevals = np.arange(xmin,xmax)*1.

    tpos, tlabels = CFtimes_plot(timevals, tunits, tkind, tformat)

    if len(tpos) > 10:
        print warnmsg
        print '  ' + fname + ': with "' + tkind + '" there are', len(tpos), 'xticks !'

    plt.xticks(tpos, tlabels)
#    plt.Axes.set_xticklabels(tlabels)


    plt.legend(loc=lloc)
    plt.xlabel(ttit)
    plt.ylabel(vtit + " (" + vunits + ")")
    plt.title(tit.replace('_','\_').replace('&','\&'))

    figname = hfileout + '_' + vtit
    
    output_kind(kfig, figname, True)

    return

#Nt = 50
#Nlines = 3

#vtvalsv = []

#valsv = np.zeros((2,Nt), dtype=np.float)
## First
#valsv[0,:] = np.arange(Nt)
#valsv[1,:] = np.arange(Nt)*180.
#vtvalsv.append(valsv)
#del(valsv)

#valsv = np.zeros((2,Nt/2), dtype=np.float)
## Second
#valsv[0,:] = np.arange(Nt/2)
#valsv[1,:] = np.arange(Nt/2)*180.*2.
#vtvalsv.append(valsv)
#del(valsv)

#valsv = np.zeros((2,Nt/4), dtype=np.float)
## Third
#valsv[0,:] = np.arange(Nt/4)
#valsv[1,:] = np.arange(Nt/4)*180.*4.
#vtvalsv.append(valsv)
#del(valsv)

#varu='mm'
#timeu='seconds'

#title='test'
#linesname = ['line 1', 'line 2', 'line 3']

#plot_TimeSeries(vtvalsv, units_lunits(varu), timeu, 'test', 'vartest', 'time', title, linesname, 'png')
#quit()

def plot_points(xval, yval, vlon, vlat, extravals, extrapar, vtit, mapv, figk, color,\
  labels, lloc, kfig, figname):
    """ plotting points
      [x/yval]: x,y values to plot
      vlon= 2D-matrix with the longitudes
      vlat= 2D-matrix with the latitudes
      extravals= extra values to be added into the plot (None for nothing)
      extrapar= [varname, min, max, cbar, varunits] of the extra variable
      vtit= title of the graph ('|' for spaces)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert-conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      figK= kind of figure
        'legend': only points in the map with the legend with the names
        'labelled',[txtsize],[txtcol]: points with the names and size, color of text
      color= color for the points/labels ('auto', for "red")
      labels= list of labels for the points (None, no labels)
      lloc = localisation of the legend
      kfig= kind of figure (jpg, pdf, png)
      figname= name of the figure

    """
    fname = 'plot_points'
# Canging line kinds every 7 pts (end of standard colors)
    ptkinds=['.','x','o','*','+','8','>','D','h','p','s']

    Npts = len(xval)
    if Npts > len(ptkinds)*7:
        print errormsg
        print '  ' + fname + ': too many',Npts,'points!!'
        print "    enlarge 'ptkinds' list"
        quit(-1)

    N7pts = 0

    if color == 'auto':
        ptcol = "red"
    else:
        ptcol = color

    dx=vlon.shape[1]
    dy=vlat.shape[0]

    plt.rc('text', usetex=True)

    if not mapv is None:
#        vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:])
#        xvala = np.array(xval)
#        xvala = np.where(xvala < 0., 360. + xvala, xvala)
#        xval = list(xvala)

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = np.min(vlon)
        xlon = np.max(vlon)
        nlat = np.min(vlat)
        xlat = np.max(vlat)

        lon2 = vlon[dy/2,dx/2]
        lat2 = vlat[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": map projecion '" + map_proj + "' not ready!!"
            print '    available: cyl, lcc'
            quit(-1)

#        lons, lats = np.meshgrid(vlon, vlat)
#        lons = np.where(lons < 0., lons + 360., lons)

        x,y = m(vlon,vlat)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])
    else:
        x = vlon
        y = vlat
#        plt.xlim(0,dx-1)
#        plt.ylim(0,dy-1)

# Extra values
    if extravals is not None:
        plt.pcolormesh(x, y, extravals, cmap=plt.get_cmap(extrapar[3]),              \
          vmin=extrapar[1], vmax=extrapar[2])
        cbar = plt.colorbar()
        cbar.set_label(extrapar[0].replace('_','\_') +'('+ units_lunits(extrapar[4])+\
          ')')

    if labels is not None:
        for iv in range(len(xval)):
            if np.mod(iv,7) == 0: N7pts = N7pts + 1
#            print iv,xval[iv],yval[iv],labels[iv],ptkinds[N7pts]
            plt.plot(xval[iv],yval[iv], ptkinds[N7pts],label=labels[iv])

        if figk[0:8] == 'labelled':
            txtsize=int(figk.split(',')[1])
            txtcol=figk.split(',')[2]
            for iv in range(len(xval)):
                plt.annotate(labels[iv], xy=(xval[iv],yval[iv]), xycoords='data',    \
                  fontsize=txtsize, color=txtcol)
        elif figk == 'legend':
            plt.legend(loc=lloc)

    else: 
        plt.plot(xval, yval, '.', color=ptcol)

    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit.replace('|', ' '))
    
    output_kind(kfig, figname, True)

    return

def plot_list_points(xval, yval, varxn, varyn, vtit, figk, color, ptk, pts, labels, lloc, kfig, figname):
    """ plotting points
      [x/yval]: x,y values to plot
      var[x/y]n: names of the x,y variables
      vtit= title of the graph ('|' for spaces)
      figK= kind of figure
        'legend': only points in the map with the legend with the names
        'labelled',[txtsize],[txtcol]: points with the names and size, color of text
      color= color for the points/labels ('auto', for "red")
      ptk= kind of point
      pts= point size
      labels= list of labels for the points (None, no labels)
      lloc = localisation of the legend
      kfig= kind of figure (jpg, pdf, png)
      figname= name of the figure

    """
    fname = 'plot_points'
# Canging line kinds every 7 pts (end of standard colors)
    ptkinds=['.','x','o','*','+','8','>','D','h','p','s']

    if color == 'auto':
        ptcol = "red"
    else:
        ptcol = color

    plt.rc('text', usetex=True)

    if labels is not None:
        for iv in range(len(xval)):
#            print iv,xval[iv],yval[iv],labels[iv],ptkinds[N7pts]
            plt.plot(xval[iv],yval[iv], ptk, markersize=pts, label=labels[iv])

        if figk[0:8] == 'labelled':
            txtsize=int(figk.split(',')[1])
            txtcol=figk.split(',')[2]
            for iv in range(len(xval)):
                plt.annotate(labels[iv], xy=(xval[iv],yval[iv]), xycoords='data',    \
                  fontsize=txtsize, color=txtcol)
        elif figk == 'legend':
            plt.legend(loc=lloc)

    else: 
        plt.plot(xval, yval, ptk, markersize=pts, color=ptcol)

    plt.xlabel(varxn)
    plt.ylabel(varyn)

    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit.replace('|', ' '))
    
    output_kind(kfig, figname, True)

    return

def plot_2Dfield(varv,dimn,colorbar,vn,vx,unit,olon,olat,ifile,vtit,zvalue,time,tk,  \
   tkt,tobj,tvals,tind,kfig,mapv,reva):
    """ Adding labels and other staff to the graph
      varv= 2D values to plot
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      vn,vm= minmum and maximum values to plot
      unit= units of the variable
      olon,olat= longitude, latitude objects
      ifile= name of the input file
      vtit= title of the variable
      zvalue= value on the z axis
      time= value on the time axis
      tk= kind of time (WRF)
      tkt= kind of time taken
      tobj= tim object 
      tvals= values of the time variable
      tind= time index
      kfig= kind of figure (jpg, pdf, png)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      reva= reverse the axes (x-->y, y-->x)
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt

    if reva:
        print '  reversing the axes of the figure (x-->y, y-->x)!!'
        varv = np.transpose(varv)
        dimn0 = []
        dimn0.append(dimn[1] + '')
        dimn0.append(dimn[0] + '')
        dimn = dimn0

    fname = 'plot_2Dfield'
    dx=varv.shape[1]
    dy=varv.shape[0]

    plt.rc('text', usetex=True)
#    plt.rc('font', family='serif')

    if not mapv is None:
        if len(olon[:].shape) == 3:
            lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
            lat0 = olat[0,]
        elif len(olon[:].shape) == 2:
            lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
            lat0 = olat[:]
        elif len(olon[:].shape) == 1:
            lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
            lat00 = olat[:]
            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

            for iy in range(len(lat00)):
                lon0[iy,:] = lon00
            for ix in range(len(lon00)):
                lat0[:,ix] = lat00

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = lon0[0,0]
        xlon = lon0[dy-1,dx-1]
        nlat = lat0[0,0]
        xlat = lat0[dy-1,dx-1]

        lon2 = lon0[dy/2,dx/2]
        lat2 = lat0[dy/2,dx/2]

        print '    lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)

        if len(olon[:].shape) == 1:
            lons, lats = np.meshgrid(olon[:], olat[:])
        else:
            lons = olon[0,:]
            lats = olat[:,0]
 
        lons = np.where(lons < 0., lons + 360., lons)

        x,y = m(lons,lats)
        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
        cbar = plt.colorbar()

        m.drawcoastlines()
#        if (nlon > 180. or xlon > 180.):
#            nlon0 = nlon
#            xlon0 = xlon
#            if (nlon > 180.): nlon0 = nlon - 360.
#            if (xlon > 180.): xlon0 = xlon - 360.
#            meridians = pretty_int(nlon0,xlon0,5)           
#            meridians = np.where(meridians < 0., meridians + 360., meridians)
#        else:
#            meridians = pretty_int(nlon,xlon,5)

        meridians = pretty_int(nlon,xlon,5)

        m.drawmeridians(meridians,labels=[True,False,False,True])
        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

    else:
        plt.xlim(0,dx-1)
        plt.ylim(0,dy-1)

        plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
        cbar = plt.colorbar()

        plt.xlabel(dimn[1].replace('_','\_'))
        plt.ylabel(dimn[0].replace('_','\_'))

# set the limits of the plot to the limits of the data
#    plt.axis([x.min(), x.max(), y.min(), y.max()])

#    plt.plot(varv)
    cbar.set_label(unit)

    figname = ifile.replace('.','_') + '_' + vtit
    graphtit = vtit.replace('_','\_').replace('&','\&')

    if zvalue != 'null':
        graphtit = graphtit + ' at z= ' + zvalue
        figname = figname + '_z' + zvalue
    if tkt == 'tstep':
        graphtit = graphtit + ' at time-step= ' + time.split(',')[1]
        figname = figname + '_t' + time.split(',')[1].zfill(4)
    elif tkt == 'CFdate':
        graphtit = graphtit + ' at ' + tobj.strfmt("%Y%m%d%H%M%S")
        figname = figname + '_t' + tobj.strfmt("%Y%m%d%H%M%S")

    if tk == 'WRF':
#        datev = str(timevals[timeind][0:9])
        datev = tvals[tind][0] + tvals[tind][1] + tvals[tind][2] + \
          timevals[timeind][3] + timevals[timeind][4] + timevals[timeind][5] +       \
          timevals[timeind][6] + timevals[timeind][7] + timevals[timeind][8] +       \
          timevals[timeind][9]
#        timev = str(timevals[timeind][11:18])
        timev = timevals[timeind][11] + timevals[timeind][12] +                      \
          timevals[timeind][13] + timevals[timeind][14] + timevals[timeind][15] +    \
          timevals[timeind][16] + timevals[timeind][17] + timevals[timeind][18]
        graphtit = vtit.replace('_','\_') + ' (' + datev + ' ' + timev + ')'

    plt.title(graphtit)
    
    output_kind(kfig, figname, True)

    return

def plot_2Dfield_easy(varv,dimxv,dimyv,dimn,colorbar,vn,vx,unit,ifile,vtit,kfig,reva):
    """ Adding labels and other staff to the graph
      varv= 2D values to plot
      dim[x/y]v = values at the axes of x and y
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      vn,vm= minmum and maximum values to plot
      unit= units of the variable
      ifile= name of the input file
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= reverse the axes (x-->y, y-->x)
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt
    fname = 'plot_2Dfield'

    if reva:
        print '  reversing the axes of the figure (x-->y, y-->x)!!'
        varv = np.transpose(varv)
        dimn0 = []
        dimn0.append(dimn[1] + '')
        dimn0.append(dimn[0] + '')
        dimn = dimn0
        if len(dimyv.shape) == 2:
            x = np.transpose(dimyv)
        else:
            if len(dimxv.shape) == 2:
                ddx = len(dimyv)
                ddy = dimxv.shape[1]
            else:
                ddx = len(dimyv)
                ddy = len(dimxv)

            x = np.zeros((ddy,ddx), dtype=np.float)
            for j in range(ddy):
                x[j,:] = dimyv

        if len(dimxv.shape) == 2:
            y = np.transpose(dimxv)
        else:
            if len(dimyv.shape) == 2:
                ddx = dimyv.shape[0]
                ddy = len(dimxv)
            else:
                ddx = len(dimyv)
                ddy = len(dimxv)

            y = np.zeros((ddy,ddx), dtype=np.float)
            for i in range(ddx):
                y[:,i] = dimxv
    else:
        if len(dimxv.shape) == 2:
            x = dimxv
        else:
            x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
            for j in range(len(dimyv)):
                x[j,:] = dimxv

        if len(dimyv.shape) == 2:
            y = dimyv
        else:
            y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
            for i in range(len(dimxv)):
                x[:,i] = dimyv

    dx=varv.shape[1]
    dy=varv.shape[0]

    plt.rc('text', usetex=True)
    plt.xlim(0,dx-1)
    plt.ylim(0,dy-1)

    plt.pcolormesh(x, y, varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
#    plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
    cbar = plt.colorbar()

    plt.xlabel(dimn[1].replace('_','\_'))
    plt.ylabel(dimn[0].replace('_','\_'))

# set the limits of the plot to the limits of the data
    plt.axis([x.min(), x.max(), y.min(), y.max()])
#    if varv.shape[1] / varv.shape[0] > 10:
#        plt.axes().set_aspect(0.001)
#    else:
#        plt.axes().set_aspect(np.float(varv.shape[0])/np.float(varv.shape[1]))

    cbar.set_label(unit)

    figname = ifile.replace('.','_') + '_' + vtit
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)

    output_kind(kfig, figname, True)
    
    return

def plot_Trajectories(lonval, latval, linesn, olon, olat, lonlatLims, gtit, kfig,    \
  mapv, obsname):
    """ plotting Trajectories
      [lon/latval]= lon,lat values to plot (as list of vectors)
      linesn: name of the lines
      o[lon/lat]= object with the longitudes and the latitudes of the map to plot
      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
      gtit= title of the graph
      kfig= kind of figure (jpg, pdf, png)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      obsname= name of the observations in graph (can be None for without). 
        Observational trajectory would be the last one
    """
    fname = 'plot_Trajectories'

    if lonval == 'h':
        print fname + '_____________________________________________________________'
        print plot_Trajectories.__doc__
        quit()

# Canging line kinds every 7 lines (end of standard colors)
    linekinds=['.-','x-','o-']

    Ntraj = len(lonval)

    if obsname is not None:
        Ntraj = Ntraj - 1

    N7lines = 0

    plt.rc('text', usetex=True)

    if not mapv is None:
        if len(olon[:].shape) == 3:
#            lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
            lon0 = olon[0,]
            lat0 = olat[0,]
        elif len(olon[:].shape) == 2:
#            lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
            lon0 = olon[:]
            lat0 = olat[:]
        elif len(olon[:].shape) == 1:
#            lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
            lon00 = olon[:]
            lat00 = olat[:]
            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

            for iy in range(len(lat00)):
                lon0[iy,:] = lon00
            for ix in range(len(lon00)):
                lat0[:,ix] = lat00

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        dx = lon0.shape[1]
        dy = lon0.shape[0]

        nlon = lon0[0,0]
        xlon = lon0[dy-1,dx-1]
        nlat = lat0[0,0]
        xlat = lat0[dy-1,dx-1]

        lon2 = lon0[dy/2,dx/2]
        lat2 = lat0[dy/2,dx/2]

        if lonlatLims is not None:
            plt.xlim(lonlatLims[0], lonlatLims[2])
            plt.ylim(lonlatLims[1], lonlatLims[3])
            if map_proj == 'cyl':
                nlon = lonlatLims[0]
                nlat = lonlatLims[1]
                xlon = lonlatLims[2]
                xlat = lonlatLims[3]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)

        if len(olon.shape) == 3:
#            lons, lats = np.meshgrid(olon[0,:,:], olat[0,:,:])
            lons = olon[0,:,:]
            lats = olat[0,:,:]

        elif len(olon.shape) == 2:
#            lons, lats = np.meshgrid(olon[:,:], olat[:,:])
            lons = olon[:,:]
            lats = olat[:,:]
        else:
            dx = olon.shape
            dy = olat.shape
#            print errormsg
#            print '  ' + fname + ': shapes of lon/lat objects', olon.shape,          \
#              'not ready!!!'

        for il in range(Ntraj):
            plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
            if il == 6: N7lines = N7lines + 1

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    else:
        if len(olon.shape) == 3:
            dx = olon.shape[2]
            dy = olon.shape[1]
        elif len(olon.shape) == 2:
            dx = olon.shape[1]
            dy = olon.shape[0]
        else:
            dx = olon.shape
            dy = olat.shape
#            print errormsg
#            print '  ' + fname + ': shapes of lon/lat objects', olon.shape,          \
#              'not ready!!!'

        if lonlatLims is not None:
            plt.xlim(lonlatLims[0], lonlatLims[2])
            plt.ylim(lonlatLims[1], lonlatLims[3])
        else:
            plt.xlim(np.min(olon[:]),np.max(olon[:]))
            plt.ylim(np.min(olat[:]),np.max(olat[:]))

        for il in range(Ntraj):
            plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
            if il == 6: N7lines = N7lines + 1

        plt.xlabel('x-axis')
        plt.ylabel('y-axis')

    figname = 'trajectories'
    graphtit = gtit

    if obsname is not None:
        plt.plot(lonval[Ntraj], latval[Ntraj], linestyle='-', color='k',             \
          linewidth=3, label= obsname)

    plt.title(graphtit.replace('_','\_').replace('&','\&'))
    plt.legend()
    
    output_kind(kfig, figname, True)

    return

def plot_topo_geogrid(varv, olon, olat, mint, maxt, lonlatLims, gtit, kfig, mapv,    \
  closeif):
    """ plotting geo_em.d[nn].nc topography from WPS files
    plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
      varv= topography values
      o[lon/lat]= longitude and latitude objects
      [min/max]t: minimum and maximum values of topography to draw
      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
      gtit= title of the graph
      kfig= kind of figure (jpg, pdf, png)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lamvbert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      closeif= Boolean value if the figure has to be closed
    """
    fname = 'plot_topo_geogrid'

    if varv == 'h':
        print fname + '_____________________________________________________________'
        print plot_topo_geogrid.__doc__
        quit()

    dx=varv.shape[1]
    dy=varv.shape[0]

    plt.rc('text', usetex=True)
#    plt.rc('font', family='serif')

    if not mapv is None:
        if len(olon[:].shape) == 3:
            lon0 = olon[0,]
            lat0 = olat[0,]
        elif len(olon[:].shape) == 2:
            lon0 = olon[:]
            lat0 = olat[:]
        elif len(olon[:].shape) == 1:
            lon00 = olon[:]
            lat00 = olat[:]
            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

            for iy in range(len(lat00)):
                lon0[iy,:] = lon00
            for ix in range(len(lon00)):
                lat0[:,ix] = lat00

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]
        dx = lon0.shape[1]
        dy = lon0.shape[0]

        if lonlatLims is not None:
            print '  ' + fname + ': cutting the domain to plot !!!!'
            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
            nlon =  lonlatLims[0]
            xlon =  lonlatLims[2]
            nlat =  lonlatLims[1]
            xlat =  lonlatLims[3]

            if map_proj == 'lcc':
                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
        else:
            nlon = lon0[0,0]
            xlon = lon0[dy-1,dx-1]
            nlat = lat0[0,0]
            xlat = lat0[dy-1,dx-1]
            lon2 = lon0[dy/2,dx/2]
            lat2 = lat0[dy/2,dx/2]

        plt.xlim(nlon, xlon)
        plt.ylim(nlat, xlat)
        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": map projection '" + map_proj + "' not ready !!"
            quit(-1)

        if len(olon[:].shape) == 1:
            lons, lats = np.meshgrid(olon[:], olat[:])
        else:
            if len(olon[:].shape) == 3:
                lons = olon[0,:,:]
                lats = olat[0,:,:]
            else:
                lons = olon[:]
                lats = olat[:]
 
        x,y = m(lons,lats)

        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
        cbar = plt.colorbar()

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')
    else:
        print emsg
        print '  ' + fname + ': A projection parameter is needed None given !!'
        quit(-1)    

    figname = 'domain'
    graphtit = gtit.replace('_','\_')
    cbar.set_label('height ($m$)')

    plt.title(graphtit.replace('_','\_').replace('&','\&'))
    
    output_kind(kfig, figname, closeif)

    return

def plot_topo_geogrid_boxes(varv, boxesX, boxesY, boxlabels, olon, olat, mint, maxt, \
  lonlatLims, gtit, kfig, mapv, gloc, closeif):
    """ plotting geo_em.d[nn].nc topography from WPS files
    plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
      varv= topography values
      boxesX/Y= 4-line sets to draw the boxes
      boxlabels= labels for the legend of the boxes
      o[lon/lat]= longitude and latitude objects
      [min/max]t: minimum and maximum values of topography to draw
      lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
      gtit= title of the graph
      kfig= kind of figure (jpg, pdf, png)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lamvbert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      gloc= location of the legend (0, autmoatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'
      closeif= Boolean value if the figure has to be closed
    """
    fname = 'plot_topo_geogrid'

    if varv == 'h':
        print fname + '_____________________________________________________________'
        print plot_topo_geogrid.__doc__
        quit()

    cols = color_lines(len(boxlabels))

    dx=varv.shape[1]
    dy=varv.shape[0]

    plt.rc('text', usetex=True)
#    plt.rc('font', family='serif')

    if not mapv is None:
        if len(olon[:].shape) == 3:
            lon0 = olon[0,]
            lat0 = olat[0,]
        elif len(olon[:].shape) == 2:
            lon0 = olon[:]
            lat0 = olat[:]
        elif len(olon[:].shape) == 1:
            lon00 = olon[:]
            lat00 = olat[:]
            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

            for iy in range(len(lat00)):
                lon0[iy,:] = lon00
            for ix in range(len(lon00)):
                lat0[:,ix] = lat00

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]
        dx = lon0.shape[1]
        dy = lon0.shape[0]

        if lonlatLims is not None:
            print '  ' + fname + ': cutting the domain to plot !!!!'
            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
            print '    limits: N-S', lonlatLims[1], lonlatLims[3]
            nlon =  lonlatLims[0]
            xlon =  lonlatLims[2]
            nlat =  lonlatLims[1]
            xlat =  lonlatLims[3]

            if map_proj == 'lcc':
                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
        else:
            nlon = np.min(lon0)
            xlon = np.max(lon0)
            nlat = np.min(lat0)
            xlat = np.max(lat0)
            lon2 = lon0[dy/2,dx/2]
            lat2 = lat0[dy/2,dx/2]

        plt.xlim(nlon, xlon)
        plt.ylim(nlat, xlat)
        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' does not exist!!"
            print '    existing ones: cyl, lcc'
            quit(-1)

        if len(olon[:].shape) == 1:
            lons, lats = np.meshgrid(olon[:], olat[:])
        else:
            if len(olon[:].shape) == 3:
                lons = olon[0,:,:]
                lats = olat[0,:,:]
            else:
                lons = olon[:]
                lats = olat[:]
 
        x,y = m(lons,lats)

        plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
        cbar = plt.colorbar()

        Nboxes = len(boxesX)/4
        for ibox in range(Nboxes):
            plt.plot(boxesX[ibox*4], boxesY[ibox*4], linestyle='-', linewidth=3,     \
              label=boxlabels[ibox], color=cols[ibox])
            plt.plot(boxesX[ibox*4+1], boxesY[ibox*4+1], linestyle='-', linewidth=3, \
              color=cols[ibox])
            plt.plot(boxesX[ibox*4+2], boxesY[ibox*4+2], linestyle='-', linewidth=3, \
              color=cols[ibox])
            plt.plot(boxesX[ibox*4+3], boxesY[ibox*4+3], linestyle='-', linewidth=3, \
              color=cols[ibox])

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')
    else:
        print emsg
        print '  ' + fname + ': A projection parameter is needed None given !!'
        quit(-1)    

    figname = 'domain_boxes'
    graphtit = gtit.replace('_','\_').replace('&','\&')
    cbar.set_label('height ($m$)')

    plt.title(graphtit)
    plt.legend(loc=gloc)
    
    output_kind(kfig, figname, closeif)

    return

def plot_2D_shadow(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,          \
  colorbar,vs,uts,vtit,kfig,reva,mapv,ifclose):
    """ Adding labels and other staff to the graph
      varsv= 2D values to plot with shading
      vnames= variable names for the figure
      dim[x/y]v = values at the axes of x and y
      dim[x/y]u = units at the axes of x and y
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      vs= minmum and maximum values to plot in shadow or:
        'Srange': for full range
        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
        'Saroundminmax@val': for min*val,max*val
        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
          percentile_(100-val)-median)
        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
           percentile_(100-val)-median)
      uts= units of the variable to shadow
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= ('|' for combination)
        * 'transpose': reverse the axes (x-->y, y-->x)
        * 'flip'@[x/y]: flip the axis x or y
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      ifclose= boolean value whether figure should be close (finish) or not
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt
    fname = 'plot_2D_shadow'

#    print dimyv[73,21]
#    dimyv[73,21] = -dimyv[73,21]
#    print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv)

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_2D_shadow.__doc__
        quit()

    if len(varsv.shape) != 2:
        print errormsg
        print '  ' + fname + ': wrong variable shape:',varsv.shape,'is has to be 2D!!'
        quit(-1)

    reva0 = ''
    if reva.find('|') != 0:
        revas = reva.split('|')
    else:
        revas = [reva]
        reva0 = reva

    for rev in revas:
        if reva[0:4] == 'flip':
            reva0 = 'flip'
            if len(reva.split('@')) != 2:
                 print errormsg
                 print '  ' + fname + ': flip is given', reva, 'but not axis!'
                 quit(-1)

        if rev == 'transpose':
            print '  reversing the axes of the figure (x-->y, y-->x)!!'
            varsv = np.transpose(varsv)
            dxv = dimyv
            dyv = dimxv
            dimxv = dxv
            dimyv = dyv

    if len(dimxv[:].shape) == 3:
        xdims = '1,2'
    elif len(dimxv[:].shape) == 2:
        xdims = '0,1'
    elif len(dimxv[:].shape) == 1:
        xdims = '0'

    if len(dimyv[:].shape) == 3:
        ydims = '1,2'
    elif len(dimyv[:].shape) == 2:
        ydims = '0,1'
    elif len(dimyv[:].shape) == 1:
        ydims = '0'

#    lon0 = dimxv
#    lat0 = dimyv
    lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)

    if not mapv is None:
        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        dx = lon0.shape[1]
        dy = lat0.shape[0]

        nlon = lon0[0,0]
        xlon = lon0[dy-1,dx-1]
        nlat = lat0[0,0]
        xlat = lat0[dy-1,dx-1]

# Thats too much! :)
#        if lonlatLims is not None:
#            print '  ' + fname + ': cutting the domain to plot !!!!'
#            plt.xlim(lonlatLims[0], lonlatLims[2])
#            plt.ylim(lonlatLims[1], lonlatLims[3])
#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]

#            if map_proj == 'cyl':
#                nlon = lonlatLims[0]
#                nlat = lonlatLims[1]
#                xlon = lonlatLims[2]
#                xlat = lonlatLims[3]
#            elif map_proj == 'lcc':
#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
#                nlon =  lonlatLims[0]
#                xlon =  lonlatLims[2]
#                nlat =  lonlatLims[1]
#                xlat =  lonlatLims[3]

        lon2 = lon0[dy/2,dx/2]
        lat2 = lat0[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
            print '    available: cyl, lcc'
            quit(-1)
 
        x,y = m(lon0,lat0)

    else:
        x = lon0
        y = lat0

    vsend = np.zeros((2), dtype=np.float)
# Changing limits of the colors
    if type(vs[0]) != type(np.float(1.)):
        if vs[0] == 'Srange':
            vsend[0] = np.min(varsv)
        elif vs[0][0:11] == 'Saroundmean':
            meanv = np.mean(varsv)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = meanv-minextrm
            vsend[1] = meanv+minextrm
        elif vs[0][0:13] == 'Saroundminmax':
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            vsend[0] = minv
            vsend[1] = maxv
        elif vs[0][0:17] == 'Saroundpercentile':
            medianv = np.median(varsv)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsv, valper)
            maxv = np.percentile(varsv, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = medianv-minextrm
            vsend[1] = medianv+minextrm
        elif vs[0][0:5] == 'Smean':
            meanv = np.mean(varsv)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:7] == 'Smedian':
            medianv = np.median(varsv)
            permedian = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permedian
            maxv = np.max(varsv)*permedian
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:11] == 'Spercentile':
            medianv = np.median(varsv)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsv, valper)
            maxv = np.percentile(varsv, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        else:
            print errormsg
            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
            quit(-1)
        print '    ' + fname + ': modified shadow min,max:',vsend
    else:
        vsend[0] = vs[0]

    if type(vs[0]) != type(np.float(1.)):
        if vs[1] == 'range':
            vsend[1] = np.max(varsv)
    else:
        vsend[1] = vs[1]

    plt.rc('text', usetex=True)

    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
    cbar = plt.colorbar()

    if not mapv is None:
        if colorbar == 'gist_gray':
            m.drawcoastlines(color="red")
        else:
            m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])
        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')
    else:
        plt.xlabel(variables_values(dimn[1])[0].replace('_','\_') + ' (' +           \
          units_lunits(dimxu) + ')')
        plt.ylabel(variables_values(dimn[0])[0].replace('_','\_') + ' (' +           \
          units_lunits(dimyu) + ')')

    txpos = pretty_int(x.min(),x.max(),5)
    typos = pretty_int(y.min(),y.max(),5)
    txlabels = list(txpos)
    for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
    tylabels = list(typos)
    for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])

# set the limits of the plot to the limits of the data

    if searchInlist(revas,'transpose'):
        x0 = y
        y0 = x
        x = x0
        y = y0
#    print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()

    if reva0 == 'flip':
        if reva.split('@')[1] == 'x':
            plt.axis([x.max(), x.min(), y.min(), y.max()])
        elif reva.split('@')[1] == 'y':
            plt.axis([x.min(), x.max(), y.max(), y.min()])
        else:
            plt.axis([x.max(), x.min(), y.max(), y.min()])
    else:
        plt.axis([x.min(), x.max(), y.min(), y.max()])

    if mapv is None:
        plt.xticks(txpos, txlabels)
        plt.yticks(typos, tylabels)

# units labels
    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')

    figname = '2Dfields_shadow'
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    
    output_kind(kfig, figname, ifclose)

    return

#Nvals=50
#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
#for j in range(Nvals):
#    for i in range(Nvals):
#      vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.)

#plot_2D_shadow(vals1, 'var1', np.arange(50)*1., np.arange(50)*1., 'ms-1',            \
#  'm', ['lat','lon'], 'rainbow', [0, Nvals], 'ms-1', 'test var1', 'pdf', 'None',     \
#  None, True)
#quit()

def plot_2D_shadow_time(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,colorbar,vs,uts,   \
  vtit,kfig,reva,taxis,tpos,tlabs,ifclose):
    """ Plotting a 2D field with one of the axes being time
      varsv= 2D values to plot with shading
      vnames= shading variable name for the figure
      dim[x/y]v= values at the axes of x and y
      dim[x/y]u= units at the axes of x and y
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      vs= minmum and maximum values to plot in shadow or:
        'Srange': for full range
        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
        'Saroundminmax@val': for min*val,max*val
        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
          percentile_(100-val)-median)
        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
           percentile_(100-val)-median)
      uts= units of the variable to shadow
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= 
        * 'transpose': reverse the axes (x-->y, y-->x)
        * 'flip'@[x/y]: flip the axis x or y
      taxis= Which is the time-axis 
      tpos= positions of the time ticks
      tlabs= labels of the time ticks
      ifclose= boolean value whether figure should be close (finish) or not
    """
    fname = 'plot_2D_shadow_time'

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_2D_shadow_time.__doc__
        quit()

# Definning ticks labels
    if taxis == 'x':
        txpos = tpos
        txlabels = tlabs
        plxlabel = dimxu
        typos = pretty_int(np.min(dimyv),np.max(dimyv),10)
        tylabels = list(typos)
        for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
        plylabel = variables_values(dimn[0])[0].replace('_','\_') + ' (' +           \
          units_lunits(dimyu) + ')'
    else:
        txpos = pretty_int(np.min(dimxv),np.max(dimxv),10)
        txlabels = list(txpos)
        plxlabel = variables_values(dimn[1])[0].replace('_','\_') + ' (' +           \
          units_lunits(dimxu) + ')'
        typos = tpos
        tylabels = tlabs
        plylabel = dimyu

# Transposing/flipping axis
    if reva.find('|') != 0:
        revas = reva.split('|')
        reva0 = ''
    else:
        revas = [reva]
        reva0 = reva

    for rev in revas:
        if rev[0:4] == 'flip':
            reva0 = 'flip'
            if len(reva.split('@')) != 2:
                 print errormsg
                 print '  ' + fname + ': flip is given', reva, 'but not axis!'
                 quit(-1)
            else:
                 print "  flipping '" + rev.split('@')[1] + "' axis !"

        if rev == 'transpose':
            print '  reversing the axes of the figure (x-->y, y-->x)!!'
# Flipping values of variable
            varsv = np.transpose(varsv)
            dxv = dimyv
            dyv = dimxv
            dimxv = dxv
            dimyv = dyv

    if len(dimxv.shape) == 3:
        dxget='1,2'
    elif len(dimxv.shape) == 2:
        dxget='0,1'
    elif len(dimxv.shape) == 1:
        dxget='0'
    else:
        print errormsg
        print '  ' + fname + ': shape of x-values:',dimxv.shape,'not ready!!'
        quit(-1)

    if len(dimyv.shape) == 3:
        dyget='1,2'
    elif len(dimyv.shape) == 2:
        dyget='0,1'
    elif len(dimyv.shape) == 1:
        dyget='0'
    else:
        print errormsg
        print '  ' + fname + ': shape of y-values:',dimyv.shape,'not ready!!'
        quit(-1)

    x,y = dxdy_lonlat(dimxv,dimyv,dxget,dyget)

    plt.rc('text', usetex=True)

    vsend = np.zeros((2), dtype=np.float)
# Changing limits of the colors
    if type(vs[0]) != type(np.float(1.)):
        if vs[0] == 'Srange':
            vsend[0] = np.min(varsv)
        elif vs[0][0:11] == 'Saroundmean':
            meanv = np.mean(varsv)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = meanv-minextrm
            vsend[1] = meanv+minextrm
        elif vs[0][0:13] == 'Saroundminmax':
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            vsend[0] = minv
            vsend[1] = maxv
        elif vs[0][0:17] == 'Saroundpercentile':
            medianv = np.median(varsv)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsv, valper)
            maxv = np.percentile(varsv, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = medianv-minextrm
            vsend[1] = medianv+minextrm
        elif vs[0][0:5] == 'Smean':
            meanv = np.mean(varsv)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permean
            maxv = np.max(varsv)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:7] == 'Smedian':
            medianv = np.median(varsv)
            permedian = np.float(vs[0].split('@')[1])
            minv = np.min(varsv)*permedian
            maxv = np.max(varsv)*permedian
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:11] == 'Spercentile':
            medianv = np.median(varsv)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsv, valper)
            maxv = np.percentile(varsv, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        else:
            print errormsg
            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
            quit(-1)
        print '    ' + fname + ': modified shadow min,max:',vsend
    else:
        vsend[0] = vs[0]

    if type(vs[0]) != type(np.float(1.)):
        if vs[1] == 'range':
            vsend[1] = np.max(varsv)
    else:
        vsend[1] = vs[1]

    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
    cbar = plt.colorbar()

#    print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()

# set the limits of the plot to the limits of the data
    if reva0 == 'flip':
        if reva.split('@')[1] == 'x':
            plt.axis([x.max(), x.min(), y.min(), y.max()])
        elif reva.split('@')[1] == 'y':
            plt.axis([x.min(), x.max(), y.max(), y.min()])
        else:
            plt.axis([x.max(), x.min(), y.max(), y.min()])
    else:
        plt.axis([x.min(), x.max(), y.min(), y.max()])

    if searchInlist(revas, 'transpose'):
        plt.xticks(typos, tylabels)
        plt.yticks(txpos, txlabels)
        plt.xlabel(plylabel)
        plt.ylabel(plxlabel)
    else:
        plt.xticks(txpos, txlabels)
        plt.yticks(typos, tylabels)
        plt.xlabel(plxlabel)
        plt.ylabel(plylabel)

# units labels
    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')

    figname = '2Dfields_shadow_time'
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    
    output_kind(kfig, figname, ifclose)

    return

def plot_2D_shadow_contour(varsv,varcv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,          \
  colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva,mapv):
    """ Adding labels and other staff to the graph
      varsv= 2D values to plot with shading
      varcv= 2D values to plot with contours
      vnames= variable names for the figure
      dim[x/y]v = values at the axes of x and y
      dim[x/y]u = units at the axes of x and y
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      ckind= contour kind
        'cmap': as it gets from colorbar
        'fixc,[colname]': fixed color [colname], all stright lines
        'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
      clabfmt= format of the labels in the contour plot (None, no labels)
      vs= minmum and maximum values to plot in shadow
      vc= vector with the levels for the contour
      uts= units of the variable [u-shadow, u-contour]
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= 
        * 'transpose': reverse the axes (x-->y, y-->x)
        * 'flip'@[x/y]: flip the axis x or y
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lamvbert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt
    fname = 'plot_2D_shadow_contour'

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_2D_shadow_contour.__doc__
        quit()

    if reva[0:4] == 'flip':
        reva0 = 'flip'
        if len(reva.split('@')) != 2:
             print errormsg
             print '  ' + fname + ': flip is given', reva, 'but not axis!'
             quit(-1)
    else:
        reva0 = reva

    if reva0 == 'transpose':
        print '  reversing the axes of the figure (x-->y, y-->x)!!'
        varsv = np.transpose(varsv)
        varcv = np.transpose(varcv)
        dxv = dimyv
        dyv = dimxv
        dimxv = dxv
        dimyv = dyv

    if not mapv is None:
        if len(dimxv[:].shape) == 3:
            lon0 = dimxv[0,]
            lat0 = dimyv[0,]
        elif len(dimxv[:].shape) == 2:
            lon0 = dimxv[:]
            lat0 = dimyv[:]
        elif len(dimxv[:].shape) == 1:
            lon00 = dimxv[:]
            lat00 = dimyv[:]
            lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
            lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

            for iy in range(len(lat00)):
                lon0[iy,:] = lon00
            for ix in range(len(lon00)):
                lat0[:,ix] = lat00

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        dx = lon0.shape[1]
        dy = lon0.shape[0]

        nlon = lon0[0,0]
        xlon = lon0[dy-1,dx-1]
        nlat = lat0[0,0]
        xlat = lat0[dy-1,dx-1]

# Thats too much! :)
#        if lonlatLims is not None:
#            print '  ' + fname + ': cutting the domain to plot !!!!'
#            plt.xlim(lonlatLims[0], lonlatLims[2])
#            plt.ylim(lonlatLims[1], lonlatLims[3])
#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]

#            if map_proj == 'cyl':
#                nlon = lonlatLims[0]
#                nlat = lonlatLims[1]
#                xlon = lonlatLims[2]
#                xlat = lonlatLims[3]
#            elif map_proj == 'lcc':
#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
#                nlon =  lonlatLims[0]
#                xlon =  lonlatLims[2]
#                nlat =  lonlatLims[1]
#                xlat =  lonlatLims[3]

        lon2 = lon0[dy/2,dx/2]
        lat2 = lat0[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)

        if len(dimxv.shape) == 1:
            lons, lats = np.meshgrid(dimxv, dimyv)
        else:
            if len(dimxv.shape) == 3:
                lons = dimxv[0,:,:]
                lats = dimyv[0,:,:]
            else:
                lons = dimxv[:]
                lats = dimyv[:]
 
        x,y = m(lons,lats)

    else:
        if len(dimxv.shape) == 2:
            x = dimxv
        else:
            if len(dimyv.shape) == 1:
                x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
                for j in range(len(dimyv)):
                    x[j,:] = dimxv
            else:
                x = np.zeros((dimyv.shape), dtype=np.float)
                if x.shape[0] == dimxv.shape[0]:
                    for j in range(x.shape[1]):
                        x[:,j] = dimxv
                else:
                    for j in range(x.shape[0]):
                        x[j,:] = dimxv

        if len(dimyv.shape) == 2:
            y = dimyv
        else:
            if len(dimxv.shape) == 1:
                y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
                for i in range(len(dimxv)):
                    y[:,i] = dimyv
            else:
                y = np.zeros((dimxv.shape), dtype=np.float)

                if y.shape[0] == dimyv.shape[0]:
                    for i in range(y.shape[1]):
                        y[i,:] = dimyv
                else:
                    for i in range(y.shape[0]):
                        y[i,:] = dimyv

    plt.rc('text', usetex=True)

    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
    cbar = plt.colorbar()

# contour
##
    contkind = ckind.split(',')[0]
    if contkind == 'cmap':
        cplot = plt.contour(x, y, varcv, levels=vc)
    elif  contkind == 'fixc':
        plt.rcParams['contour.negative_linestyle'] = 'solid'
        coln = ckind.split(',')[1]
        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
    elif  contkind == 'fixsigc':
        coln = ckind.split(',')[1]
        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
    else:
        print errormsg
        print '  ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
        quit(-1)

    if clabfmt is not None:
        plt.clabel(cplot, fmt=clabfmt)
        mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
        maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
    else:
        mincntS = '{:g}'.format(vc[0])
        maxcntS = '{:g}'.format(vc[len(vc)-1])        

    if not mapv is None:
        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])
        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')
    else:
        plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
        plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')

        txpos = pretty_int(x.min(),x.max(),5)
        typos = pretty_int(y.min(),y.max(),5)
        txlabels = list(txpos)
        for i in range(len(txlabels)): txlabels[i] = '{:.1f}'.format(txlabels[i])
        tylabels = list(typos)
        for i in range(len(tylabels)): tylabels[i] = '{:.1f}'.format(tylabels[i])
        plt.xticks(txpos, txlabels)
        plt.yticks(typos, tylabels)

# set the limits of the plot to the limits of the data
    if reva0 == 'flip':
        if reva.split('@')[1] == 'x':
            plt.axis([x.max(), x.min(), y.min(), y.max()])
        else:
            plt.axis([x.min(), x.max(), y.max(), y.min()])
    else:
        plt.axis([x.min(), x.max(), y.min(), y.max()])


# units labels
    cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
    plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' +  \
      mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction',    \
      color=coln)

    figname = '2Dfields_shadow-contour'
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    
    output_kind(kfig, figname, True)

    return

#Nvals=50
#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
#vals2 = np.zeros((Nvals,Nvals), dtype= np.float)
#for j in range(Nvals):
#    for i in range(Nvals):
#      vals1[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.)
#      vals2[j,i]=np.sqrt((j-Nvals/2)**2. + (i-Nvals/2)**2.) - Nvals/2

#prettylev=pretty_int(-Nvals/2,Nvals/2,10)

#plot_2D_shadow_contour(vals1, vals2, ['var1', 'var2'], np.arange(50)*1.,             \
#  np.arange(50)*1., ['x-axis','y-axis'], 'rainbow', 'fixc,b', "%.2f", [0, Nvals],    \
#  prettylev, ['$ms^{-1}$','$kJm^{-1}s^{-1}$'], 'test var1 & var2', 'pdf', False)

def plot_2D_shadow_contour_time(varsv,varcv,vnames,valv,timv,timpos,timlab,valu,     \
  timeu,axist,dimn,colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva):
    """ Adding labels and other staff to the graph
      varsv= 2D values to plot with shading
      varcv= 2D values to plot with contours
      vnames= variable names for the figure
      valv = values at the axes which is not time
      timv = values for the axis time
      timpos = positions at the axis time
      timlab = labes at the axis time
      valu = units at the axes which is not time
      timeu = units at the axes which is time
      axist = which is the axis time
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      ckind= contour kind
        'cmap': as it gets from colorbar
        'fixc,[colname]': fixed color [colname], all stright lines
        'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
      clabfmt= format of the labels in the contour plot (None, no labels)
      vs= minmum and maximum values to plot in shadow
      vc= vector with the levels for the contour
      uts= units of the variable [u-shadow, u-contour]
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= 
        * 'transpose': reverse the axes (x-->y, y-->x)
        * 'flip'@[x/y]: flip the axis x or y
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt
    fname = 'plot_2D_shadow_contour_time'

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_2D_shadow_contour_time.__doc__
        quit()

    # Usually axis > x must be the time, thus...
    dimxv0 = valv.copy()
    dimyv0 = timv.copy()
    dimxt0 = np.array(pretty_int(valv.min(),valv.max(),10))
    dimyt0 = np.array(timpos)
    dimxl0 = []
    for i in range(len(dimxt0)): dimxl0.append(str(dimxt0[i]))
    dimyl0 = list(timlab)
    dimxT0 = variables_values(dimn[0])[0] + ' (' + units_lunits(valu) + ')'
    dimyT0 = timeu

    if axist == 'x':
        reva = 'transpose'

    if reva is not None:
        varcv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT =              \
          transform(varcv, reva)
        varsv, dimxv, dimyv, dimxt, dimyt, dimxl, dimyl, dimxT, dimyT =              \
          transform(varsv, reva, dxv=dimxv0, dyv=dimyv0, dxt=dimxt0, dyt=dimyt0,     \
          dxl=dimxl0, dyl=dimyl0, dxtit=dimxT0, dytit=dimyT0)
    else:
        dimxv = dimxv0
        dimyv = dimyv0
        dimxt = dimxt0
        dimyt = dimyt0
        dimxl = dimxl0
        dimyl = dimyl0
        dimxT = dimxT0
        dimyT = dimyT0

    dx=varsv.shape[1]
    dy=varsv.shape[0]

    x,y = gen.lonlat2D(dimxv,dimyv)

    plt.rc('text', usetex=True)

    plt.xticks(dimxt, dimxl)
    plt.yticks(dimyt, dimyl)

    plt.xlim(0,dx)
    plt.ylim(0,dy)

    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
    cbar = plt.colorbar()
    
# contour
##
    contkind = ckind.split(',')[0]
    if contkind == 'cmap':
        cplot = plt.contour(x, y, varcv, levels=vc)
    elif  contkind == 'fixc':
        plt.rcParams['contour.negative_linestyle'] = 'solid'
        coln = ckind.split(',')[1]
        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
    elif  contkind == 'fixsigc':
        coln = ckind.split(',')[1]
        cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
    else:
        print errormsg
        print '  ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
        quit(-1)

    if clabfmt is not None:
        plt.clabel(cplot, fmt=clabfmt)
        mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
        maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
    else:
        mincntS = '{:g}'.format(vc[0])
        maxcntS = '{:g}'.format(vc[len(vc)-1])        

    plt.xlabel(dimxT)
    plt.ylabel(dimyT)

# set the limits of the plot to the limits of the data
    plt.axis([dimxv.min(), dimxv.max(), dimyv.min(), dimyv.max()])

# units labels
    cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
    plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' +  \
      mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction',    \
      color=coln)

    figname = '2Dfields_shadow-contour'
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    
    output_kind(kfig, figname, True)

    return

def dxdy_lonlat(dxv,dyv,ddx,ddy):
    """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
    dxdy_lonlat(dxv,dyv,Lv,lv)
      dx: values for the x
      dy: values for the y
      ddx: ',' list of which dimensions to use from values along x
      ddy: ',' list of which dimensions to use from values along y
    """

    fname = 'dxdy_lonlat'

    if ddx.find(',') > -1:
        dxk = 2
        ddxv = ddx.split(',')
        ddxy = int(ddxv[0])
        ddxx = int(ddxv[1])
    else:
        dxk = 1
        ddxy = int(ddx)
        ddxx = int(ddx)

    if ddy.find(',') > -1:
        dyk = 2
        ddyv = ddy.split(',')
        ddyy = int(ddyv[0])
        ddyx = int(ddyv[1])
    else:
        dyk = 1
        ddyy = int(ddy)
        ddyx = int(ddy)

    ddxxv = dxv.shape[ddxx]
    ddxyv = dxv.shape[ddxy]
    ddyxv = dyv.shape[ddyx]
    ddyyv = dyv.shape[ddyy]

    slicex = []
    if len(dxv.shape) > 1:
        for idim in range(len(dxv.shape)):
            if idim == ddxx or idim == ddxy:
                slicex.append(slice(0,dxv.shape[idim]))
            else:
                slicex.append(0)
    else:
        slicex.append(slice(0,len(dxv)))

    slicey = []
    if len(dyv.shape) > 1:
        for idim in range(len(dyv.shape)):
            if idim == ddyx or idim == ddyy:
                slicey.append(slice(0,dyv.shape[idim]))
            else:
                slicey.append(0)
    else:
        slicey.append(slice(0,len(dyv)))

    if dxk == 2 and dyk == 2:
        if ddxxv != ddyxv:
            print errormsg
            print '  ' + fname + ': wrong dx dimensions! ddxx=',ddxxv,'ddyx=',ddyxv
            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
            quit(-1)
        if ddxyv != ddyyv:
            print errormsg
            print '  ' + fname + ': wrong dy dimensions! ddxy=',ddxyv,'ddyy=',ddyv
            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
            quit(-1)
        dx = ddxxv
        dy = ddxyv

        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
        lonv = np.zeros((dy,dx), dtype=np.float)
        latv = np.zeros((dy,dx), dtype=np.float)


        lonv = dxv[tuple(slicex)]
        latv = dyv[tuple(slicey)]

    elif dxk == 2 and dyk == 1:
        if not ddxxv == ddyxv and not ddxyv == ddyyv:
            print errormsg
            print '  ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv,    \
              'ddyx=',ddyxv,'ddyy=',ddyyv
            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
            quit(-1)
        dx = ddxvv
        dy = ddxyv

        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
        lonv = np.zeros((dy,dx), dtype=np.float)
        latv = np.zeros((dy,dx), dtype=np.float)
        lonv = dxv[tuple(slicex)]

        if ddxxv == ddyxv: 
            for iy in range(dy):
                latv[iy,:] = dyv[tuple(slicey)]
        else:
            for ix in range(dx):
                latv[:,ix] = dyv[tuple(slicey)]

    elif dxk == 1 and dyk == 2:
        if not ddxxv == ddyxv and not ddxyv == ddyyv:
            print errormsg
            print '  ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv,    \
              'ddyx=',ddyxv,'ddyy=',ddyyv
            print '    choose another for x:',dxv.shape,'or y:',dyv.shape
            quit(-1)
        dx = ddyxv
        dy = ddyyv
 
        print '  ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
        lonv = np.zeros((dy,dx), dtype=np.float)
        latv = np.zeros((dy,dx), dtype=np.float)

        latv = dyv[tuple(slicey)]

        if ddyxv == ddxxv: 
            for iy in range(dy):
                lonv[iy,:] = dxv[tuple(slicex)]
        else:
            for ix in range(dx):
                lonv[:,ix] = dxv[tuple(slicex)]


    elif dxk == 1 and dyk == 1:
        dx = ddxxv
        dy = ddyyv
 
#        print 'dx:',dx,'dy:',dy 

        lonv = np.zeros((dy,dx), dtype=np.float)
        latv = np.zeros((dy,dx), dtype=np.float)

        for iy in range(dy):
            lonv[iy,:] = dxv[tuple(slicex)]
        for ix in range(dx):
            latv[:,ix] = dyv[tuple(slicey)]

    return lonv,latv

def dxdy_lonlatDIMS(dxv,dyv,dnx,dny,dd):
    """ Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values for a given
      list of values
    dxdy_lonlat(dxv,dyv,Lv,lv)
      dxv: values for the x
      dyv: values for the y
      dnx: mnames of the dimensions for values on x
      dny: mnames of the dimensions for values on y
      dd: list of [dimname]|[val] for the dimensions use
        [dimname]: name of the dimension
        [val]: value (-1 for all the range)
    """
    fname = 'dxdy_lonlatDIMS'

    slicex = []
    ipos=0
    for dn in dnx:
        for idd in range(len(dd)):
            dname = dd[idd].split('|')[0]
            dvalue = dd[idd].split('|')[1]
            if dn == dname:
                if dvalue.find('@') != -1:
                    slicex.append(slice(int(dvalue.split('@')[0]),                   \
                      int(dvalue.split('@')[1])))
                else:
                    if int(dvalue) == -1:
                        slicex.append(slice(0,dxv.shape[ipos]))
                    elif int(dvalue) == -9:
                        slicex.append(dxv.shape[ipos]-1)
                    else:
                        slicex.append(int(dvalue))
                    break
        ipos = ipos + 1

    slicey = []
    ipos=0
    for dn in dny:
        for idd in range(len(dd)):
            dname = dd[idd].split('|')[0]
            dvalue = dd[idd].split('|')[1]
            if dn == dname:
                if dvalue.find('@') != -1:
                    slicey.append(slice(int(dvalue.split('@')[0]),                   \
                      int(dvalue.split('@')[1])))
                else:
                    if int(dvalue) == -1:
                        slicey.append(slice(0,dyv.shape[ipos]))
                    elif int(dvalue) == -9:
                        slicey.append(dyv.shape[ipos]-1)
                    else:
                        slicey.append(int(dvalue))
                    break
        ipos = ipos + 1

    lonv = dxv[tuple(slicex)]
    latv = dyv[tuple(slicey)]

    if len(lonv.shape) != len(latv.shape):
        print '  ' + fname + ': dimension size on x:', len(lonv.shape), 'and on y:', \
          len(latv.shape),'do not coincide!!'
        quit(-1)

    return lonv,latv

def plot_2D_shadow_line(varsv,varlv,vnames,vnamel,dimxv,dimyv,dimxu,dimyu,dimn,             \
  colorbar,colln,vs,uts,utl,vtit,kfig,reva,mapv,ifclose):
    """ Plotting a 2D field with shadows and another one with a line
      varsv= 2D values to plot with shading
      varlv= 1D values to plot with line
      vnames= variable names for the shadow variable in the figure
      vnamel= variable names for the line varibale in the figure
      dim[x/y]v = values at the axes of x and y
      dim[x/y]u = units at the axes of x and y
      dimn= dimension names to plot
      colorbar= name of the color bar to use
      colln= color for the line
      vs= minmum and maximum values to plot in shadow
      uts= units of the variable to shadow
      utl= units of the variable to line
      vtit= title of the variable
      kfig= kind of figure (jpg, pdf, png)
      reva= 
        * 'transpose': reverse the axes (x-->y, y-->x)
        * 'flip'@[x/y]: flip the axis x or y
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      ifclose= boolean value whether figure should be close (finish) or not
    """
##    import matplotlib as mpl
##    mpl.use('Agg')
##    import matplotlib.pyplot as plt
    fname = 'plot_2D_shadow_line'

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_2D_shadow_line.__doc__
        quit()

    if reva[0:4] == 'flip':
        reva0 = 'flip'
        if len(reva.split('@')) != 2:
             print errormsg
             print '  ' + fname + ': flip is given', reva, 'but not axis!'
             quit(-1)
    else:
        reva0 = reva

    if reva0 == 'transpose':
        print '  reversing the axes of the figure (x-->y, y-->x)!!'
        varsv = np.transpose(varsv)
        dxv = dimyv
        dyv = dimxv
        dimxv = dxv
        dimyv = dyv

    if len(dimxv[:].shape) == 3:
        lon0 = dimxv[0,]
    elif len(dimxv[:].shape) == 2:
        lon0 = dimxv[:]

    if len(dimyv[:].shape) == 3:
        lat0 = dimyv[0,]
    elif len(dimyv[:].shape) == 2:
        lat0 = dimyv[:]

    if len(dimxv[:].shape) == 1 and len(dimyv[:].shape) == 1:
        lon00 = dimxv[:]
        lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )

        for iy in range(len(lat00)):
            lon0[iy,:] = lon00
        for ix in range(len(lon00)):
            lat0[:,ix] = lat00

    if not mapv is None:
        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        dx = lon0.shape[1]
        dy = lat0.shape[0]

        nlon = lon0[0,0]
        xlon = lon0[dy-1,dx-1]
        nlat = lat0[0,0]
        xlat = lat0[dy-1,dx-1]

# Thats too much! :)
#        if lonlatLims is not None:
#            print '  ' + fname + ': cutting the domain to plot !!!!'
#            plt.xlim(lonlatLims[0], lonlatLims[2])
#            plt.ylim(lonlatLims[1], lonlatLims[3])
#            print '    limits: W-E', lonlatLims[0], lonlatLims[2]
#            print '    limits: N-S', lonlatLims[1], lonlatLims[3]

#            if map_proj == 'cyl':
#                nlon = lonlatLims[0]
#                nlat = lonlatLims[1]
#                xlon = lonlatLims[2]
#                xlat = lonlatLims[3]
#            elif map_proj == 'lcc':
#                lon2 = (lonlatLims[0] + lonlatLims[2])/2.
#                lat2 = (lonlatLims[1] + lonlatLims[3])/2.
#                nlon =  lonlatLims[0]
#                xlon =  lonlatLims[2]
#                nlat =  lonlatLims[1]
#                xlat =  lonlatLims[3]

        lon2 = lon0[dy/2,dx/2]
        lat2 = lat0[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": map projection '" + map_proj + "' not defined!!!"
            print '    available: cyl, lcc'
            quit(-1)

        if len(dimxv.shape) == 1:
            lons, lats = np.meshgrid(dimxv, dimyv)
        else:
            if len(dimxv.shape) == 3:
                lons = dimxv[0,:,:]
            else:
                lons = dimxv[:]

            if len(dimyv.shape) == 3:
                lats = dimyv[0,:,:]
            else:
                lats = dimyv[:]
 
        x,y = m(lons,lats)

    else:
        if len(dimxv.shape) == 3:
            x = dimxv[0,:,:]
        elif len(dimxv.shape) == 2:
            x = dimxv
        else:
# Attempt of simplier way...
#            x = np.zeros((lon0.shape), dtype=np.float)
#            for j in range(lon0.shape[0]):
#                x[j,:] = dimxv

## This way is too complicated and maybe not necessary ? (assuming dimxv.shape == dimyv.shape)
            if len(dimyv.shape) == 1:
                x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
                for j in range(len(dimxv)):
                    x[j,:] = dimxv
            else:
                x = np.zeros((dimyv.shape), dtype=np.float)
                if x.shape[0] == dimxv.shape[0]:
                    for j in range(x.shape[1]):
                        x[:,j] = dimxv
                else:
                    for j in range(x.shape[0]):
                        x[j,:] = dimxv

        if len(dimyv.shape) == 3:
            y = dimyv[0,:,:]
        elif len(dimyv.shape) == 2:
            y = dimyv
        else:
#            y = np.zeros((lat0.shape), dtype=np.float)
#            for i in range(lat0.shape[1]):
#                x[:,i] = dimyv

# Idem
            if len(dimxv.shape) == 1:
                y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
                for i in range(len(dimxv)):
                    y[:,i] = dimyv
            else:
                y = np.zeros((dimxv.shape), dtype=np.float)
                if y.shape[0] == dimyv.shape[0]:
                    for i in range(y.shape[1]):
                        y[:,i] = dimyv
                else:
                    for j in range(y.shape[0]):
                        y[j,:] = dimyv

    plt.rc('text', usetex=True)

    plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
    cbar = plt.colorbar()

    if not mapv is None:
        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])
        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])

        plt.xlabel('W-E')
        plt.ylabel('S-N')
    else:
        plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
        plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')

# Line
##

    if reva0 == 'flip' and reva.split('@')[1] == 'y':
        b=-np.max(y[0,:])/np.max(varlv)
        a=np.max(y[0,:])
    else:
        b=np.max(y[0,:])/np.max(varlv)
        a=0.

    newlinv = varlv*b+a
    if reva0 == 'transpose':
        plt.plot(newlinv, x[0,:], '-', color=colln, linewidth=2)
    else:
        plt.plot(x[0,:], newlinv, '-', color=colln, linewidth=2)

    txpos = pretty_int(x.min(),x.max(),10)
    typos = pretty_int(y.min(),y.max(),10)
    txlabels = list(txpos)
    for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
    tylabels = list(typos)
    for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])

    tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
    for it in range(len(tllabels)):
        yval = (tllabels[it]*b+a)
        plt.plot([x.max()*0.97, x.max()], [yval, yval], '-', color='k')
        plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)),             \
          xycoords='axes fraction')

# set the limits of the plot to the limits of the data
    if reva0 == 'flip':
        if reva.split('@')[1] == 'x':
            plt.axis([x.max(), x.min(), y.min(), y.max()])
        else:
            plt.axis([x.min(), x.max(), y.max(), y.min()])
    else:
        plt.axis([x.min(), x.max(), y.min(), y.max()])

    plt.tick_params(axis='y',right='off')
    if mapv is None:
        plt.xticks(txpos, txlabels)
        plt.yticks(typos, tylabels)

    tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
    for it in range(len(tllabels)):
        plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)), xycoords='axes fraction')

# units labels
    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')

    plt.annotate(vnamel +' (' + units_lunits(utl) + ')', xy=(0.75,0.04), 
      xycoords='figure fraction', color=colln)
    figname = '2Dfields_shadow_line'
    graphtit = vtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    
    output_kind(kfig, figname, ifclose)

    return

def plot_Neighbourghood_evol(varsv, dxv, dyv, vnames, ttits, tpos, tlabels, colorbar, \
  Nng, vs, uts, gtit, kfig, ifclose):
    """ Plotting neighbourghood evolution
      varsv= 2D values to plot with shading
      vnames= shading variable name for the figure
      d[x/y]v= values at the axes of x and y
      ttits= titles of both time axis
      tpos= positions of the time ticks
      tlabels= labels of the time ticks
      colorbar= name of the color bar to use
      Nng= Number of grid points of the full side of the box (odd value)
      vs= minmum and maximum values to plot in shadow or:
        'Srange': for full range
        'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
        'Saroundminmax@val': for min*val,max*val
        'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
          percentile_(100-val)-median)
        'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
        'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
        'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
           percentile_(100-val)-median)
      uts= units of the variable to shadow
      gtit= title of the graph
      kfig= kind of figure (jpg, pdf, png)
      ifclose= boolean value whether figure should be close (finish) or not
    """
    import numpy.ma as ma

    fname = 'plot_Neighbourghood_evol'

    if varsv == 'h':
        print fname + '_____________________________________________________________'
        print plot_Neighbourghood_evol.__doc__
        quit()

    if len(varsv.shape) != 2:
        print errormsg
        print '  ' + fname + ': wrong number of dimensions of the values: ',         \
          varsv.shape
        quit(-1)

    varsvmask = ma.masked_equal(varsv,fillValue)

    vsend = np.zeros((2), dtype=np.float)
# Changing limits of the colors
    if type(vs[0]) != type(np.float(1.)):
        if vs[0] == 'Srange':
            vsend[0] = np.min(varsvmask)
        elif vs[0][0:11] == 'Saroundmean':
            meanv = np.mean(varsvmask)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsvmask)*permean
            maxv = np.max(varsvmask)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = meanv-minextrm
            vsend[1] = meanv+minextrm
        elif vs[0][0:13] == 'Saroundminmax':
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsvmask)*permean
            maxv = np.max(varsvmask)*permean
            vsend[0] = minv
            vsend[1] = maxv
        elif vs[0][0:17] == 'Saroundpercentile':
            medianv = np.median(varsvmask)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsvmask, valper)
            maxv = np.percentile(varsvmask, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = medianv-minextrm
            vsend[1] = medianv+minextrm
        elif vs[0][0:5] == 'Smean':
            meanv = np.mean(varsvmask)
            permean = np.float(vs[0].split('@')[1])
            minv = np.min(varsvmask)*permean
            maxv = np.max(varsvmask)*permean
            minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:7] == 'Smedian':
            medianv = np.median(varsvmask)
            permedian = np.float(vs[0].split('@')[1])
            minv = np.min(varsvmask)*permedian
            maxv = np.max(varsvmask)*permedian
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        elif vs[0][0:11] == 'Spercentile':
            medianv = np.median(varsvmask)
            valper = np.float(vs[0].split('@')[1])
            minv = np.percentile(varsvmask, valper)
            maxv = np.percentile(varsvmask, 100.-valper)
            minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
            vsend[0] = -minextrm
            vsend[1] = minextrm
        else:
            print errormsg
            print '  ' + fname + ": range '" + vs[0] + "' not ready!!!"
            quit(-1)
        print '    ' + fname + ': modified shadow min,max:',vsend
    else:
        vsend[0] = vs[0]

    if type(vs[0]) != type(np.float(1.)):
        if vs[1] == 'range':
            vsend[1] = np.max(varsv)
    else:
        vsend[1] = vs[1]

    plt.rc('text', usetex=True)

#    plt.pcolormesh(dxv, dyv, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
    plt.pcolormesh(varsvmask, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
    cbar = plt.colorbar()

    newtposx = (tpos[0][:] - np.min(dxv)) * len(dxv) * Nng / (np.max(dxv) - np.min(dxv))
    newtposy = (tpos[1][:] - np.min(dyv)) * len(dyv) * Nng / (np.max(dyv) - np.min(dyv))

    plt.xticks(newtposx, tlabels[0])
    plt.yticks(newtposy, tlabels[1])
    plt.xlabel(ttits[0])
    plt.ylabel(ttits[1])

    plt.axes().set_aspect('equal')
# From: http://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib
    plt.axes().xaxis.tick_top
    plt.axes().xaxis.set_ticks_position('top')

# units labels
    cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')

    figname = 'Neighbourghood_evol'
    graphtit = gtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit, position=(0.5,1.05))
    
    output_kind(kfig, figname, ifclose)

    return

def plot_lines(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, gtit, gloc, kfig):
    """ Function to plot a collection of lines
      vardv= list of set of dimension values
      varvv= list of set of values
      vaxis= which axis will be used for the values ('x', or 'y')
      dtit= title for the common dimension
      linesn= names for the legend
      vtit= title for the vaxis
      vunit= units of the vaxis
      gtit= main title
      gloc= location of the legend (0, autmoatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'
      kfig= kind of figure
      plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)',      \
  ['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
    """
    fname = 'plot_lines'

    if vardv == 'h':
        print fname + '_____________________________________________________________'
        print plot_lines.__doc__
        quit()

# Canging line kinds every 7 lines (end of standard colors)
    linekinds=['.-','x-','o-']
    pointkindsauto=['.', ',', 'x', 'o', '*', '+', '<', '|', '_', '>', '1', '8', 's', \
      'p', 'h', 'D']

    Ntraj = len(vardv)

    N7lines = 0

    xmin = 100000.
    xmax = -100000.
    ymin = 100000.
    ymax = -100000.
    for il in range(Ntraj):
        minv = np.min(varvv[il])
        maxv = np.max(varvv[il])
        mind = np.min(vardv[il])
        maxd = np.max(vardv[il])

        if minv < xmin: xmin = minv
        if maxv > xmax: xmax = maxv
        if mind < ymin: ymin = mind
        if maxd > ymax: ymax = maxd

    print 'x:',xmin,',',xmax,'y:',ymin,ymax

    plt.rc('text', usetex=True)

    if vaxis == 'x':
        for il in range(Ntraj):
            plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
            if il == 6: N7lines = N7lines + 1

        plt.xlabel(vtit + ' (' + vunit + ')')
        plt.ylabel(dtit)
        plt.xlim(xmin,xmax)
        plt.ylim(ymin,ymax)

    else:
        for il in range(Ntraj):
            plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
            if il == 6: N7lines = N7lines + 1

        plt.xlabel(dtit)
        plt.ylabel(vtit + ' (' + vunit + ')')
        
        plt.xlim(ymin,ymax)
        plt.ylim(xmin,xmax)

    figname = 'lines'
    graphtit = gtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    plt.legend(loc=gloc)
    
    output_kind(kfig, figname, True)

    return

def ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq):
    """ Function to provide the colors, lines & points styles
      Nstyles= total number of styles
      colors= list of colors, None for no value
      lines= list of lines, None for no value
      points= list of points, None for no value
      lwidths= list of line widths, None for no value
      psizes= list of point sizes, None for no value
      ptfreq= frequency for the points, None for all values
    >>> Nstyles = 3
    >>> colors = ['blue']
    >>> lines = ['-']
    >>> points = ['x', '*', 'o']
    >>> lwidths = ['1']
    >>> psizes = ['2']
    >>> ptfreq = 2
    >>> colors, lines, points = ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq)
    ['blue', 'blue', 'blue']
    ['-', '-', '-']
    ['x', '*', 'o']
    ['2', '2', '2']
    ['1', '1', '1']
    """
    fname = 'ColorsLinesPointsStyles'

    usecolors = []
    uselines = []
    usepoints = []
    usewlines = []
    usespoints = []

# Colors
    if colors is None:
        Ncols = len(colorsauto)
        for ic in range(Nstyles):
            iic = np.mod(ic,Ncols) - 1
            if iic == 0: iic = Ncols - 1
            usecolors.append(colorsauto[iic])
    else:
        Ncols = len(colors)
        if Ncols == 1:
            for ic in range(Nstyles):
                usecolors.append(colors[0])
        else:
            if Ncols != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided colors:', Ncols,          \
                  'and required:', Nstyles,'differ !!'
                quit(-1)
            usecolors = colors

# Lines
    if lines is None:
        Nklns = len(linekindsauto)
        for il in range(Nstyles):
            iil = np.mod(il,Nklns) - 1
            if iil == 0: iil = Nklns - 1
            uselines.append(linekindsauto[iil])
    else:
        Nklns = len(lines)
        if Nklns == 1:
            for il in range(Nstyles):
                uselines.append(lines[0])
        else:
            if Nklns != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided lines:', Nklns,           \
                  'and required:', Nstyles,'differ !!'
                quit(-1)
            uselines = lines

# Points
    if points is None:
        Nkpts = len(pointkindsauto)
        for ip in range(Nstyles):
            iip = np.mod(ip,Nkpts) - 1
            if iip == 0: iip = Nkpts - 1
            usepoints.append(pointkindsauto[iip])
    else:
        Nkpts = len(points)
        if Nkpts == 1:
            for ip in range(Nstyles):
                usepoints.append(points[0])
        else:
            if Nkpts != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided points:', Nkpts,          \
                  'and required:', Nstyles,'differ !!'
                quit(-1)
            usepoints = points

# Line widths
    if lwidths is None:
        Nwlns = len(linewidthsauto)
        for il in range(Nstyles):
            iil = np.mod(il,Nwlns) - 1
            if iil == 0: iil = Nwlns - 1
            usewlines.append(linewidthsauto[iil])
    else:
        Nwlns = len(lwidths)
        if Nwlns == 1:
            for il in range(Nstyles):
                usewlines.append(lwidths[0])
        else:
            if Nwlns != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided line widthss:', Nwlns,    \
                  'and required:', Nstyles, 'differ !!'
                quit(-1)
            usewlines = lwidths

# Point sizes
    if psizes is None:
        Nspts = len(pointsizesauto)
        for ip in range(Nstyles):
            iip = np.mod(ip,Nspts) - 1
            if iip == 0: iip = Nspts - 1
            usespoints.append(pointsizesauto[iip])
    else:
        Nspts = len(psizes)
        if Nspts == 1:
            for ip in range(Nstyles):
                usespoints.append(psizes[0])
        else:
            if Nspts != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided point sizes:', Nspts,     \
                  'and required:', Nstyles, 'differ !!'
                quit(-1)
            usespoints = psizes

    Ncols = len(usecolors)
    Nlins = len(uselines)
    Npnts = len(usepoints)

    lkinds = []
    pkinds = []

    lcolors = usecolors

# Old way
#    if ptfreq is not None:
#        lkinds = uselines
#        pkinds = usepoints
#    else:
#        pkinds = usepoints
#        lkinds = uselines
#        for ilp in range(Nstyles):
#            lkinds.append(usepoints[ilp] + uselines[ilp])
            
    return lcolors, uselines, usepoints, usewlines, usespoints

#Nstyles = 3
#colors = ['blue']
#lines = ['-']
#points = ['x', '*', 'o']
#lwidths = ['2']
#psizes = ['1']
#ptfreq = 2
#print ColorsLinesPointsStyles(Nstyles, colors, lines, points, lwidths, psizes, ptfreq)

def LinesPointsStyles(Nstyles, lines, points, ptfreq):
    """ Function to provide the lines & points styles
      Nstyles= total number of styles
      lines= list of lines, None for no value
      points= list of points, None for no value
      ptfreq= frequency for the points, None for all values
    >>> Nstyles = 3
    >>> lines = ['-']
    >>> points = ['x', '*', 'o']
    >>> ptfreq = 2
    >>> lines, points = LinesPointsStyles(Nstyles, lines, points, ptfreq)
    ['-', '-', '-']
    ['x', '*', 'o']
    """
    fname = 'LinesPointsStyles'

# Canging line kinds every 7 lines (end of standard colors)
    uselines = []
    usepoints = []
    if lines is None:
        Nklns = len(linekindsauto)
        for il in range(Nstyles):
            iil = np.mod(il,Nklns)
            uselines.append(linekindsauto[iil])
    else:
        Nklns = len(lines)
        if Nklns == 1:
            for il in range(Nstyles):
                uselines.append(lines[0])
        else:
            if Nklns != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided lines:', Nklns,           \
                  'and required:', Nstyles,'differ !!'
                quit(-1)
            uselines = lines

    if points is None:
        Nkpts = len(pointkindsauto)
        for ip in range(Nstyles):
            iip = np.mod(ip,Nkpts)
            usepoints.append(pointkindsauto[iip])
    else:
        Nkpts = len(points)
        if Nkpts == 1:
            for ip in range(Nstyles):
                usepoints.append(points[0])
        else:
            if Nkpts != Nstyles:
                print errormsg
                print '  ' + fname + ': number of provided points:', Nkpts,          \
                  'and required:', Nstyles,'differ !!'
                quit(-1)
            usepoints = points

    Nlins = len(uselines)
    Npnts = len(usepoints)

    lkinds = []
    pkinds = []
    if ptfreq is not None:
        lkinds = uselines
        pkinds = usepoints
    else:
        pkinds = usepoints
        for ilp in range(Nstyles):
            lkinds.append(usepoints[ilp] + uselines[ilp])
            
    return lkinds, pkinds

def plot_lines_time(vardv, varvv, vaxis, dtit, linesn0, vtit, vunit, tpos, tlabs,    \
  gtit, gloc, gsize, kfig, vmin, vmax, lsl, coll, ptl, lwidth, psize, ptf):
    """ Function to plot a collection of lines with a time axis
      vardv= list of set of dimension values
      varvv= list of set of values
      vaxis= which axis will be used for the time values ('x', or 'y')
      dtit= title for the common dimension
      linesn= names for the legend (None, no legend)
      vtit= title for the vaxis
      vunit= units of the vaxis
      tpos= positions of the time ticks
      tlabs= labels of the time ticks
      gtit= main title
      vmin= minimum value for the axis-value 'auto' for use the extreme
      vmax= maximum value for the axis-value 'auto' for use the extreme
      gloc= location of the legend (0, autmoatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'
      gsize= fontsize of the legend
      kfig= kind of figure
      lsl= ',' list of line styles
      coll= ',' list of colors for the lines, None for automatic, single
          value all the same
      ptl= ',' list of type of points for the lines, None for automatic, single
          value all the same
      lwidth= ',' list of line widths
      psize= ',' list of point sizes
      ptf= frequency of point plotting, 'all' for all time steps

      plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)',      \
  ['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
    """
    fname = 'plot_lines_time'

    if vardv == 'h':
        print fname + '_____________________________________________________________'
        print plot_lines.__doc__
        quit()

    Ntraj = len(vardv)

    colvalues, linekinds, pointkinds, lwidths, psizes = ColorsLinesPointsStyles(     \
      Ntraj, coll, lsl, ptl, lwidth, psize,  ptf)

    plt.rc('text', usetex=True)
    xtrmvv = [fillValueF,-fillValueF]
    xtrmdv = [fillValueF,-fillValueF]

# Do we have legend?
##
    if linesn0 is None:
        linesn = []
        for itrj in range(Ntraj):
            linesn.append(str(itrj))
    else:
        linesn = linesn0

    if vaxis == 'x':
        for il in range(Ntraj):
            plt.plot(varvv[il], vardv[il], linekinds[il], marker=pointkinds[il],     \
              label=linesn[il], color=colvalues[il], linewidth=lwidths[il],          \
              markersize=psizes[il], markevery=ptf)

# Old way
#            if ptf is None:
#                plt.plot(varvv[il], vardv[il], linekinds[il], label= linesn[il],     \
#                  color=colvalues[il], linewidth=lwidths[il], markersize=psizes[il])
#            else:
#                plt.plot(varvv[il], vardv[il], linekinds[il], color=colvalues[il],   \
#                  linewidth=lwidths[il])
#                plt.plot(varvv[il][::ptf], vardv[il][::ptf], pointkinds[il],         \
#                  label= linesn[il], color=colvalues[il], linewidth=lwidths[il],     \
#                  markersize=psizes[il])

            minvv = np.min(varvv[il])
            maxvv = np.max(varvv[il])
            mindv = np.min(vardv[il])
            maxdv = np.max(vardv[il])

            if minvv < xtrmvv[0]: xtrmvv[0] = minvv
            if maxvv > xtrmvv[1]: xtrmvv[1] = maxvv
            if mindv < xtrmdv[0]: xtrmdv[0] = mindv
            if maxdv > xtrmdv[1]: xtrmdv[1] = maxdv

        plt.xlabel(vtit + ' (' + vunit + ')')
        plt.ylabel(dtit)
#        plt.xlim(np.min(varTvv),np.max(varTvv))
#        plt.ylim(np.min(varTdv),np.max(varTdv))

        if vmin == 'auto':
            xmin = xtrmvv[0]
        else:
            xmin = vmin
        if vmax == 'auto':
            xmax = xtrmvv[1]
        else:
            xmax = vmax

        plt.xlim(xmin,xmax)
        plt.ylim(xtrmdv[0],xtrmdv[1])

        plt.yticks(tpos, tlabs)
    else:
        for il in range(Ntraj):
            Nvals = len(vardv[il])
            plt.plot(vardv[il], varvv[il], linekinds[il], marker=pointkinds[il],     \
              label=linesn[il], color=colvalues[il], linewidth=lwidths[il],          \
              markersize=psizes[il], markevery=ptf)

# Old way
#            if ptf is None:
#                plt.plot(vardv[il], varvv[il], linekinds[il], label= linesn[il],     \
#                  color=colvalues[il], linewidth=lwidths[il], markersize=psizes[il], markevery=ptf)
#            else:
#                plt.plot(vardv[il], varvv[il], linekinds[il], color=coll[il],        \
#                  linewidth=lwidths[il])
#                plt.plot(vardv[il][::ptf], varvv[il][::ptf], pointkinds[il],         \
#                  label= linesn[il], color=colvalues[il], linewidth=lwidths[il],     \
#                  markersize=psizes[il])

            minvv = np.min(varvv[il])
            maxvv = np.max(varvv[il])
            mindv = np.min(vardv[il])
            maxdv = np.max(vardv[il])

            if minvv < xtrmvv[0]: xtrmvv[0] = minvv
            if maxvv > xtrmvv[1]: xtrmvv[1] = maxvv
            if mindv < xtrmdv[0]: xtrmdv[0] = mindv
            if maxdv > xtrmdv[1]: xtrmdv[1] = maxdv

        plt.xlabel(dtit)
        plt.ylabel(vtit + ' (' + vunit + ')')

        if vmin == 'auto':
            ymin = xtrmvv[0]
        else:
            ymin = vmin
        if vmax == 'auto':
            ymax = xtrmvv[1]
        else:
            ymax = vmax

        plt.xlim(xtrmdv[0],xtrmdv[1])
        plt.ylim(ymin,ymax)

#        plt.xlim(np.min(varTdv),np.max(varTdv))
#        plt.ylim(np.min(varTvv),np.max(varTvv))
        plt.xticks(tpos, tlabs)

    figname = 'lines_time'
    graphtit = gtit.replace('_','\_').replace('&','\&')

    plt.title(graphtit)
    if linesn0 is not None:
        if Ntraj < 10:
            plt.legend(loc=gloc, prop={'size':gsize})
        elif 10 < Ntraj < 20:
            plt.legend(loc=gloc, prop={'size':int(gsize*0.8)})
        else:
            plt.legend(loc=gloc, prop={'size':int(gsize*0.6)})

    print plt.xlim(),':', plt.ylim()
    
    output_kind(kfig, figname, True)

    return

def plot_barbs(xvals,yvals,uvals,vvals,vecfreq,veccolor,veclength,windn,wuts,mapv,graphtit,kfig,figname):
    """ Function to plot wind barbs
      xvals= x position of the values
      yvals= y position of the values
      uvals= values for the x-wind
      uvals= values for the y-wind
      vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 
        'auto', computed automatically to have 20 vectors along each axis)
      veccolor= color of the vectors (None, for 'red')
      veclength= length of the wind barbs (None, for 9)
      windn= name of the wind variable in the graph
      wuts= units of the wind variable in the graph
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      graphtit= title of the graph ('|', for spaces)
      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_barbs'
 
    dx=xvals.shape[1]
    dy=xvals.shape[0]

# Frequency of vectors
    if vecfreq is None:
        xfreq = 1
        yfreq = 1
    elif vecfreq == 'auto':
        xfreq = dx/20
        yfreq = dy/20
    else:
        xfreq=int(vecfreq.split('@')[0])
        yfreq=int(vecfreq.split('@')[1])

    if veccolor == 'auto':
        vcolor = "red"
    else:
        vcolor = veccolor

    if veclength == 'auto':
        vlength = 9
    else:
        vlength = veclength

    plt.rc('text', usetex=True)

    if not mapv is None:
        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
        lat00 = yvals[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = np.min(xvals[::yfreq,::xfreq])
        xlon = np.max(xvals[::yfreq,::xfreq])
        nlat = np.min(yvals[::yfreq,::xfreq])
        xlat = np.max(yvals[::yfreq,::xfreq])

        lon2 = xvals[dy/2,dx/2]
        lat2 = yvals[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    plt.barbs(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq], uvals[::yfreq,::xfreq],\
      vvals[::yfreq,::xfreq], color=vcolor, pivot='tip')

    plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',           \
      xy=(0.85,-0.10), xycoords='axes fraction', color=vcolor)

    plt.title(graphtit.replace('|',' ').replace('&','\&'))

## NOT WORKING ##

# No legend so it is imposed
##    windlabel=windn.replace('_','\_') +' (' + units_lunits(wuts[1]) + ')'
##    vecpatch = mpatches.Patch(color=vcolor, label=windlabel)

##    plt.legend(handles=[vecpatch])

##    vecline = mlines.Line2D([], [], color=vcolor, marker='.', markersize=10, label=windlabel)
##    plt.legend(handles=[vecline], loc=1)

    output_kind(kfig, figname, True)

    return

def plot_ptZvals(vname,vunits,points,ptype,ptsize,graphlims,minmax,figtitle,cbar,    \
  mapv,kfig):
    """ Function to plot a given list of points and values 
      vname= name of the variable in the graph
      vunits= units of the variable
      points= [lon,lat,val] matrix of values
      ptype= type of the point
      ptsize= size of the point
      graphlims= minLON,minLAT,maxLON,maxLAT limits of the graph, None for the full size
      minmax= minimum and maximum type
        'auto': values taken from the extrems of the data
        [min],[max]: given minimum and maximum values
      figtitle= title of the figure
      cbar= color bar
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert-conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      kfig= kind of figure
    """
    fname = 'plot_ptZvals'

    figname = 'pointsZval'

    minlon = points[:,0].min()
    maxlon = points[:,0].max() 

    minlat = points[:,1].min()
    maxlat = points[:,1].max()

    minval = points[:,2].min()
    maxval = points[:,2].max()

#    print 'min/max val;',minval,maxval

    lonrange = (points[:,0] - minlon)/(maxlon - minlon)
    latrange = (points[:,1] - minlat)/(maxlat - minlat)
    colorrange = (points[:,2] - minval)/(maxval - minval)

    plt.rc('text', usetex=True)

    if mapv is not None:
        vlon = points[:,0]
        vlat = points[:,1]
        dx = len(vlon)
        dy = len(vlat)

#        vlon = np.where(vlon[:] < 0., 360. + vlon[:], vlon[:])
#        xvala = np.array(xval)
#        xvala = np.where(xvala < 0., 360. + xvala, xvala)
#        xval = list(xvala)

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        if graphlims is not None:
            nlon = graphlims[0]
            xlon = graphlims[2]
            nlat = graphlims[1]
            xlat = graphlims[3]
        else:
            nlon = np.min(vlon)
            xlon = np.max(vlon)
            nlat = np.min(vlat)
            xlat = np.max(vlat)

        lon2 = vlon[dy/2]
        lat2 = vlat[dy/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": map projecion '" + map_proj + "' not ready!!"
            print '    available: cyl, lcc'
            quit(-1)

#        lons, lats = np.meshgrid(vlon, vlat)
#        lons = np.where(lons < 0., lons + 360., lons)

        x,y = m(vlon,vlat)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True])

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False])
#    else:
#        x = vlon
#        y = vlat
#        plt.xlim(0,dx-1)
#        plt.ylim(0,dy-1)
 
    if minmax == 'auto':
        plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar,    \
          marker=ptype)
    else:
        minv = np.float(minmax.split(',')[0])
        maxv = np.float(minmax.split(',')[1])

        plt.scatter(points[:,0], points[:,1], c=points[:,2], s=ptsize, cmap=cbar,    \
          marker=ptype, vmin=minv, vmax=maxv)

    cbar = plt.colorbar()
    cbar.set_label(vname.replace('_','\_') +' ('+ units_lunits(vunits) + ')')

    plt.title(figtitle)
    if graphlims is not None:
        plt.xlim(graphlims[0], graphlims[2])
        plt.ylim(graphlims[1], graphlims[3])

    output_kind(kfig, figname, True)

    return

#pts = np.zeros((10,3), dtype=np.float)
#pts[:,0] = np.arange(10,20)*1.
#pts[:,1] = np.arange(30,40)*1.
#pts[:,2] = np.arange(-5,5)*1.

#plot_ptZvals('vals','kgm-2',pts,'.',300, 'values of values', 'seismic', 'cyl,l', 'pdf')

def plot_ZQradii(Zmeans, graphtit, kfig, figname):
    """ Function to plot following radial averages only at exact grid poins
      Zmeans= radial means
      radii= values of the taken radii
      graphtit= title of the graph ('|', for spaces)
      kfig= kind of figure
      figname= name of the figure
    """

    fname = 'plot_ZQradii'

    output_kind(kfig, figname, True)

    return

def plot_vector(xvals,yvals,uvals,vvals,vecfreq,vecoln,veccolor,veclength,windn,wuts,\
  mapv,graphtit,kfig,figname):
    """ Function to plot vectors
      xvals= values for the x-axis
      yvals= values for the y-axis
      uvals= values for the x-wind
      vvals= values for the y-wind
      vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 
        'auto', computed automatically to have 20 vectors along each axis)
      veccoln= name for the color of the vectors
        'singlecol'@[colorn]: all the vectors same color ('auto': for 'red')
        'wind'@[colorbar]: color of the vectors according to wind speed sqrt(u^2+v^2) and given [colorbar]
        '3rdvar'@[colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
      veccolor= color of the vectors
      veclength= length of the wind vectors:
        'singlecol': 'auto', for 9
        'wind' and '3rdvar': 'auto' length as wind speed, otherwise fix length
      windn= name of the wind variable in the graph
      wuts= units of the wind variable in the graph
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      graphtit= title of the graph ('|', for spaces)
      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_vector'

    dx=xvals.shape[1]
    dy=xvals.shape[0]

# Frequency of vectors
    if vecfreq is None:
        xfreq = 1
        yfreq = 1
    elif vecfreq == 'auto':
        xfreq = dx/20
        yfreq = dy/20
    else:
        xfreq=int(vecfreq.split('@')[0])
        yfreq=int(vecfreq.split('@')[1])

# Vector length
    if veclength == 'auto':
        vlength = 9
    else:
        vlength = veclength

# Colors
    VecN = vecoln.split('@')[0]
    if VecN == 'singlecol':
        if veccolor == 'auto':
            vcolor = "red"
        else:
            vcolor = veccolor
    elif VecN == 'wind' or VecN == '3rdvar':
        vcolor = vecoln.split('@')[1]

    plt.rc('text', usetex=True)

    if not mapv is None:
        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
        lat00 = yvals[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = np.min(xvals[::yfreq,::xfreq])
        xlon = np.max(xvals[::yfreq,::xfreq])
        nlat = np.min(yvals[::yfreq,::xfreq])
        xlat = np.max(yvals[::yfreq,::xfreq])

        lon2 = xvals[dy/2,dx/2]
        lat2 = yvals[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    if VecN == 'singlecol':
        plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq],                   \
          uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], color=vcolor, pivot='middle')
        plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',       \
          xy=(0.80,-0.15), xycoords='axes fraction', color=vcolor)
    else:
        if veclength != 'auto':
            wind = np.sqrt(uvals**2 + vvals**2)
            uvals = np.where(wind == 0., 0., uvals)
            vvals = np.where(wind == 0., 0., vvals)
            uvals = np.float(veclength)*uvals/wind
            vvals = np.float(veclength)*vvals/wind

        plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq],                   \
          uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], 
          veccolor[::yfreq,::xfreq], cmap=plt.get_cmap(vcolor), pivot='middle')
        cbar = plt.colorbar()

        if VecN == 'wind':
            cbar.set_label('$\sqrt{u^{2} + v^{2}}$ (' + units_lunits(wuts) + ')')
        else:
            vN = vecoln.split('@')[2]
            vU = vecoln.split('@')[3]
            cbar.set_label(vN + ' (' + units_lunits(vU) + ')')

        plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',       \
          xy=(0.80,-0.15), xycoords='axes fraction', color='black')

    plt.title(graphtit.replace('|',' ').replace('&','\&'))

    output_kind(kfig, figname, True)

    return

def plot_basins(xvals,yvals,fvals,vecfreq,vecoln,veccolor,veclength,vcolmin,vcolmax,windn,wuts,\
  mapv,drawcountry,basinid,graphtit,kfig,figname):
    """ Function to plot vectors
      xvals= values for the x-axis
      yvals= values for the y-axis
      fvals= values for the flow (1-8: N, NE, E, ..., NW; 97: sub-basin; 98: small to sea; 99: large to sea)
      vecfreq= [xfreq],[yfreq] frequency of values allong each axis (None, all grid points; 
        'auto', computed automatically to have 20 vectors along each axis)
      veccoln= name for the color of the vectors (as '3rdvar@[basinsvar]@[varn]@[units]' from plot_vector)
      veccolor= color of the vectors
      veclength= length of the wind vectors:
        'singlecol': 'auto', for 9
        'wind' and '3rdvar': 'auto' length as wind speed, otherwise fix length
      windn= name of the wind variable in the graph
      wuts= units of the wind variable in the graph
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      drawcountry= whether country lines should be plotted or not
      basinid= whether  basins id should be plotted or not
      graphtit= title of the graph ('|', for spaces)
      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_basins'

    dx=xvals.shape[1]
    dy=xvals.shape[0]

# Frequency of vectors
    if vecfreq is None:
        xfreq = 1
        yfreq = 1
    elif vecfreq == 'auto':
        xfreq = dx/20
        yfreq = dy/20
    else:
        xfreq=int(vecfreq.split('@')[0])
        yfreq=int(vecfreq.split('@')[1])

# Vector length
    if veclength == 'auto':
        vlength = 9
    else:
        vlength = veclength

# flow direction
    angle = (fvals[::yfreq,::xfreq] - 1)*np.pi/4
    uvals = np.where(fvals[::yfreq,::xfreq] < 9, np.float(veclength)*np.sin(angle), 0.)
    vvals = np.where(fvals[::yfreq,::xfreq] < 9, np.float(veclength)*np.cos(angle), 0.)

# Colors
    vcolor = vecoln.split('@')[0]

    plt.rc('text', usetex=True)

    ddx = np.abs(xvals[dy/2+1,dx/2] - xvals[dy/2,dx/2])
    ddy = np.abs(yvals[dy/2,dx/2+1] - yvals[dy/2,dx/2])
    fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3}

# Setting up colors for each label
#   From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar
    my_cmap = plt.cm.get_cmap(vcolor)
#    vcolmin = np.min(veccolor[::yfreq,::xfreq])
#    vcolmax = np.max(veccolor[::yfreq,::xfreq])
#    vcolmin = 0.
#    vcolmax = 2500.
    norm = mpl.colors.Normalize(vcolmin, vcolmax)
    print 'min col:',vcolmin,'max col:',vcolmax

    xlabpos = []
    ylabpos = []
    labels = []
    labcol = []
    flow = []
    flowvals = []
    
    for j in range(0,dy,yfreq):
        for i in range(0,dx,xfreq):
            if veccolor[j,i] != '--':
                xlabpos.append(xvals[j,i])
                ylabpos.append(yvals[j,i])
                labels.append(int(veccolor[j,i]))
                labcol.append(my_cmap(norm(veccolor[j,i])))
                flowvals.append(fvals[j,i])

    if not mapv is None:
        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
        lat00 = yvals[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = np.min(xvals[::yfreq,::xfreq])
        xlon = np.max(xvals[::yfreq,::xfreq])
        nlat = np.min(yvals[::yfreq,::xfreq])
        xlat = np.max(yvals[::yfreq,::xfreq])

        lon2 = xvals[dy/2,dx/2]
        lat2 = yvals[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()
        drawcountry = True
        if drawcountry: m.drawcountries()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    if veclength != 'auto':
        wind = np.sqrt(uvals**2 + vvals**2)
        uvals = np.where(wind == 0., 0., uvals)
        vvals = np.where(wind == 0., 0., vvals)
        uvals = np.float(veclength)*uvals/wind
        vvals = np.float(veclength)*vvals/wind

    vecolorvals = veccolor[::yfreq,::xfreq]
    vecolorvals = np.where(vecolorvals < vcolmin, vcolmin - 1, vecolorvals)
    vecolorvals = np.where(vecolorvals > vcolmax, vcolmax + 1, vecolorvals)

#    plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq],                       \
#      uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], veccolor[::yfreq,::xfreq],     \
#       cmap=plt.get_cmap(vcolor), pivot='middle')
    plt.quiver(xvals[::yfreq,::xfreq], yvals[::yfreq,::xfreq],                       \
      uvals[::yfreq,::xfreq], vvals[::yfreq,::xfreq], vecolorvals,                   \
       cmap=plt.get_cmap(vcolor), pivot='middle')
    cbar = plt.colorbar()

    vN = vecoln.split('@')[1]
    vU = vecoln.split('@')[2]
    cbar.set_label(vN + ' (' + units_lunits(vU) + ')')

    if basinid:
        for i in range(len(xlabpos)):
            plt.text(xlabpos[i]+0.5*ddx, ylabpos[i]+0.95*ddy, labels[i],             \
              color=labcol[i], fontdict=fontcharac)

# Sea-flow
    for i in range(len(xlabpos)):
        if flowvals[i] == 97:
            plt.plot(xlabpos[i], ylabpos[i], 'x', color=labcol[i])
        elif flowvals[i] == 98:
            plt.plot(xlabpos[i], ylabpos[i], '*', color=labcol[i])
        elif flowvals[i] == 99:
            plt.plot(xlabpos[i], ylabpos[i], 'h', color=labcol[i])

    plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',           \
      xy=(0.80,-0.15), xycoords='axes fraction', color='black')

    plt.title(graphtit.replace('|',' ').replace('&','\&'))

    output_kind(kfig, figname, True)

    return

def plot_river_pattern(xvals, yvals, fvals, veccolor, descid, desclon, desclat,      \
  stn, stlon, stlat, drawcountry, drawbasinid, colordescid, coloriver, colornoriver, \
  colorst, stsize, graphtit, mapv, kfig, figname):
    """ Function to plot vectors
      xvals= values for the x-axis
      yvals= values for the y-axis
      fvals= values for the flow (1-8: N, NE, E, ..., NW; 97: sub-basin; 98: small to sea; 99: large to sea)
      veccolor= basin id
      descid= description id
, desclon, desclat,
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      drawcountry= whether country lines should be plotted or not
      drawbasinid= whether basins id should be plotted or not
      graphtit= title of the graph ('|', for spaces)
      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_river_pattern'

    dx=xvals.shape[1]
    dy=xvals.shape[0]
    xmin = np.min(xvals)
    xmax = np.max(xvals)
    ymin = np.min(yvals)
    ymax = np.max(yvals)

# flow direction

    plt.rc('text', usetex=True)

    ddx = np.abs(xvals[dy/2+1,dx/2] - xvals[dy/2,dx/2])
    ddy = np.abs(yvals[dy/2,dx/2+1] - yvals[dy/2,dx/2])
    fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 3}

    xlabpos = []
    ylabpos = []
    labels = []
    labcol = []
    flow = []
    flowvals = []
    xtrack = []
    ytrack = []
    colortrack = []
    trackid = []
    stnlL = []
    stlonlL = []
    stlatlL = []

# Setting up colors for each label
#   From: http://stackoverflow.com/questions/15235630/matplotlib-pick-up-one-color-associated-to-one-value-in-a-colorbar
    my_cmap = plt.cm.get_cmap(colordescid)
    vcolmin = np.min(descid)
    vcolmax = np.max(descid)

    norm = mpl.colors.Normalize(vcolmin, vcolmax)

# Vector angles
    lengthtrac = np.float(xvals[0,1] - xvals[0,0])
    lengthtrac2 = lengthtrac*np.sqrt(2.)

    for j in range(0,dy):
        for i in range(0,dx):
            if veccolor[j,i] != '--':
#                ibeg = xvals[j,i]-lengthtrac/2.
#                jbeg = yvals[j,i]-lengthtrac/2.
                ibeg = xvals[j,i]
                jbeg = yvals[j,i]
                labels.append(int(veccolor[j,i]))
                flowvals.append(fvals[j,i])
                angle = (fvals[j,i] - 1)*np.pi/4
                if gen.searchInlist([2,4,6,8], fvals[j,i]):
                    iend = ibeg + lengthtrac2*np.sin(angle)
                    jend = jbeg + lengthtrac2*np.cos(angle)
                elif gen.searchInlist([1,3,5,7], fvals[j,i]):
                    iend = ibeg + lengthtrac*np.sin(angle)
                    jend = jbeg + lengthtrac*np.cos(angle)
                else:
                    ibeg = xvals[j,i]
                    jbeg = yvals[j,i]
                    iend = None
                    jend = None

                xlabpos.append(ibeg)
                ylabpos.append(jbeg)
                xtrack.append(ibeg)
                xtrack.append(iend)
                xtrack.append(None)
                ytrack.append(jbeg)
                ytrack.append(jend)
                ytrack.append(None)
                if len(desclon.shape) == 2:
                    difflonlat = np.sqrt((desclon-xvals[j,i])**2 + (desclat-yvals[j,i])**2)
                    mindiffLl = np.min(difflonlat)
                    ilatlon = gen.index_mat(difflonlat, mindiffLl)
                else:
                    ilatlon = np.zeros((2), dtype=int)
                    difflon = np.abs(desclon - xvals[j,i])
                    mindiffl = np.min(difflon)
                    ilatlon[1] = gen.index_vec(difflon,mindiffl)
                    difflat = np.abs(desclat - yvals[j,i])
                    mindiffL = np.min(difflat)
                    ilatlon[0] = gen.index_vec(difflat,mindiffL)

                if descid.mask[ilatlon[0],ilatlon[1]]: 
                    labcol.append(colornoriver)
                    colortrack.append(colornoriver)
                else:
                    if veccolor[j,i] != 6:
                        labcol.append(colornoriver)
                        colortrack.append(colornoriver)
#                        print 'Lluis: veccol:', veccolor[j,i], 'mindiffl', mindiffl,'mindiffL:',mindiffL
                    else:
#                        print 'Lluis right! mindiffl', mindiffl,'mindiffL:',mindiffL
                        labcol.append(coloriver)
                        colortrack.append(my_cmap(norm(descid[ilatlon[0],ilatlon[1]])))

                trackid.append(descid[ilatlon[0],ilatlon[1]])

    totvals = len(flowvals)

    if not mapv is None:
        lon00 = np.where(xvals[:] < 0., 360. + xvals[:], xvals[:])
        lat00 = yvals[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        nlon = np.min(xvals)
        xlon = np.max(xvals)
        nlat = np.min(yvals)
        xlat = np.max(yvals)

        lon2 = xvals[dy/2,dx/2]
        lat2 = yvals[dy/2,dx/2]

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()
        drawcountry = True
        if drawcountry: m.drawcountries()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    j = 0
    for i in range(len(xlabpos)):
        plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i])
        j = j + 3

# Sea-flow
    for i in range(len(xlabpos)):
        if flowvals[i] == 97:
            plt.plot(xlabpos[i], ylabpos[i], 'x', color=colortrack[i])
        elif flowvals[i] == 98:
            plt.plot(xlabpos[i], ylabpos[i], '*', color=colortrack[i])
        elif flowvals[i] == 99:
            plt.plot(xlabpos[i], ylabpos[i], 'h', color=colortrack[i])

    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)

    if drawbasinid:
        for i in range(len(xlabpos)):
            plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, labels[i],    \
              color=labcol[i], fontdict=fontcharac)
#            plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, trackid[i],   \
#              color=colortrack[i], fontdict=fontcharac)
#            plt.text(xlabpos[i]+0.1*lengthtrac, ylabpos[i]+0.1*lengthtrac, str(flowvals[i]),   \
#              color=colortrack[i], fontdict=fontcharac)

    for ist in range(len(stlon)):
        if stlon[ist] > xmin and stlon[ist] < xmax and stlat[ist] > ymin and stlat[ist] < ymax:
            stname = '*'
            for ic in range(len(stn[ist,:])):
                stname  = stname + stn[ist,ic]
#            plt.annotate(stname, xy=(stlon[ist],stlat[ist]), color=colorst, )
            plt.text(stlon[ist], stlat[ist], stname, color=colorst, fontsize=stsize)

#    cbar = plt.colorbar()
    plt.title(graphtit.replace('&','\&'))

    output_kind(kfig, figname, True)

    return

def var_3desc(ovar):
    """ Function to provide std_name, long_name and units from an object variable
      ovar= object variable
    """
    fname = 'var_desc'
    varattrs = ovar.ncattrs()

    if searchInlist(varattrs,'std_name'):
        stdn = ovar.getncattr('std_name')
    else:
        vvalues = variables_values(ovar._name)
        stdn = vvalues[1]

    if searchInlist(varattrs,'long_name'):
        lonn = ovar.getncattr('long_name')
    else:
        lonn = vvalues[4].replace('|',' ')

    if searchInlist(varattrs,'units'):
        un = ovar.getncattr('units')
    else:
        un = vvalues[5]

    return stdn, lonn, un

def plot_river_desc(lons, lats, rns, rss, rus, ros, bcolor, ucolor, uuts, mapv,      \
  graphtit, kfig, lloc, figname):
    """ Function to plot rivers from 'river_desc.nc' ORCDHIEE
      lons= values for the x-axis
      lats= values for the y-axis
      rns= list of the name of the rivers to plot
      rss= dictionary with the lon,lats of the subbasins of each river
      rus= dictionary with the lon,lats of the upstream values of each river
      ros= dictionary with the lon,lats of the outflow points of each river
      bcolor= color of the lines for the subbasins
      ucolor= bar color for the upstreams
      uuts= units of the upstream
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      graphtit= title of the graph ('|', for spaces)
      lloc= location of the legend (0, automatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_basins'

    colvalues, linekinds, pointkinds, lwidths, psizes = ColorsLinesPointsStyles(     \
      len(rns), bcolor, ['-'], ['.'], [1.], [1.],  None)

    dx=lons.shape[1]
    dy=lats.shape[0]

    plt.rc('text', usetex=True)

    nlon = np.min(lons)
    xlon = np.max(lons)
    nlat = np.min(lats)
    xlat = np.max(lats)

    dlon = xlon - nlon
    dlat = xlat - nlat

# Making bigger the area to map
    nlon = np.min(lons)-dlon*0.1
    xlon = np.max(lons)+dlon*0.1
    nlat = np.min(lats)-dlat*0.1
    xlat = np.max(lats)+dlat*0.1 

    plt.xlim(nlon,xlon)
    plt.ylim(nlat,xlat)

    if not mapv is None:
        lon00 = np.where(lons[:] < 0., 360. + lons[:], lons[:])
        lat00 = lats[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        lon2 = (nlon + xlon)/2.
        lat2 = (nlat + xlat)/2.

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    umin = 10.e20
    umax = -umin
    for rn in rns:
        nrns = np.min(rus[rn])
        xrns = np.max(rus[rn])

        if umin > nrns: umin = nrns
        if umax < xrns: umax = xrns

    i = 0
    for rn in rns:
        plt.pcolormesh(lons, lats, rus[rn], vmin=umin, vmax=umax,                    \
          cmap=plt.get_cmap(ucolor), norm=matplotlib.colors.LogNorm())
#        plt.pcolormesh(lons, lats, rss[rn], cmap=plt.get_cmap(ucolor))
        if rn == rns[0]: cbar = plt.colorbar()
        riverarea = rss[rn]

        riverarea.mask = ma.nomask       
        rarea = np.where(riverarea == 999999999, -5, riverarea)
        rarea = np.where(rarea > 0, 1, 0)
        plt.contour(lons, lats, rarea, levels = [0.], colors=colvalues[i],           \
          linewidths=1.5)

        oflow = ros[rn]
        plt.plot(oflow[0], oflow[1], 'h', color=colvalues[i], label=rn)
#        plt.annotate(rn, xy=(0., i*0.05), xycoords='axes fraction', color=colvalues[i])
        i = i+1


# Attempts to plot sub-basins

# Gradient subbasins (not working)
#        gradriver = np.zeros((dy,dx), dtype=int)
#        for j in range(dy-2):
#            for i in range(dx-2):
#                gradriver[j,i] = riverarea[j,i+1]-riverarea[j,i]+riverarea[j+1,i]-  \
#                  riverarea[j,i]
#        garea = np.where(gradriver != 0, 1, gradriver)
#        print garea
#        plt.contour(lons, lats, garea, levels = [0.], colors='red', linewidths=1.)

# Not working
#        subbasins = np.sort(np.unique(rarea))
#        for subb in subbasins:
#            if subb > 0.:
#                sarea = np.where(rarea != subb, 0, rarea)
#                plt.contour(lons, lats, sarea, levels = [1.], colors='gray', linewidths=1.)

    cbar.set_label('upstream (' + units_lunits(uuts) + ')')
    plt.legend(loc=lloc)

#    plt.annotate(windn.replace('_','\_') +' (' + units_lunits(wuts) + ')',           \
#      xy=(0.80,-0.15), xycoords='axes fraction', color='black')

    plt.title(graphtit.replace('|',' ').replace('&','\&'))

    output_kind(kfig, figname, True)

    return

def plot_vertical_lev(vertz, vertp, zlog, dzlog, plog, dplog, gtit, kfig, lloc):
    """ plotting vertical levels distribution
    plot_vertical_lev(vertz, gtit, kfig, lloc)
      vertz= distribution of vertical heights [z]
      vertp= distribution of vertical pressures [Pa]
      zlog: to use logarithmic scale on the height axis ('true/false')
      dzlog: to use logarithmic scale on the difference of height between levels axis ('true/false')
      plog: to use logarithmic scale on the pressure axis ('true/false')
      pzlog: to use logarithmic scale on the difference of pressure between levels axis ('true/false')
      gtit= title of the graph ('!' for spaces)
      kfig= kind of figure (jpg, pdf, png)
      lloc= location of the legend (0, automatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'      kfig= kind of figure
    """
    fname = 'plot_vertical_lev'

    figname = 'vertical_lev'

    if vertz == 'h':
        print fname + '_____________________________________________________________'
        print plot_vertical_lev.__doc__
        quit()

    Nlev = len(vertz)
    dvertz = vertz[1:Nlev] - vertz[0:Nlev-1]

# From: http://stackoverflow.com/questions/14762181/adding-a-y-axis-label-to-secondary-y-axis-in-matplotlib

    plt.rc('text', usetex=True)
# Height plot
##
    if vertp is not None: 
        print plt.subplots.__doc__
        fig, (ax1, ax3) = plt.subplots(nrows=2, ncols=1, sharex=True)
    else:
        fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()

    plt.xlim(0,Nlev)
 
    l1 = ax1.plot(range(Nlev), vertz, 'r-x', label='height')
    l2 = ax2.plot(range(1,Nlev), dvertz, 'b-x', label='dheight')

#    ax1.set_xlabel('level (\#)')
    ax1.set_ylabel('height (m)', color='r')
    ax1.set_ylim(1,np.max(vertz))
    ax1.grid()
    ax2.set_ylabel('difference between levels (m)', color='b')

    if vertp is not None: ax1.set_title('height')
    if zlog: ax1.set_yscale('log')
    if dzlog: ax2.set_yscale('log')

# Pressure plot
##
    if vertp is not None:

        Nlev = len(vertp)
        dvertp = vertp[0:Nlev-1] - vertp[1:Nlev]

        ax4 = ax3.twinx()

        l3 = ax3.plot(range(Nlev), vertp, 'r-o', label='pressure')
        l4 = ax4.plot(range(1,Nlev), dvertp, 'b-o', label='dpressure')

        ax3.set_xlim(0, Nlev)
        ax3.set_xlabel('level (\#)')
        ax3.set_ylabel('pressure (Pa)', color='r')
        ax3.set_ylim(np.min(vertp), np.max(vertp))
        ax3.grid()
        ax4.set_ylabel('difference between levels (Pa)', color='b')

        if plog: ax3.set_yscale('log')
        if dplog: ax4.set_yscale('log')

    if vertp is not None:
        fig.suptitle(gtit)
        ax3.set_title('pressure')
    else:
        plt.title(gtit)

    output_kind(kfig, figname, True)

    return

def plot_subbasin(subname, lons, lats, rsf, rLlf, rof, rcolors, mapv, basinwidth,    \
  drawsubid, graphtit, kfig, lloc, figname):
    """ Function to plot rivers from 'river_desc.nc' ORCDHIEE
      subname= name of the subbasin
      lons= values for the x-axis
      lats= values for the y-axis
      rsf= list of the name of the sub-flows to plot
      rLlf= dictionary with the lat,lon of the subflows
      rof= dictionary with the outflows of the subflows
      rcolors= base of colors of the lines for the subflows (first level)
      mapv= map characteristics: [proj],[res]
        see full documentation: http://matplotlib.org/basemap/
        [proj]: projection
          * 'cyl', cilindric
          * 'lcc', lambert conformal
        [res]: resolution:
          * 'c', crude
          * 'l', low
          * 'i', intermediate
          * 'h', high
          * 'f', full
      basinwidth = with of the lines of the basin
      drawsubid= wehther sub-flow ids should be plot or not
      graphtit= title of the graph ('|', for spaces)
      lloc= location of the legend (0, automatic)
        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
        9: 'upper center', 10: 'center'      kfig= kind of figure
      figname= name of the figure
    """
    fname = 'plot_subbasin'

    dx=lons.shape[1]
    dy=lats.shape[0]

    lengthtrac = np.float(lons[0,1] - lons[0,0])
    lengthtrac2 = lengthtrac*np.sqrt(2.)

    ddx = np.abs(lons[dy/2+1,dx/2] - lons[dy/2,dx/2])
    ddy = np.abs(lats[dy/2,dx/2+1] - lats[dy/2,dx/2])
    fontcharac = {'family': 'serif', 'weight': 'normal', 'size': 10}

    subhr = gen.chainSnumHierarchy(rsf)

# Getting colors
    Nsubs = len(rsf)
    maxNgrids = -1
    
    collev = {}
    if rcolors is None:
# 1st level colors from colorsauto
        for isub in range(subhr.xlevs[1]):
            collev[gen.num_chainSnum(isub+1)] = gen.colhex_dec(colorsauto[isub])
    else:
        if len(rcolors) != subhr.Llev1:
            print errormsg
            print '  ' + fname + ': number of first subbasins:', subhr.Llev1,        \
              'and number of provided colors:', len(rcolors), 'differ !!'
            print '    provided colors:', rcolors
            quit(-1)
        for isub in range(subhr.xlevs[1]):
            collev[gen.num_chainSnum(isub+1)] = rcolors[isub]

# All the colors should finish white...
    endcol = [1., 1., 1.]
#    endcol = [0.9, 0.9, 0.9]
    colorsub = {}
    xlabpos = []
    ylabpos = []
    xtrack = []
    ytrack = []
    colortrack = []
    outflows = []

    for isub in rsf:
        levs = gen.num_split(isub)
        Nlevs = len(levs)
        if Nlevs == 1:
            begcol = [0., 0., 0.]
        else:
            begcol = collev[levs[1]]
        for ilev in range(2,Nlevs):
# We do not want a last white value...
            if Nlevs < subhr.Nlevs: 
                maxNcol = subhr.xlevs[ilev]+1
            else:
                maxNcol = subhr.xlevs[ilev]+1
            colors = gen.color_deg([begcol, endcol], maxNcol, 'std')
            begcol = colors[1:4,gen.chainSnum_num(levs[ilev])-1]

        colorsub[isub] = tuple(begcol)
# Avoiding too clear colors
        cols = colorsub[isub]
        if (cols[0] >= 0.95)*(cols[1] >= 0.95)*(cols[2] >= 0.95):
            print warnmsg
            print '  ' + fname + ': avoiding to clear color:', colorsub[isub], '!!'
            maxcol = np.max(cols)
            imaxcol = gen.index_vec(cols, maxcol)
            clear = np.ones((3), dtype=np.float)*0.9
            clear[imaxcol] = 0.95
            print '    keeping it clear:', clear
            colorsub[isub] = tuple(clear)

        Llfs = rLlf[isub]
        outfs = rof[isub]
        for igrid in range(len(outfs)):
            Llf = Llfs[igrid]
            ibeg = Llf[1]
            jbeg = Llf[0]
            outflows.append(outfs[igrid])
            angle = (outfs[igrid] - 1)*np.pi/4
            if gen.searchInlist([2,4,6,8], outfs[igrid]):
                iend = ibeg + lengthtrac2*np.sin(angle)
                jend = jbeg + lengthtrac2*np.cos(angle)
            elif gen.searchInlist([1,3,5,7], outfs[igrid]):
                iend = ibeg + lengthtrac*np.sin(angle)
                jend = jbeg + lengthtrac*np.cos(angle)
            else:
                ibeg = xvals[j,i]
                jbeg = yvals[j,i]
                iend = None
                jend = None

            xtrack.append(ibeg)
            xtrack.append(iend)
            xtrack.append(None)
            ytrack.append(jbeg)
            ytrack.append(jend)
            ytrack.append(None)
            colortrack.append(colorsub[isub])

        xlabpos.append(ibeg)
        ylabpos.append(jbeg)

    plt.rc('text', usetex=True)

    nlon = np.min(xlabpos)
    xlon = np.max(xlabpos)
    nlat = np.min(ylabpos)
    xlat = np.max(ylabpos)

    dlon = xlon - nlon
    dlat = xlat - nlat

# Making bigger the area to map
#    nlon = nlon-dlon-lengthtrac2
#    xlon = xlon+dlon+lengthtrac2
#    nlat = nlat-dlat-lengthtrac2
#    xlat = xlat+dlat+lengthtrac2
    nlon = nlon-2.*lengthtrac*dlon/np.abs(dlon)
    xlon = xlon+2.*lengthtrac*dlon/np.abs(dlon)
    nlat = nlat-2.*lengthtrac*dlat/np.abs(dlat)
    xlat = xlat+2.*lengthtrac*dlat/np.abs(dlat)

    plt.xlim(nlon,xlon)
    plt.ylim(nlat,xlat)

    if not mapv is None:
        lon00 = np.where(lons[:] < 0., 360. + lons[:], lons[:])
        lat00 = lats[:]

        map_proj=mapv.split(',')[0]
        map_res=mapv.split(',')[1]

        lon2 = (nlon + xlon)/2.
        lat2 = (nlat + xlat)/2.

        print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:',     \
          xlon, ',', xlat

        if map_proj == 'cyl':
            m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat,         \
              urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        elif map_proj == 'lcc':
            m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
              llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
        else:
            print errormsg
            print '  ' + fname + ": projection '" + map_proj + "' not ready!!"
            print '    projections available: cyl, lcc'
            quit(-1)

        m.drawcoastlines()

        meridians = pretty_int(nlon,xlon,5)
        m.drawmeridians(meridians,labels=[True,False,False,True],color="black")

        parallels = pretty_int(nlat,xlat,5)
        m.drawparallels(parallels,labels=[False,True,True,False],color="black")

        plt.xlabel('W-E')
        plt.ylabel('S-N')

    j = 0
    for i in range(len(xtrack)/3):
        plt.plot(xtrack[j:j+2], ytrack[j:j+2], color=colortrack[i], linewidth=basinwidth)
        j = j + 3

# Sea-flow
    for i in range(len(outflows)):
        if outflows[i] == 97:
            plt.plot(xlabpos[i], ylabpos[i], 'x', color=colorsub[i])
        elif outflows[i] == 98:
            plt.plot(xlabpos[i], ylabpos[i], '*', color=colorsub[i])
        elif outflows[i] == 99:
            plt.plot(xlabpos[i], ylabpos[i], 'h', color=colorsub[i])

    if drawsubid:
        for i in range(len(xlabpos)):
            plt.text(xlabpos[i]+0.05*lengthtrac, ylabpos[i]+0.05*lengthtrac, rsf[i], \
              color=colorsub[rsf[i]], fontdict=fontcharac)

    plt.title(graphtit)
    output_kind(kfig, figname, True)

    return