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

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

fillValue = 1.e20

####### Funtions
# searchInlist:
# datetimeStr_datetime:
# dateStr_date:
# numVector_String:
# timeref_datetime:
# slice_variable:
# interpolate_locs:
# datetimeStr_conversion:
# percendone:
# netCDFdatetime_realdatetime:
# file_nlines:
# variables_values:
# check_colorBar:
# units_lunits:
# ASCII_LaTeX:
# pretty_int:
# DegGradSec_deg:
# intT2dt:
# 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

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

    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]
    """
    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 = []

    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(int(varobj.shape[idd])-1)
                    else:
                        varvalsdim.append(int(dimcutv))
                found = True
                break
        if not found and not searchInlist(dimnslice,vardims[idd]):
            varvalsdim.append(slice(0,varobj.shape[idd]))
            dimnslice.append(vardims[idd])

    varvalues = varobj[tuple(varvalsdim)]

    return varvalues, dimnslice

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 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', '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 == '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 == '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 == '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 == '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/(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 == 'm2/s': lu='$m2s^{-1}$'
        elif u == 'm2s-1': lu='$m2s^{-1}$'
        elif u == 'm2/s2': lu='$m2s^{-2}$'
        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='\%'
        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.

    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]
            yeari = firstTvec[0]

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

            datev = firstTdt
            while datev < lastTdt:
                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]

            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)
            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]

            if np.sum(firstTvec[4:5]) > 0 or np.mod(hour,Nunits) != 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:
                    firstTdt = dt.datetime(yr, mon, day, Nunits, 0, 0)
                    datev = dt.datetime(yr, mon, day, 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]

            if firstTvec[5] > 0 or np.mod(minu,Nunits) != 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:
                    firstTdt = dt.datetime(yr, mon, day, hour, Nunits, 0)
                    datev = dt.datetime(yr, mon, day, hour, 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
            it=0
            while datev <= lastTdt:
                datev = firstTdt + dt.timedelta(seconds=Nunits)
                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,Nargs,args,char,expectargs):
    """ 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
      Nargs= theoretical number of arguments
      args= passed arguments
      char= character used to split the arguments
    """

    fname = 'check_arguments'

    Nvals = len(args.split(char))
    if Nvals != Nargs:
        print errormsg
        print '  ' + fname + ': wrong number of arguments:',Nvals," passed to  '",   \
          funcname, "' which requires:",Nargs,'!!'
        print '    given arguments:',args.split(char)
        print '    expected arguments:',expectargs
        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 == '1' or val == 'yes': 
        boolv = True
    elif 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'

    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
    tuA = tunitA.split(' ')[0]
    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 == 'dayss':
            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):
    """ 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 (-1, 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
    """
    fname = 'plot_TimeSeries'

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


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

    Nlines = len(valtimes)

    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,:])
        plt.plot(array[1,:],array[0,:], linekinds[N7lines], label= linesn[il])
        if il == 6: N7lines = N7lines + 1

    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[:])

        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)

#        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:
#        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_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 points
      [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, 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
      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=0)
    
    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:',varv.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,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
      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 not 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
      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 axist == 'x':
        dimxv = timv.copy()
        dimyv = valv.copy()
    else:
        dimxv = valv.copy()
        dimyv = timv.copy()

    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':
        if axist == 'x': 
            axist = 'y'
        else:
            axist = 'x'

    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
        if reva0 == 'transpose':
            print '  reversing the axes of the figure (x-->y, y-->x)!!'
            varsv = np.transpose(varsv)
            varcv = np.transpose(varcv)
            lon0 = np.transpose(lon0)
            lat0 = np.transpose(lat0)

        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 reva0  == 'transpose':
            print '  reversing the axes of the figure (x-->y, y-->x)!!'
            varsv = np.transpose(varsv)
            varcv = np.transpose(varcv)
            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:
                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

    dx=varsv.shape[1]
    dy=varsv.shape[0]
    
    plt.rc('text', usetex=True)

    if axist == 'x':
        valpos = pretty_int(y.min(),y.max(),10)
        vallabels = list(valpos)
        for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])
    else:
        valpos = pretty_int(x.min(),x.max(),10)
        vallabels = list(valpos)
        for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])

    if reva0 == 'flip':
        if reva.split('@')[1] == 'x':
            varsv[:,0:dx-1] = varsv[:,dx-1:0:-1]
            varcv[:,0:dx-1] = varcv[:,dx-1:0:-1]
            plt.xticks(valpos, vallabels[::-1])
        else:
            varsv[0:dy-1,:] = varsv[dy-1:0:-1,:]
            varcv[0:dy-1,:] = varcv[dy-1:0:-1,:]
            plt.yticks(valpos, vallabels[::-1])
    else:
        plt.xlim(0,dx-1)
        plt.ylim(0,dy-1)

    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:
        if axist == 'x':
            plt.xlabel(timeu)
            plt.xticks(timpos, timlab)
            plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(valu) + ')')
            plt.yticks(valpos, vallabels)
        else:
            plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(valu) + ')')
            plt.xticks(valpos, vallabels)
            plt.ylabel(timeu)
            plt.yticks(timpos, timlab)

# set the limits of the plot to the limits of the data
    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

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 = int(dd[idd].split('|')[1])
            if dn == dname:
                if dvalue == -1:
                    slicex.append(slice(0,dxv.shape[ipos]))
                else:
                    slicex.append(dvalue)
                break
        ipos = ipos + 1

    slicey = []
    ipos=0
    for dn in dny:
        for idd in range(len(dd)):
            dname = dd[idd].split('|')[0]
            dvalue = int(dd[idd].split('|')[1])
            if dn == dname:
                if dvalue == -1:
                    slicey.append(slice(0,dyv.shape[ipos]))
                else:
                    slicey.append(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 (-1, 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-']

    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 plot_lines_time(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, tpos, tlabs,     \
  gtit, gloc, kfig, coll):
    """ 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
      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
      gloc= location of the legend (-1, 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
      coll= ',' list of colors for the lines or None for automatic
      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-']

    Ntraj = len(vardv)

    N7lines = 0

    plt.rc('text', usetex=True)
    varTvv = []
    varTdv = []

    if vaxis == 'x':
        for il in range(Ntraj):
            if coll is None:
                plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
            else:
                plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il],\
                  color=coll[il])

            varTvv = varTvv + list(varvv[il])
            varTdv = varTdv + list(vardv[il])
            if il == 6: N7lines = N7lines + 1

        plt.xlabel(vtit + ' (' + vunit + ')')
        plt.ylabel(dtit)
        plt.xlim(np.min(varTvv),np.max(varTvv))
        plt.ylim(np.min(varTdv),np.max(varTdv))
        plt.yticks(tpos, tlabs)
    else:
        for il in range(Ntraj):
            if coll is None:
                plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
            else:
                plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il],\
                  color=coll[il])

            varTvv = varTvv + list(varvv[il])
            varTdv = varTdv + list(vardv[il])
            if il == 6: N7lines = N7lines + 1

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

        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)
    plt.legend(loc=gloc)
    
    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= values for the 'x-axis'
      yvals= values for the 'y-axis'
      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