import numpy as np
from netCDF4 import Dataset as NetCDFFile
import os
import re
import numpy.ma as ma

main = 'nc_var_tools.py'

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

fillValue = 1.e20
fillValueF = 1.e20
fillValueC = '-'
fillValueI = -99999
fillValueI16 = -99999
fillValueF64 = 1.e20
fillValueI32 = -99999

def values_line(line, splitv, chars):
    """ Function to retrieve values from a line of ASCII text removing a series of characters
      line= line of ASCII text to work with
      splitv= character to use to split with the line
      chars= list of characters/strings to remove
    >>> values_line('0.0 0.0 1 [ 0 0 ]', ' ', ['[', ']'])
    ['0.0', '0.0', '1', '0', '0']
    """
    fname = 'values_line'

    vals = line.replace('\n', '')
    for ic in chars:
        vals = vals.replace(ic, '')

    vals0 = vals.split(splitv)
    vals = []
    for val in vals0:
        if not len(val) == 0: vals.append(val)

    return vals

def filexist(filen, emsg, fmsg):
    """ Check existence of a file 
    filen= file name
    emsg= general error message
    fmsg= message before generic text of inexistence of file
    """
    if not os.path.isfile(filen):
        print emsg
        print '  ' + fmsg + ' file "' + filen + '" does not exist !!'
        print emsg
        quit(-1)    

    return

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 reduce_tabulars(string):
    """ Function to give words of a line of text removing any extra space and tabular separation
    """
    values = string.replace('\n','').split('\t')
    vals = []
    for val in values:
         if len(val) > 0:
             vals.append(val)

    return vals

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

    fname = 'check_arguments'

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

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

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

    return

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

    fname = 'Str_Bool'

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

    return boolv

def Nchar(N,string):
    """ Function to provide 'N' times concatenated 'string'
      N= number of copies
      strin= string to reproduce
    >>> Nchar(4,'#')
    ####
    """
    fname = 'Nchar'

    newstring=''
    for it in range(N):
        newstring=newstring + string

    return newstring

def string_dicvals_char(dictionary, string, char):
    """ Function to provide a string taking values from a given [string] where each single character 
      coincide with the key of the dictionary [dictionary], separated by a given set of characters [char]
      dictionary= dictionary of the values dim([c]) = {[value]}
      string = [c1][...[cN]] String with the combination of [c]s
      >>> dictvals = {'i': 'dimx', 'j': 'dimy', 'k': 'dimz', 't': 'dimt'}
      >>> string_dicvals_char(dictvals, 'ijkt', ', ')
      dimx, dimy, dimz, dimt
    """
    fname = 'string_dicvals_char'

    Lstring = len(string)
    newstring = ''
    for ic in range(Lstring):
        cc = string[ic:ic+1]
        if not dictionary.has_key(cc):
            print errormsg
            print '  ' + fname + ": dictionary of values does not have character '" +\
              cc + "' !!"
            print '    it has:',dictionary

        if ic == 0:
            newstring = dictionary[cc]
        else:
            newstring = newstring + char + dictionary[cc]

    return newstring

def substring(stringv,pos,char):
    """ Function to substitute a character of a given string
      stringv= string to change
      pos= position to change
      char= characters to use as substitution
    >>> substring('0123456',2,'ii')
    01ii3456
    """
    fname = 'substring'

# From: http://stackoverflow.com/questions/1228299/change-one-character-in-a-string-in-python
    newstring = stringv[:pos] + char + stringv[pos+1:]

    return newstring

def period_information(idate, edate, totunits):
    """ Function to provide the information of a given period idate, edate (dates in [YYYY][MM][DD][HH][MI][SS] format)
      [idate]= initial date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
      [edate]= end date of the period (in [YYYY][MM][DD][HH][MI][SS] format)
      [totunits]= latest number of entire temporal units to compute: 'year', 'month', 'day', 'hour', 'minute', 'second'

    resultant information given as ',' list of:
      [dT]: sign of the temporal increment ('pos', 'neg')
      [idate]: initial date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
      [edate]: end date as [YYYY]:[MM]:[DD]:[HH]:[MI]:[SS]
      [Nsecs]: total seconds between dates
      [Nyears]: Number of years taken as entire units ([iYYYY]0101000000 to [eYYYY+1]0101000000)
      [ExactYears]: Exact year-dates of the period as a '@' separated list of  [YYYY]0101000000
      [Nmonths]: Number of months taken as entire units ([iYYYY][iMM]01000000 to [eYYYY][eMM+1]01000000)
      [ExactMonths]: Exact mon-dates of the period as  a '@' separated list of [YYYY][MM]01000000
      [Ndays]: Number of days taken as entire units ([iYYYY][iMM][iDD]000000 to [eYYYY][eMM][eDD+1]000000)
      [ExactDays]: Exact day-dates of the period as  a '@' separated list of [YYYY][MM][DD]000000
      [Nhours]: Number of hours taken as entire units ([iYYYY][iMM][iDD][iHH]0000 to [eYYYY][eMM][eDD][eHH+1]0000)
      [ExactHours]: Exact hour-dates of the period as a '@' separated list of [YYYY][MM][DD][HH]0000
      [Nminutes]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI]00 to [eYYYY][eMM][eDD][eHH][eMI+1]00)
      [ExactMinutes]: Exact minutes-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI]00
      [Nsecs]: Number of minutes taken as entire units ([iYYYY][iMM][iDD][iHH][iMI][iS] to [eYYYY][eMM][eDD][eHH][eMI][eSE+1])
      [ExactSeconds]: Exact seconds-dates of the period as a '@' separated list of [YYYY][MM][DD][HH][MI][SS]

    >>> period_information( '19760217083110', '19770828000000', 'month')
    pos,1976:02:17:08:31:10,1977:08:28:00:00:00,48180530.0,19,19760201000000@19760301000000@19760401000000@
    19760501000000@19760601000000@19760701000000@19760801000000@19760901000000@19761001000000@19761101000000@
    19761201000000@19770101000000@19770201000000@19770301000000@19770401000000@19770501000000@19770601000000@
    19770701000000@19770801000000@19770901000000
    """
    import datetime as dt

    fname = 'period_information'

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

    expectargs = '[idate],[edate],[totunits]'
 
    check_arguments(fname,str(idate)+','+str(edate)+','+str(totunits),expectargs,',')

# Checking length
    Lidate = len(idate)
    Ledate = len(edate)

    if Lidate != Ledate:
        print errormsg
        print '   ' + fname + ': string dates of the period have different length !!'
        print "    idate: '" + idate + "' (L=", Lidate, ") edate: '" + edate +       \
          "' (L=", Ledate, ')'
        quit(-1)
    else:
        if Lidate != 14:
            print errormsg
            print '  ' + fname + ": wrong initial date: '" + idate + "' must " +     \
              'have 14 characters!!'
            print '    and it has:', Lidate
            quit(-1)
        if Ledate != 14:
            print errormsg
            print '  ' + fname + ": wrong end date: '" + edate + "' must " +         \
              'have 14 characters!!'
            print '    and it has:', Ledate
            quit(-1)

# Checking for right order of dates
    if int(idate) > int(edate):
        print warnmsg
        print '  ' + fname + ": initial date '" + idate + "' after end date '" +     \
          edate + "' !!"
        print "  re-sorting them!"
        i1 = edate
        e1 = idate
        idate = i1 
        edate = e1
        oinf = 'neg'
    else:
        oinf = 'pos'

# Year, month, day, hour, minute, second initial date
    iyrS = idate[0:4]
    imoS = idate[4:6]
    idaS = idate[6:8]
    ihoS = idate[8:10]
    imiS = idate[10:12]
    iseS = idate[12:14]

    oinf = oinf + ',' + iyrS+':'+imoS+':'+idaS+':'+ihoS+':'+imiS+':'+iseS

# Year, month, day, hour, minute, second end date
    eyrS = edate[0:4]
    emoS = edate[4:6]
    edaS = edate[6:8]
    ehoS = edate[8:10]
    emiS = edate[10:12]
    eseS = edate[12:14]

    oinf = oinf + ',' + eyrS+':'+emoS+':'+edaS+':'+ehoS+':'+emiS+':'+eseS

    idateT = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
    edateT = dt.datetime(int(eyrS),int(emoS),int(edaS),int(ehoS),int(emiS),int(eseS))

# Seconds between dates
    DT = edateT - idateT
    Nsecs = DT.total_seconds()

    oinf = oinf + ',' + str(Nsecs)

# Looking for number of years tacking exact beginning of the units [iYYYY]0101000000, [eYYYY+1]0101000000
##
    if totunits == 'year':
        itime = int(iyrS)
        if edate[4:15] != '0101000000':
            etime = int(eyrS)+1
        else:
            etime = int(eyrS)

        itimeT = dt.datetime(itime,1,1,0,0,0)

        Nyears = 1
        ExactYears = itimeT.strftime("%Y%m%d%H%M%S")
        it = int(iyrS)
        while it + 1 <= etime:
            it = it + 1
            Nyears = Nyears + 1
            itT = dt.datetime(it,1,1,0,0,0)
            ExactYears = ExactYears + '@' + itT.strftime("%Y%m%d%H%M%S")

        oinf = oinf + ',' + str(Nyears) + ',' + ExactYears

# Looking for number of months tacking exact beginning of the units [iYYYY][iMM]01000000, [eYYYY][eMM+1]01000000
##
    if totunits == 'month':
        itime = int(iyrS)*100 + int(imoS)
        if edate[6:15] != '01000000':
            emo1 = int(emoS)+1
            if emo1 > 13:
                eyr1 = str(int(eyrS) + 1)
                emo1 = 01
            else:
                eyr1 = eyrS
            
        else:
            eyr1 = eyrS
            emo1 = emoS
    
        etime = int(eyr1)*100 + int(emo1) 
        itimeT = dt.datetime(int(iyrS),int(imoS),1,0,0,0)
        
        Nmonths = 1
        ExactMonths = itimeT.strftime("%Y%m%d%H%M%S")
        it = itime
        while it + 1 <= etime:
            it = it + 1
            Nmonths = Nmonths + 1
            ityr = it/100
            itmo = it - ityr*100
            if itmo > 12:
                ityr = ityr + 1
                itmo = 1
                it = ityr * 100 + itmo
    
            itT = dt.datetime(ityr,itmo,1,0,0,0)
            ExactMonths = ExactMonths + '@' + itT.strftime("%Y%m%d%H%M%S")
    
        oinf = oinf + ',' + str(Nmonths) + ',' + ExactMonths

# Looking for number of days tacking exact beginning of the units [iYYYY][iMM][iDD]000000, [eYYYY][eMM][eDD+1]000000
##
    if totunits == 'day':
        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),0,0,0)
        if edate[8:15] != '000000':
            etime = dt.datetime(int(eyrS), int(emoS), int(edaS)+1,0,0,0)
        else:
            etime = edateT
    
        Ndays = 1
        ExactDays = itime.strftime("%Y%m%d%H%M%S")
        it = itime
        while it + dt.timedelta(days=1) <= etime:
            it = it + dt.timedelta(days=1)
            Ndays = Ndays + 1
            ExactDays = ExactDays + '@' + it.strftime("%Y%m%d%H%M%S")
    
        oinf = oinf + ',' + str(Ndays) + ',' + ExactDays

# Looking for number of hours tacking exact beginning of the units [iYYYY][iMM][iDD][iHH]0000, [eYYYY][eMM][eDD][iHH+1]0000
##
    if totunits == 'hour':
        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),0,0)
        if edate[10:15] != '0000':
            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS)+1,0,0)
        else:
            etime = edateT
    
        Nhours = 1
        ExactHours = itime.strftime("%Y%m%d%H%M%S")
        it = itime
        while it  + dt.timedelta(hours=1) <= etime:
            it = it + dt.timedelta(hours=1)
            Nhours = Nhours + 1
            ExactHours = ExactHours + '@' + it.strftime("%Y%m%d%H%M%S")
    
        oinf = oinf + ',' + str(Nhours) + ',' + ExactHours

# Looking for number of minutes tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI]00, [eYYYY][eMM][eDD][iHH][iMI+1]00
##
    if totunits == 'minute':
        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),0)
        if edate[12:15] != '00':
            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS)+1,0)
        else:
            etime = edateT

        Nminutes = 1
        ExactMinutes = itime.strftime("%Y%m%d%H%M%S")
        it = itime
        while it + dt.timedelta(minutes=1) <= etime:
            it = it + dt.timedelta(minutes=1)
            Nminutes = Nminutes + 1
            ExactMinutes = ExactMinutes + '@' + it.strftime("%Y%m%d%H%M%S")
    
        oinf = oinf + ',' + str(Nminutes) + ',' + ExactMinutes

# Looking for number of seconds tacking exact beginning of the units [iYYYY][iMM][iDD][iHH][iMI][iSE], 
#   [eYYYY][eMM][eDD][iHH][iMI][iSE+1]
##
    if totunits == 'second':
        itime = dt.datetime(int(iyrS),int(imoS),int(idaS),int(ihoS),int(imiS),int(iseS))
        print 'Lluis ', edate[12:15], '00'
        if edate[12:15] != '00':
            etime = dt.datetime(int(eyrS), int(emoS), int(edaS), int(ehoS), int(emiS), int(eseS)+1)
        else:
            etime = edateT

        Nseconds = 1
        ExactSeconds = itime.strftime("%Y%m%d%H%M%S")
        it = itime
        while it + dt.timedelta(seconds=1) < etime:
            it = it + dt.timedelta(seconds=1)
            Nseconds = Nseconds + 1
            ExactSeconds = ExactSeconds + '@' + it.strftime("%Y%m%d%H%M%S")
    
        oinf = oinf + ',' + str(Nseconds) + ',' + ExactSeconds

    return oinf

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':
        varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150.,         \
          'upward|sensible|heat|flux|at|the|surface', '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 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 index_mat(mat,val):
    """ Function to provide the coordinates of a given value inside a matrix
    index_mat(mat,val)
      mat= matrix with values
      val= value to search
    >>> index_mat(np.arange(27).reshape(3,3,3),22)
    [2 1 1]
    """

    fname = 'index_mat'

    matshape = mat.shape

    matlist = list(mat.flatten())
    ifound = matlist.index(val)

    Ndims = len(matshape)
    valpos = np.zeros((Ndims), dtype=int)
    baseprevdims = np.zeros((Ndims), dtype=int)

    for dimid in range(Ndims):
        baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
        if dimid == 0:
            alreadyplaced = 0
        else:
            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])

    return valpos

def multi_index_mat(mat,val):
    """ Function to provide the multiple coordinates of a given value inside a matrix
    index_mat(mat,val)
      mat= matrix with values
      val= value to search
    >>> vals = np.ones((24), dtype=np.float).reshape(2,3,4)
    vals[:,:,2] = 0.
    vals[1,:,:] = np.pi
    vals[:,2,:] = -1.
    multi_index_mat(vals,1.)
    [array([0, 0, 0]), array([0, 0, 1]), array([0, 0, 3]), array([0, 1, 0]), array([0, 1, 1]), array([0, 1, 3])]
    """
    fname = 'multi_index_mat'

    matshape = mat.shape

    ivalpos = []
    matlist = list(mat.flatten())
    Lmatlist = len(matlist)
    
    val0 = val - val
    if val != val0:
        valdiff = val0
    else:
        valdiff = np.ones((1), dtype = type(val))
    
    ifound = 0
    while ifound < Lmatlist:
        if matlist.count(val) == 0:
            ifound = Lmatlist + 1
        else:
            ifound = matlist.index(val)

            Ndims = len(matshape)
            valpos = np.zeros((Ndims), dtype=int)
            baseprevdims = np.zeros((Ndims), dtype=int)

            for dimid in range(Ndims):
                baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
                if dimid == 0:
                    alreadyplaced = 0
                else:
                    alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
                valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
            matlist[ifound] = valdiff
            ivalpos.append(valpos)

    return ivalpos

def addfileInfile(origfile,destfile,addfile,addsign):
    """ Function to add the content of a file [addfile] to another one [origfile] at 
      a given position [addsign] creating a new file called [destfile]
      origfile= original file
      destfile= destination file
      addfile= content of the file to add
      addsign= sign where to add the file
    """
    fname = 'addfileInfile'

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

    if not os.path.isfile(addfile):
        print errormsg
        print '  ' + fname + ": adding file '" + addfile + "' does not exist !!"
        quit()
    
    ofile = open(origfile, 'r')

# Inspecting flag
##
    Nfound = 0
    for line in ofile:
        if line == addsign + '\n': Nfound = Nfound + 1

    if Nfound == 0:
        print errormsg
        print '  ' + fname + ": adding sign '" + addsign + "' not found !!"
        quit()
    print '  '+ fname + ': found', Nfound," times sign '" + addsign + "' "

    ofile.seek(0)
    dfile = open(destfile, 'w')

    for line in ofile:
        if line != addsign + '\n':
            dfile.write(line)
        else:
            afile = open(addfile,'r')
            for aline in afile:
                print aline
                dfile.write(aline)

            afile.close()

    ofile.close()
    dfile.close()

    print main + " successful writting of file '" + destfile + "' !!"
    return


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

def valmodoper(varVal, valuesS):
    """ Function to run the modification of a variable
      varVAl= matrix with the values
      valuesS= [modins],[[modval1],...,[modvalN]] modification instruction, value with which modify
        [modins]: Modifiers:
          sumc: add [modval1]
          subc: remove [modval1]
          mulc: multiply by [modval1]
          divc: divide by [modval1]
          lowthres: if [val] < [modval1]; val = [modval2]
          upthres: if [val] > [modval1]; val = [modval2]
          lowthres@oper: if [val] < [modval1]; val = [oper] (operation as [modval2],[modval3])
          upthres@oper: if [val] > [modval1]; val = [oper] (operation as [modval2],[modval3])
          potc: [val] ** [modval1]
    """
    fname='valmodoper'

    if valuesS == 'h':
        print fname + '_____________________________________________________________'
        print valmodoper.__doc__
        quit()

    valsS = valuesS.split(',')
    modins = valsS[0]
    modval = float(valsS[1])

    if modins == 'sumc':
        varVal[:] = varVal[:] + modval
    elif modins == 'subc':
        varVal[:] = varVal[:] - modval
    elif modins == 'mulc':
        varVal[:] = varVal[:] * modval
    elif modins == 'divc':
        varVal[:] = varVal[:] / modval
    elif modins == 'lowthres':
        varVal2 = np.where(varVal[:] < float(valsS[1]), float(valsS[2]), varVal[:])
        varVal[:] = varVal2
    elif modins == 'upthres':
        varVal2 = np.where(varVal[:] > float(valsS[1]), float(valsS[2]), varVal[:])
        varVal[:] = varVal2
    elif modins == 'lowthres@oper':
        if valsS[2] == 'sumc': 
           varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
             varVal[:] + float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'subc': 
           varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
             varVal[:] - float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'mulc': 
           varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
             varVal[:] * float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'divc': 
           varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
             varVal[:] / float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'potc': 
           varVal2 = np.where(varVal[:] < float(valsS[1]),                           \
             varVal[:] ** float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        else:
            print errormsg
            print '  ' + fname + ": Operation to modify values '" + modins +         \
              "' is not defined !!"
            quit(-1)
    elif modins == 'upthres@oper':
        if valsS[2] == 'sumc': 
           varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
             varVal[:] + float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'subc': 
           varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
             varVal[:] - float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'mulc': 
           varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
             varVal[:] * float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'divc': 
           varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
             varVal[:] / float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        elif valsS[2] == 'potc': 
           varVal2 = np.where(varVal[:] > float(valsS[1]),                           \
             varVal[:] ** float(valsS[3]), varVal[:])
           varVal[:] = varVal2
        else:
            print errormsg
            print '  ' + fname + ": Operation to modify values '" + modins +         \
              "' is not defined !!"
            quit(-1)
    elif modins == 'potc':
        varVal[:] = varVal[:] ** modval
    else: 
        print errormsg
        print '  ' + fname + ": Operation to modify values '" + modins +             \
          "' is not defined !!"
        quit(-1)

    return varVal

def rangedim(end, shape):
    """Gives the instruction to retrieve values from a dimension of a variable
    >>> print rangedim(-1, 15)
    15
    """
    if end == -1:
      return shape
    else:
      return end

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
    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 timeStr_time(StringDate):
  """ Function to transform a string date ([HH]:[MI]:[SS] format) to a time object
  >>> datetimeStr_datetime('04:32:54')
  04:32:54
  """
  import datetime as dt

  timeT = StringDate.split(':')
  if len(timeT) == 3:
    newtime = dt.time(int(timeT[0]), int(timeT[1]), int(timeT[2]))
  else:
    newtime = dt.time(int(timeT[0]), int(timeT[1]), 0)

  return newtime

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 timeref_datetime_mat(refd, timeval, tu):
    """ Function to transform from a [timeval] in [tu] units from the time referece [tref] to matrix with: year, day, month, hour, minute, second
    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    2   17    8   36   21]
    """
    import datetime as dt
    import numpy as np


    realdates = np.zeros(6, dtype=int)
## 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 tunits == 'seconds':
        realdate = refd + dt.timedelta(seconds=float(timeval))
    elif tunits == 'milliseconds':
        realdate = refd + dt.timedelta(milliseconds=float(timeval))
    else:
          print errormsg
          print '    timeref_datetime: time units "' + tu + '" not ready!!!!'
          quit(-1)

    realdates[0] = int(realdate.year)
    realdates[1] = int(realdate.month)
    realdates[2] = int(realdate.day)
    realdates[3] = int(realdate.hour)
    realdates[4] = int(realdate.second)
    realdates[5] = int(realdate.minute)

    return realdates

def realdatetime_CFcompilant(times, Srefdate, tunits):
    """ Function to transform a matrix with real time values ([year, month, day, hour, minute, second]) to a netCDF one
    times= matrix with times
    Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
    tunits= units of time respect to Srefdate
    >>> realdatetime_CFcompilant(np.array([ [1976, 2, 17, 8, 20, 0], [1976, 2, 18, 8, 20, 0]], dtype=int), '19491201000000', 'hours')
    [ 229784.33333333  229808.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)

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

    return cfdates

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(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    else:
        refdate = dateStr_date(Srefdate)

    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 date_juliandate(refyr, TOTsecs):
    """ Function to transform from a total quantity of seconds since the 
      beginning of a year to 360 days / year calendar
      TOTsecs= total number of seconds
    """
    fname = 'date_juliandate'

    secsYear = 3600*24*30*12.
    secsMonth = 3600*24*30.
    secsDay = 3600*24.
    secsHour = 3600.
    secsMinute = 60.

    rldate = np.zeros((6), dtype=int)
    rldate[0] = refyr 
    if TOTsecs > secsYear:
        rldate[0] = refyr + int(TOTsecs/secsYear) + 1
        TOTsecs = TOTsecs - int(TOTsecs/secsYear)*secsYear

    if np.mod(TOTsecs,secsMonth) == 0:
        rldate[1] = int(TOTsecs/secsMonth)
    else:
        rldate[1] = int(TOTsecs/secsMonth) + 1
    TOTsecs = TOTsecs - int(TOTsecs/secsMonth)*secsMonth

    if np.mod(TOTsecs,secsDay) == 0:
        rldate[2] = int(TOTsecs/secsDay)
    else:
        rldate[2] = int(TOTsecs/secsDay) + 1
    TOTsecs = TOTsecs - int(TOTsecs/secsDay)*secsDay

    if np.mod(TOTsecs,secsHour) == 0:
        rldate[3] = int(TOTsecs/secsHour)
    else:
        rldate[3] = int(TOTsecs/secsHour) + 1
    TOTsecs = TOTsecs - int(TOTsecs/secsHour)*secsHour

    if np.mod(TOTsecs,secsMinute) == 0:
        rldate[4] = int(TOTsecs/secsMinute)
    else:
        rldate[4] = int(TOTsecs/secsMinute) + 1
    TOTsecs = TOTsecs - int(TOTsecs/secsMinute)*secsMinute

    rldate[5] = TOTsecs

#    print refyr,TOTsecs,':',rldate
#    quit()

    return rldate

def CFtimes_datetime(ncfile, tname):
    """ Provide date/time array from a file with a series of netCDF CF-compilant time variable
    ncfile = netCDF file name
    tname = name of the variable time in [ncfile]
    output:
      array(dimt, 0) = year
      array(dimt, 1) = month
      array(dimt, 2) = day
      array(dimt, 3) = hour
      array(dimt, 4) = minute
      array(dimt, 5) = second
    """
    import datetime as dt
    fname = 'CFtimes_datetime'

    times = ncfile.variables[tname]
    timevals = times[:]

    attvar = times.ncattrs()
    if not searchInlist(attvar, 'units'):
        print errormsg
        print '  ' + fname + ": '" + tname + "' does not have attribute: 'units'"
        quit(-1)
    else:
        units = times.getncattr('units')
  
    txtunits = units.split(' ')
    tunits = txtunits[0]
    Srefdate = txtunits[len(txtunits) - 1]
# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
##
    timeval = Srefdate.find(':')

    if not timeval == -1:
#        print '  refdate with time!'
        refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    else:
        refdate = dateStr_date(Srefdate)

    dimt = len(timevals)
    realdates = np.zeros((dimt, 6), dtype=int)

    secsDay=3600*24.

# Checking calendar! 
## 
    y360 = False
    daycal360 = ['earth_360d', '360d', '360days', '360_day']
    if searchInlist(attvar, 'calendar'):
        calendar = times.getncattr('calendar')
        if searchInlist(daycal360,calendar):
            print warnmsg
            print '  ' + fname + ': calendar of 12 months of 30 days !!'
            y360 = True

## 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)
    if y360:
        if tunits == 'weeks':
            for it in range(dimt):
                deltat = dt.timedelta(weeks=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        elif tunits == 'days':
            for it in range(dimt):
                deltat = dt.timedelta(days=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        elif tunits == 'hours':
           for it in range(dimt):
                realdate = dt.timedelta(hours=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        elif tunits == 'minutes':
           for it in range(dimt):
                realdate = dt.timedelta(minutes=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        elif tunits == 'seconds':
           for it in range(dimt):
                realdate = dt.timedelta(seconds=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        elif tunits == 'miliseconds':
           for it in range(dimt):
                realdate = dt.timedelta(miliseconds=float(times[it]))
                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
        else:
              print errormsg
              print '    CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
              quit(-1)
    else:
        if tunits == 'weeks':
            for it in range(dimt):
                realdate = refdate + dt.timedelta(weeks=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        elif tunits == 'days':
            for it in range(dimt):
                realdate = refdate + dt.timedelta(days=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        elif tunits == 'hours':
           for it in range(dimt):
                realdate = refdate + dt.timedelta(hours=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        elif tunits == 'minutes':
           for it in range(dimt):
                realdate = refdate + dt.timedelta(minutes=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        elif tunits == 'seconds':
           for it in range(dimt):
                realdate = refdate + dt.timedelta(seconds=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        elif tunits == 'milliseconds':
           for it in range(dimt):
                realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
                realdates[it,0] = int(realdate.year)
                realdates[it,1] = int(realdate.month)
                realdates[it,2] = int(realdate.day)
                realdates[it,3] = int(realdate.hour)
                realdates[it,4] = int(realdate.second)
                realdates[it,5] = int(realdate.minute)
        else:
              print errormsg
              print '    CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
              quit(-1)
    
    return realdates

class statsValWeigthed(object):
  """Weigthed Statistics class providing:
  vals = values (can be a matrix)
  wgs = weights (can be a matrix)
  self.meanv: mean weigthed value
  self.mean2v: mean quadratic weigthed value
  self.stdv: weigthed standard deviation
  self.Nokvalue non None values of a list of values 
  self.meanwgt: mean of the weigths
  self.mean2wgt: cuadratic mean of the weigths
  self.stdwgt: standard deviation of the weigths
  """

  def __init__(self, vals, wgs):
    if vals is None:
      self.Nv = None
      self.meanv = None
      self.mean2v = None
      self.stdv = None
      self.Nokvalues = None
      self.meanwgt = None
      self.mean2wgt = None
      self.stdwgt = None
    else:
      values = vals.flatten() 
      weights = wgs.flatten()
      self.Nv=len(values)
      self.meanv=0.
      self.mean2v=0.
      self.stdv=0.
      self.meanwgt = 0.
      self.mean2wgt = 0.
      self.stdwgt = 0.
      self.Nokvalues = 0

      for inum in range(self.Nv):
        if not values[inum] is None:
          self.Nokvalues = self.Nokvalues + 1
          self.meanv = self.meanv+values[inum]*weights[inum]
          self.mean2v = self.mean2v+values[inum]*weights[inum]*values[inum]
          self.meanwgt = self.meanwgt+weights[inum]
          self.mean2wgt = self.mean2wgt+weights[inum]*weights[inum]

      self.meanv = self.meanv/float(self.meanwgt)
      self.mean2v = self.mean2v/float(self.meanwgt)
      self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
      self.meanwgt = self.meanwgt/float(self.Nokvalues)
      self.mean2wgt = self.mean2wgt/float(self.Nokvalues)
      self.stdwgt = np.sqrt(self.mean2wgt-self.meanwgt*self.meanwgt)

    return

class statsValWeighted_missVal(object):
  """Weighted Statistics taking into account a missing value class providing:
  vals = values (can be a matrix)
  wgs = weights (can be a matrix)
  missVal= missing value
  self.meanv: mean weigthed value
  self.mean2v: mean quadratic weigthed value
  self.stdv: weigthed standard deviation
  self.Nokvalue non None values of a list of values 
  self.meanwgt: mean of the weigths
  self.mean2wgt: cuadratic mean of the weigths
  self.stdwgt: standard deviation of the weigths
  self.quantilesv: quantiles of the weighted values
  """

  def __init__(self, vals, wgs, missVal):

    fname='statsValWeigthed_missVal'
    if vals is None:
      self.Nv = None
      self.meanv = None
      self.mean2v = None
      self.stdv = None
      self.Nokvalues = None
      self.meanwgt = None
      self.mean2wgt = None
      self.stdwgt = None
      self.quantilesv = None
    else:    
      Npt = 1
      for idim in range(len(vals.shape)):
        Npt = Npt * vals.shape[idim]
      if np.sum(vals >= missVal) == Npt:
          print errormsg
          print '  ' + fname + ' all values get missing!!'
          print errormsg
          quit(-1)
      vals0 = np.where(vals >= missVal, None, vals)
      values = vals0.flatten() 
      weights = wgs.flatten()
      self.Nv=Npt
      self.meanv=0.
      self.mean2v=0.
      self.stdv=0.
      self.meanwgt = 0.
      self.mean2wgt = 0.
      self.stdwgt = 0.
      self.Nokvalues = 0

      valswgt = values
      for inum in range(self.Nv):
        if not values[inum] is None:
          self.Nokvalues = self.Nokvalues + 1
          valswgt[inum] = valswgt[inum]*weights[inum]
          self.meanv = self.meanv+values[inum]*weights[inum]
          self.mean2v = self.mean2v+values[inum]*weights[inum]*values[inum]
          self.meanwgt = self.meanwgt+weights[inum]
          self.mean2wgt = self.mean2wgt+weights[inum]*weights[inum]

      self.meanv = self.meanv/float(self.meanwgt)
      self.mean2v = self.mean2v/float(self.meanwgt)
      self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
      valswgt = valswgt/np.float(self.meanwgt)
      self.meanwgt = self.meanwgt/float(self.Nokvalues)
      self.mean2wgt = self.mean2wgt/float(self.Nokvalues)
      self.stdwgt = np.sqrt(self.mean2wgt-self.meanwgt*self.meanwgt)
      valsq=Quantiles(valswgt, 20)
      self.quantilesv=valsq.quantilesv

    return

class stats2Val(object):
  """two variables Statistics class providing:
  vals1 = variable 1
  vals2 = variable 2
  power = power of the polynomial fitting to apply between both variables
  self.min[var], self.max[var], self.mean[var], self.mean2[var], self.std[var] of 
    [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
  self.Nokvalues1: number of correct values of variable 1
  self.Nokvalues2: number of correct values of variable 2
  self.Nokvalues12: number of correct coincident values of variable 1 and variable 2
  self.mae=mean(abs(var1-var2)) 
  self.rmse=sqrt((var1-var2)**2) 
  self.correlation (and p-value) 
  self.linRegress: linear regression [trend, intercept, regression coefficient, p_value, standard error]
  self.polRegress: polinomial Regresion  of degree [power] [coef**[power], coef**[power-1], ...., coef**0]
  """

  def __init__(self, vals1, vals2, power):
    import numpy as np
    from scipy import stats as sts

    if vals1 is None:
      self.Nv = None
      self.Nokvalues1 = None
      self.Nokvalues2 = None
      self.Nokvalues12 = None
      self.NDvalNone = None
      self.minv1Av2 = None
      self.maxv1Av2 = None
      self.meanv1Av2 = None
      self.mean2v1Av2 = None
      self.stdv1Av2 = None
      self.minv1Sv2 = None
      self.maxv1Sv2 = None
      self.meanv1Sv2 = None
      self.mean2v1Sv2 = None
      self.stdv1Sv2 = None
      self.minv1Dv2 = None
      self.maxv1Dv2 = None
      self.meanv1Dv2 = None
      self.mean2v1Dv2 = None
      self.stdv1Dv2 = None
      self.minv1Pv2 = None
      self.maxv1Pv2 = None
      self.meanv1Pv2 = None
      self.mean2v1Pv2 = None
      self.stdv1Pv2 = None
      self.mae = None
      self.rmse = None
      self.corr = None
      self.linRegress = None
      self.polRegress = None
      self.polRegressResidual = None
      self.polRegressRes = None
      self.polRegressSingVal = None
    else:
      values1 = vals1.flatten() 
      values2 = vals2.flatten() 

      if not len(values1) == len(values2):
        print errormsg
        print '    stats2Val: lengths of variables differ!! Lvar1: ', len(values1), ' Lvar2: ',len(values2),' statistics between them can not be computed!'
        quit(-1)

      self.Nv=len(values1)
      self.minv1Av2=10000000000.
      self.maxv1Av2=-self.minv1Av2
      self.meanv1Av2=0.
      self.mean2v1Av2=0.
      self.stdv1Av2=0.
      self.minv1Sv2=self.minv1Av2
      self.maxv1Sv2=-self.minv1Av2
      self.meanv1Sv2=0.
      self.mean2v1Sv2=0.
      self.stdv1Sv2=0.
      self.minv1Dv2=self.minv1Av2
      self.maxv1Dv2=-self.minv1Av2
      self.meanv1Dv2=0.
      self.mean2v1Dv2=0.
      self.stdv1Dv2=0.
      self.minv1Pv2=self.minv1Av2
      self.maxv1Pv2=-self.minv1Av2
      self.meanv1Pv2=0.
      self.mean2v1Pv2=0.
      self.stdv1Pv2=0.
      self.mae = 0.
      self.rmse = 0.
      self.corr = np.array([0., 0.])
      self.linRegress = np.zeros(5, float)
      self.polRegress = np.zeros(power+1, float)
      self.polRegressResidual = 0.
      self.polRegressSingVal = np.zeros(power+1, float)

# v1 [+ / - / / / *] v2
##
      self.Nokvalues1 = 0
      self.Nokvalues2 = 0
      self.Nokvalues12 = 0
      self.NDvalNone = 0
      for inum in range(self.Nv):
        if not values1[inum] is None:
          self.Nokvalues1 = self.Nokvalues1 + 1
        if not values2[inum] is None:
          self.Nokvalues2 = self.Nokvalues2 + 1
        if not values1[inum] is None and not values2[inum] is None:
          self.Nokvalues12 = self.Nokvalues12 + 1
          Aval = values1[inum] + values2[inum]
          Sval = values1[inum] - values2[inum]
          Pval = values1[inum] * values2[inum]
          if np.isinf(values1[inum] / values2[inum]) or np.isnan(values1[inum] / values2[inum]):
            if self.NDvalNone < 1:
               print warnmsg
               print '      stats2Val: val1/val2 inf or Nan!!!!'
            Dval = None
            self.NDvalNone = self.NDvalNone + 1
          else:
            Dval = values1[inum] / values2[inum]

          self.mae = self.mae + abs(Sval)
          self.rmse = self.rmse + Sval**2

          if Aval < self.minv1Av2:
            self.minv1Av2 = Aval
          if Aval > self.maxv1Av2:
            self.maxv1Av2 = Aval
          if Sval < self.minv1Sv2:
            self.minv1Sv2 = Sval
          if Sval > self.maxv1Sv2:
            self.maxv1Sv2 = Sval
          if not Dval is None and Dval < self.minv1Dv2:
            self.minv1Dv2 = Dval
          if not Dval is None and  Dval > self.maxv1Dv2:
            self.maxv1Dv2 = Dval
          if Pval < self.minv1Pv2:
            self.minv1Pv2 = Pval
          if Pval > self.maxv1Pv2:
            self.maxv1Pv2 = Pval

          self.meanv1Av2 = self.meanv1Av2+Aval
          self.mean2v1Av2 = self.mean2v1Av2+Aval*Aval
          self.meanv1Sv2 = self.meanv1Sv2+Sval
          self.mean2v1Sv2 = self.mean2v1Sv2+Sval*Sval
          if not Dval is None:
            self.meanv1Dv2 = self.meanv1Dv2+Dval
            self.mean2v1Dv2 = self.mean2v1Dv2+Dval*Dval
          self.meanv1Pv2 = self.meanv1Pv2+Pval
          self.mean2v1Pv2 = self.mean2v1Pv2+Pval*Pval

##      print 'Nokvalues1: ', self.Nokvalues1, 'Nokvalues2: ', self.Nokvalues2, 'Nokvalues12: ', float(self.Nokvalues12), 'NDvalNone: ',self.NDvalNone
      self.meanv1Av2 = self.meanv1Av2/float(self.Nokvalues12)
      self.mean2v1Av2 = self.mean2v1Av2/float(self.Nokvalues12)
      self.stdv1Av2 = np.sqrt(self.mean2v1Av2-self.meanv1Av2*self.meanv1Av2)
      self.meanv1Sv2 = self.meanv1Sv2/float(self.Nokvalues12)
      self.mean2v1Sv2 = self.mean2v1Sv2/float(self.Nokvalues12)
      self.stdv1Sv2 = np.sqrt(self.mean2v1Sv2-self.meanv1Sv2*self.meanv1Sv2)
      if self.Nokvalues12 - self.NDvalNone == 0:
          self.meanv1Dv2 = None
          self.mean2v1Dv2 = None
          self.stdv1Dv2 = None
          print warnmsg
          print '      stats2Val: all values of val1/val2 are None!'
      else:
          self.meanv1Dv2 = self.meanv1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
          self.mean2v1Dv2 = self.mean2v1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
          self.stdv1Dv2 = np.sqrt(self.mean2v1Dv2-self.meanv1Dv2*self.meanv1Dv2)
      self.meanv1Pv2 = self.meanv1Pv2/float(self.Nokvalues12)
      self.mean2v1Pv2 = self.mean2v1Pv2/float(self.Nokvalues12)
      self.stdv1Pv2 = np.sqrt(self.mean2v1Pv2-self.meanv1Pv2*self.meanv1Pv2)

      self.mae = self.mae/self.Nokvalues12
      self.rmse = np.sqrt(self.rmse/self.Nokvalues12)

      self.corr = sts.pearsonr(values1, values2)

      self.linRegress[0], self.linRegress[1], self.linRegress[2], self.linRegress[3], self.linRegress[4] = sts.linregress(values1, values2)

      polyfitvals=np.polyfit(values1, values2, power, full = True)

      self.polRegress = polyfitvals[0]
      self.polRegressRes = polyfitvals[1]
      self.polRegressSingVal = polyfitvals[3]

    return

class stats2Val_missVal(object):
  """two variables Statistics taking into account a missing value class providing:
  vals1 = variable 1
  vals2 = variable 2
  missVal = missing value
  power = power of the polynomial fitting to apply between both variables
  self.min[var], self.max[var], self.mean[var], self.mean2[var], self.std[var] of 
    [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]
  self.Nokvalues1: number of correct values of variable 1
  self.Nokvalues2: number of correct values of variable 2
  self.Nokvalues12: number of correct coincident values of variable 1 and variable 2
  self.mae=mean(abs(var1-var2)) 
  self.rmse=sqrt((var1-var2)**2) 
  self.correlation (and p-value) 
  self.linRegress: linear regression [trend, intercept, regression coefficient, p_value, standard error]
  self.polRegress: polinomial Regresion  of degree [power] [coef**[power], coef**[power-1], ...., coef**0]
  """

  def __init__(self, vals1, vals2, power, missVal):
    import numpy as np
    from scipy import stats as sts

    fname='stats2Val_missVal'
    if vals1 is None:
      self.Nv = None
      self.Nokvalues1 = None
      self.Nokvalues2 = None
      self.Nokvalues12 = None
      self.NDvalNone = None
      self.minv1Av2 = None
      self.maxv1Av2 = None
      self.meanv1Av2 = None
      self.mean2v1Av2 = None
      self.stdv1Av2 = None
      self.minv1Sv2 = None
      self.maxv1Sv2 = None
      self.meanv1Sv2 = None
      self.mean2v1Sv2 = None
      self.stdv1Sv2 = None
      self.minv1Dv2 = None
      self.maxv1Dv2 = None
      self.meanv1Dv2 = None
      self.mean2v1Dv2 = None
      self.stdv1Dv2 = None
      self.minv1Pv2 = None
      self.maxv1Pv2 = None
      self.meanv1Pv2 = None
      self.mean2v1Pv2 = None
      self.stdv1Pv2 = None
      self.mae = None
      self.rmse = None
      self.corr = None
      self.linRegress = None
      self.polRegress = None
      self.polRegressResidual = None
      self.polRegressRes = None
      self.polRegressSingVal = None
    else:
      Npt1 = 1
      for idim in range(len(vals1.shape)):
        Npt1 = Npt1 * vals1.shape[idim]
      if np.sum(vals1 >= missVal) == Npt1:
          print errormsg
          print '  ' + fname + ' all values 1 get missing!!'
          print errormsg
          quit(-1)
      Npt2 = 1
      for idim in range(len(vals2.shape)):
        Npt2 = Npt2 * vals2.shape[idim]
      if np.sum(vals2 >= missVal) == Npt2:
          print errormsg
          print '  ' + fname + ' all values 2 get missing!!'
          print errormsg
          quit(-1)
      vals10 = np.where(abs(vals1) >= missVal, None, vals1)
      vals20 = np.where(abs(vals2) >= missVal, None, vals2)
      values1 = vals10.flatten() 
      values2 = vals20.flatten() 

      if not len(values1) == len(values2):
        print errormsg
        print '    stats2Val: lengths of variables differ!! Lvar1: ', len(values1), ' Lvar2: ',len(values2),' statistics between them can not be computed!'
        quit(-1)

      self.Nv=Npt1
      self.minv1Av2=10000000000.
      self.maxv1Av2=-self.minv1Av2
      self.meanv1Av2=0.
      self.mean2v1Av2=0.
      self.stdv1Av2=0.
      self.minv1Sv2=self.minv1Av2
      self.maxv1Sv2=-self.minv1Av2
      self.meanv1Sv2=0.
      self.mean2v1Sv2=0.
      self.stdv1Sv2=0.
      self.minv1Dv2=self.minv1Av2
      self.maxv1Dv2=-self.minv1Av2
      self.meanv1Dv2=0.
      self.mean2v1Dv2=0.
      self.stdv1Dv2=0.
      self.minv1Pv2=self.minv1Av2
      self.maxv1Pv2=-self.minv1Av2
      self.meanv1Pv2=0.
      self.mean2v1Pv2=0.
      self.stdv1Pv2=0.
      self.mae = 0.
      self.rmse = 0.
      self.corr = np.array([0., 0.])
      self.linRegress = np.zeros(5, float)
      self.polRegress = np.zeros(power+1, float)
      self.polRegressResidual = 0.
      self.polRegressSingVal = np.zeros(power+1, float)

# v1 [+ / - / / / *] v2
##
      self.Nokvalues1 = 0
      self.Nokvalues2 = 0
      self.Nokvalues12 = 0
      self.NDvalNone = 0
      for inum in range(self.Nv):
        if not values1[inum] is None:
          self.Nokvalues1 = self.Nokvalues1 + 1
        if not values2[inum] is None:
          self.Nokvalues2 = self.Nokvalues2 + 1
        if not values1[inum] is None and not values2[inum] is None:
          self.Nokvalues12 = self.Nokvalues12 + 1
          Aval = values1[inum] + values2[inum]
          Sval = values1[inum] - values2[inum]
          Pval = values1[inum] * values2[inum]
          if np.isinf(values1[inum] / values2[inum]) or np.isnan(values1[inum] / values2[inum]):
            if self.NDvalNone < 1:
               print warnmsg
               print '      stats2Val: val1/val2 inf or Nan!!!!'
            Dval = None
            self.NDvalNone = self.NDvalNone + 1
          else:
            Dval = values1[inum] / values2[inum]

          self.mae = self.mae + abs(Sval)
          self.rmse = self.rmse + Sval**2

          if Aval < self.minv1Av2:
            self.minv1Av2 = Aval
          if Aval > self.maxv1Av2:
            self.maxv1Av2 = Aval
          if Sval < self.minv1Sv2:
            self.minv1Sv2 = Sval
          if Sval > self.maxv1Sv2:
            self.maxv1Sv2 = Sval
          if not Dval is None and Dval < self.minv1Dv2:
            self.minv1Dv2 = Dval
          if not Dval is None and  Dval > self.maxv1Dv2:
            self.maxv1Dv2 = Dval
          if Pval < self.minv1Pv2:
            self.minv1Pv2 = Pval
          if Pval > self.maxv1Pv2:
            self.maxv1Pv2 = Pval

          self.meanv1Av2 = self.meanv1Av2+Aval
          self.mean2v1Av2 = self.mean2v1Av2+Aval*Aval
          self.meanv1Sv2 = self.meanv1Sv2+Sval
          self.mean2v1Sv2 = self.mean2v1Sv2+Sval*Sval
          if not Dval is None:
            self.meanv1Dv2 = self.meanv1Dv2+Dval
            self.mean2v1Dv2 = self.mean2v1Dv2+Dval*Dval
          self.meanv1Pv2 = self.meanv1Pv2+Pval
          self.mean2v1Pv2 = self.mean2v1Pv2+Pval*Pval

##      print 'Nokvalues1: ', self.Nokvalues1, 'Nokvalues2: ', self.Nokvalues2, 'Nokvalues12: ', float(self.Nokvalues12), 'NDvalNone: ',self.NDvalNone
      self.meanv1Av2 = self.meanv1Av2/float(self.Nokvalues12)
      self.mean2v1Av2 = self.mean2v1Av2/float(self.Nokvalues12)
      self.stdv1Av2 = np.sqrt(self.mean2v1Av2-self.meanv1Av2*self.meanv1Av2)
      self.meanv1Sv2 = self.meanv1Sv2/float(self.Nokvalues12)
      self.mean2v1Sv2 = self.mean2v1Sv2/float(self.Nokvalues12)
      self.stdv1Sv2 = np.sqrt(self.mean2v1Sv2-self.meanv1Sv2*self.meanv1Sv2)
      if self.Nokvalues12 - self.NDvalNone == 0:
          self.meanv1Dv2 = None
          self.mean2v1Dv2 = None
          self.stdv1Dv2 = None
          print warnmsg
          print '      stats2Val: all values of val1/val2 are None!'
      else:
          self.meanv1Dv2 = self.meanv1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
          self.mean2v1Dv2 = self.mean2v1Dv2/(float(self.Nokvalues12 - self.NDvalNone))
          self.stdv1Dv2 = np.sqrt(self.mean2v1Dv2-self.meanv1Dv2*self.meanv1Dv2)
      self.meanv1Pv2 = self.meanv1Pv2/float(self.Nokvalues12)
      self.mean2v1Pv2 = self.mean2v1Pv2/float(self.Nokvalues12)
      self.stdv1Pv2 = np.sqrt(self.mean2v1Pv2-self.meanv1Pv2*self.meanv1Pv2)

      self.mae = self.mae/self.Nokvalues12
      self.rmse = np.sqrt(self.rmse/self.Nokvalues12)

      vals1Nomiss = np.ones(len(values1), dtype=bool)
      vals2Nomiss = np.ones(len(values1), dtype=bool)

      for i in range(len(values1)):
          if values1[i] is None:
              vals1Nomiss[i] = False

      for i in range(len(values2)):
          if values2[i] is None:
              vals2Nomiss[i] = False

      v1 = np.array(values1[vals1Nomiss], dtype=float)
      v2 = np.array(values2[vals2Nomiss], dtype=float)

      if not v1.shape == v2.shape:
          print errormsg
          print '  ' + fname + ': variables without missing values v1: ',v1.shape , ' and v2: ',v2.shape ,' do not have the same shape! '
          print errormsg
          quit(-1)

      self.corr = sts.pearsonr(v1, v2)

      self.linRegress[0], self.linRegress[1], self.linRegress[2], self.linRegress[3], self.linRegress[4] = sts.linregress(v1, v2)

      polyfitvals=np.polyfit(v1, v2, power, full = True)

      self.polRegress = polyfitvals[0]
      self.polRegressRes = polyfitvals[1]
      self.polRegressSingVal = polyfitvals[3]

    return

def mask_2masked(vals1, vals2):
    """ Function to provide the boolean matrix (in opposite way as it is in the mask) as combination of mask from to masked matrices
    """
    import numpy.ma as ma
    fname = 'mask_2masked'

#    if len(vals1.shape) != len(vals2.shape):
#        print errormsg
#        print '  ' + fname + ' matrix 1 :', len(vals1.shape), ' and matrix 2 ', len(vals2.shape), ' have different size!'
#        print errormsg
#        quit(-1)

#    for idim in range(len(vals1.shape)):
#        if vals1.shape[idim] != vals2.shape[idim]:
#            print errormsg
#            print '  ' + fname + ' dimension ', idim,' from matrix 1 :', vals1.shape[idim], ' and matrix 2 ', \
#                vals2.shape[idim], ' have different size!'
#            print errormsg
#            quit(-1)

    if type(vals1) == type(ma.array(1)):
        mask1array=np.where(ma.getmaskarray(vals1) == False, True, False)
    else:
        mask1array=np.ones(vals1.shape, dtype=bool)

    if type(vals2) == type(ma.array(1)):
        mask2array=np.where(ma.getmaskarray(vals2) == False, True, False)
    else:
        mask2array=np.ones(vals2.shape, dtype=bool)

    mask12 = mask1array*mask2array

    return mask12

def mask_pearsonr(xvals, yvals):
    """ Function to compute a pearson correlation from mask matrices
    """
    from scipy import stats as sts
    fillVal = 1.e20
    maskxy = mask_2masked(xvals, yvals)

    if np.sum(maskxy) > 1:
        pearsonr = sts.pearsonr(xvals[maskxy], yvals[maskxy])
        if np.isnan(pearsonr[0]) or np.isnan(pearsonr[1]): 
           pearsonr = ( fillVal, fillVal)
    else:
        pearsonr = (fillVal, fillVal)

    return pearsonr

def mask_quantiles(maskmat, Nquants):
    """ Function to provide the quantiles of a masked array 20 for %5 bins (21 in total)
    """
    import numpy.ma as ma

    fillValue = 1.e20

    sortmat = maskmat.flatten().copy()
    sortmat.sort()
    quants = np.zeros(Nquants+1, dtype=type(maskmat[0]))
    Nv = ma.size(maskmat)
    NoMask=maskmat.count()

    if NoMask < Nquants:
        quants[:] = fillValue
    else:
        for iq in range(Nquants):
            quants[iq] = sortmat[int((NoMask-1)*iq/(Nquants))]

        quants[Nquants] = sortmat[NoMask-1]

    return quants

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 mask_linregres(vals1, vals2):
    """ Function to compute a linear regression from masked data
    vals1: x-values for the regresion
    vals2: y-values for the regresion
    """
    import numpy.ma as ma

    fname = 'mask_linregres'

    missval1 = vals1.get_fill_value()    
    missval2 = vals2.get_fill_value()    

    vals10 = np.where(abs(vals1) >= abs(missval1*0.9), None, vals1)
    vals20 = np.where(abs(vals2) >= abs(missval2*0.9), None, vals2)

    values1 = vals10.flatten() 
    values2 = vals20.flatten() 

    vals1Nomiss = np.ones(len(values1), dtype=bool)
    vals2Nomiss = np.ones(len(values2), dtype=bool)

    for i in range(len(values1)):
        if values1[i] is None:
            vals1Nomiss[i] = False
    for i in range(len(values2)):
        if values2[i] is None:
            vals2Nomiss[i] = False

    v1 = np.array(values1[vals1Nomiss], dtype=float)
    v2 = np.array(values2[vals2Nomiss], dtype=float)

    if len(v1) != len(v2):
        print errormsg
        print fname + ': length of masked matrices mat1:',len(v1),'and mat2:',len(v2),'does not match!'
        print errormsg
        quit(-1)

    linregres = np.array(sts.linregress(v1, v2), dtype= np.float64)

    return linregres

def mask_space_stats(maskmat,statsvals,dim):
    """ Function to give back the multi-dimensional statisitcs of a given masked array
    maskmat=multidimensional masked array
    statsvals=[statn]:[values] 
      [statn]: statistics to do: 
        'quant', quantiles
      [values]: value for the statistics: Nquantiles
    dim= dimension to run the statistics
    """
    from sys import stdout

    fname = 'mask_space_stats'

    statn=statsvals.split(':')[0]
    if len(statsvals.split(':')) > 1:
        values=statsvals.split(':')[1]

    maskshape = maskmat.shape
    Ndims = len(maskshape)
    if statn == 'quant':
        if len(statsvals.split(':')) == 1:
            print errormsg
            print fname + ': statistics "' + statn + '" requires a value!!!'
            print errormsg
            quit(-1)
        Nquants=int(values)
        if Ndims == 2:
            if dim == 0:
                statval = np.ones((21, maskshape[1]), dtype=np.float64)*fillValue
                for i in range(maskshape[1]):
                   percendone(i, maskshape[1], 5, 'quantiles')
                   statval[:,i] = mask_quantiles(maskmat[:,i],Nquants)
            if dim == 1:
                statval = np.ones((21, maskshape[0]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                   percendone(i, maskshape[0], 5, 'quantiles')
                   statval[:,i] = mask_quantiles(statval[i,:],Nquants)
        elif Ndims == 3:
            if dim == 0:
                statval = np.ones((21, maskshape[1], maskshape[2]), dtype=np.float64)*fillValue
                for i in range(maskshape[1]):
                    for j in range(maskshape[2]):
                        percendone(i*maskshape[2] + j, maskshape[1]*maskshape[2], 5, 'quantiles')
                        statval[:,i,j] = mask_quantiles(maskmat[:,i,j],Nquants)
            if dim == 1:
                statval = np.ones((21, maskshape[0], maskshape[2]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                    for j in range(maskshape[2]):
                        percendone(i*maskshape[2] + j, maskshape[0]*maskshape[2], 5, 'quantiles')
                        statval[:,i,j] = mask_quantiles(maskmat[i,:,j],Nquants)
            if dim == 2:
                statval = np.ones((21, maskshape[0], maskshape[1]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                    for j in range(maskshape[1]):
                        percendone(i*maskshape[1] + j, maskshape[0]*maskshape[1], 5, 'quantiles')
                        statval[:,i,j] = mask_quantiles(maskmat[i,j,:],Nquants)
        elif Ndims == 4:
            if dim == 0:
                statval = np.ones((21, maskshape[1], maskshape[2], maskshape[3]), dtype=np.float64)*fillValue
                for i in range(maskshape[1]):
                    for j in range(maskshape[2]):
                        for k in range(maskshape[3]):
                            percendone(i*maskshape[1]*maskshape[2] + j*maskshape[2] + k, maskshape[1]*maskshape[2]*maskshape[3], 5, 'quantiles')
                            statval[:,i,j,k] = mask_quantiles(maskmat[:,i,j,k],Nquants)
            if dim == 1:
                statval = np.ones((21, maskshape[0], maskshape[2], maskshape[3]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                    for j in range(maskshape[2]):
                        for k in range(maskshape[3]):
                            percendone(i*maskshape[0]*maskshape[2] + j*maskshape[2] + k, maskshape[0]*maskshape[2]*maskshape[3], 5, 'quantiles')
                            statval[:,i,j,k] = mask_quantiles(maskmat[i,:,j,k],Nquants)
            if dim == 2:
                statval = np.ones((21, maskshape[0], maskshape[1], maskshape[3]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                    for j in range(maskshape[1]):
                        for k in range(maskshape[3]):
                            percendone(i*maskshape[0]*maskshape[1] + j*maskshape[1] + k, maskshape[0]*maskshape[1]*maskshape[3], 5, 'quantiles')
                            statval[:,i,j,k] = mask_quantiles(maskmat[i,j,:,k],Nquants)
            if dim == 3:
                statval = np.ones((21, maskshape[0], maskshape[1], maskshape[2]), dtype=np.float64)*fillValue
                for i in range(maskshape[0]):
                    for j in range(maskshape[1]):
                        for k in range(maskshape[2]):
                            percendone(i*maskshape[0]*maskshape[1] + j*maskshape[1] + k, maskshape[0]*maskshape[1]*maskshape[2], 5, 'quantiles')
                            statval[:,i,j,k] = mask_quantiles(maskmat[i,j,k,:],Nquants)
        else:
            print errormsg
            print fname + ': size of matrix ', Ndims,'not ready!!!'
            print errormsg
            quit(-1)

    else:
        print errormsg
        print fname + ':  statistics "' + statn + '" not ready!!!!'
        print errormsg
        quit(-1)

    print stdout.write("\n")
    
    return statval

class statsVal(object):
  """Statistics class providing
  vals = variable
  self.Nv = number of values
  self.minv = minimum value
  self.maxv = maximum value
  self.meanv = mean value
  self.mean2v = cuadratic mean value
  self.stdv = standard deviation value
  self.Nokvalues = number of correct values of variable
  self.quantilesv = quantiles (%5 bins) of the variable
  """

  def __init__(self, vals):
    if vals is None:
      self.Nv = None
      self.minv = None
      self.maxv = None
      self.meanv = None
      self.mean2v = None
      self.stdv = None
      self.Nokvalues = None
      self.quantilesv = None
    else:
      values = vals.flatten() 
      self.Nv=len(values)
      self.minv=10000000000.
      self.maxv=-100000000.
      self.meanv=0.
      self.mean2v=0.
      self.stdv=0.

      sortedvalues = sorted(values)

      self.Nokvalues = 0
      for inum in range(self.Nv):
        if not values[inum] is None:
          self.Nokvalues = self.Nokvalues + 1
          if values[inum] < self.minv:
            self.minv = values[inum]
          if values[inum] > self.maxv:
            self.maxv = values[inum]

          self.meanv = self.meanv+values[inum]
          self.mean2v = self.mean2v+values[inum]*values[inum]

      self.meanv = self.meanv/float(self.Nokvalues)
      self.mean2v = self.mean2v/float(self.Nokvalues)
      self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
      self.quantilesv = []
      for iq in range(20):
        self.quantilesv.append(sortedvalues[int((self.Nv-1)*iq/20)])

      self.quantilesv.append(sortedvalues[self.Nv-1])
      self.medianv = self.quantilesv[10]

    return

class Quantiles(object):
    """ Class to provide quantiles from a given arrayof values
    """

    def __init__(self, values, Nquants):
        import numpy.ma as ma

        if values is None:
            self.quantilesv = None

        else:
            self.quantilesv = []

            vals0 = values.flatten()
            Nvalues = len(vals0)
            vals = ma.masked_equal(vals0, None)
            Nvals=len(vals.compressed())

            sortedvals = sorted(vals.compressed())
            for iq in range(Nquants):
                self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])

            self.quantilesv.append(sortedvals[Nvals-1])

        return

class statsVal_missVal(object):
  """Statistics class tacking into account a missing value providing
  vals = variable
  missval = missing value
  self.Nv = number of values
  self.minv = minimum value
  self.maxv = maximum value
  self.meanv = mean value
  self.mean2v = cuadratic mean value
  self.stdv = standard deviation value
  self.Nokvalues = number of correct values of variable
  self.quantilesv = quantiles (%5 bins) of the variable
  """

  def __init__(self, vals, missVal):
    fname='statsVal_missVal'
  
    if vals is None:
      self.Nv = None
      self.minv = None
      self.maxv = None
      self.meanv = None
      self.mean2v = None
      self.stdv = None
      self.Nokvalues = None
      self.quantilesv = None
    else:
      Npt = 1
      for idim in range(len(vals.shape)):
        Npt = Npt * vals.shape[idim]
      if np.sum(vals >= missVal) == Npt:
          print errormsg
          print '  ' + fname + ' all values get missing!!'
          print errormsg
          quit(-1)

      vals1 = np.where(abs(vals) >= missVal, None, vals)
      vals0 = np.where(np.isnan(vals1), None, vals1)

      values = vals0.flatten() 
      self.Nv=Npt
      self.minv=10000000000.
      self.maxv=-100000000.
      self.meanv=0.
      self.mean2v=0.
      self.stdv=0.

      sortedvalues = sorted(values)

      self.Nokvalues = 0
      for inum in range(self.Nv):
        if not values[inum] is None:
          self.Nokvalues = self.Nokvalues + 1
          if values[inum] < self.minv:
            self.minv = values[inum]
          if values[inum] > self.maxv:
            self.maxv = values[inum]

          self.meanv = self.meanv+values[inum]
          self.mean2v = self.mean2v+values[inum]*values[inum]

      self.meanv = self.meanv/float(self.Nokvalues)
      self.mean2v = self.mean2v/float(self.Nokvalues)
      self.stdv = np.sqrt(self.mean2v-self.meanv*self.meanv)
      self.quantilesv = []
      for iq in range(20):
        self.quantilesv.append(sortedvalues[int((self.Nv-1)*iq/20)])

      self.quantilesv.append(sortedvalues[self.Nv-1])
      self.medianv = self.quantilesv[10]

    return

def printing_class(classobj):
    """ Function to print all the values of a given class
    """

    valscls = vars(classobj)
    for attrcls in valscls:
        print attrcls, ':', valscls[attrcls]

    return

def fmtprinting_class(classobj):
    """ Function to print all the values of a given class
    """

    valscls = vars(classobj)
    for attrcls in valscls:
        print '@' + attrcls + '@', ':', valscls[attrcls]

    return

def cycl_incr(cyclval,cycle,ninc):
    """ Function to increment a cyclic value [cyclval] with a cycle [cycle] a given number [ninc] of times
    >>> cycl_incr(1,4,1)
    2
    >>> cycl_incr(3,4,1)
    0
    >>> cycl_incr(1,4,10)
    3
    """

    if ninc >= cycle:
        print 'cycl_incr: WARNING -- warning -- WARNING -- warning'
        print '    given increment: ', ninc,' is larger than the cycle !!'
        ninc = ninc - cycle*int(ninc/cycle)
        print '    reducing it to: ', ninc

    val=cyclval + ninc
    if val >= cycle:
        val=cyclval + ninc - cycle

    return val

def times_4seasons(tvals, integrity):
    """ Function to split a time series in matricial date format ([:,year,mon,day,hour,minute,second]) in the four stations DJF,MAM,JJA,SON
    tvals= matrix withe times as [:,year,mon,day,hour,minute,second]
    integrity= only give values for entire seasons [True, 3 months for the season], or just split the values by seasons [False] 
    """
    fillVal=1.e20

#    print tvals

    dt=tvals.shape[0]
    seasons=np.ones((dt,4),dtype=bool)
    seasons=seasons*False

    monseas=[12,3,6,9]
    firstseas=False

    for it in range(dt):
        if not firstseas:
            if integrity:
                for iseas in range(4):
                    if tvals[it,1] == monseas[iseas]:
                        nseas=iseas
                        seasons[it,nseas]=True
                        firstseas=True
                        begseas=it
            else:
                for iseas in range(4):
                    for imon in range(3):
                        if tvals[it,1] == cycl_incr(monseas[iseas],12,imon):
                            nseas=iseas
                            seasons[it,nseas]=True
                            firstseas=True                    
        else:
            newseas=cycl_incr(nseas,4,1)
            if tvals[it,1] == monseas[newseas]:
                seasons[it,newseas] = True
                nseas=newseas
                begseas=it
            else:
                seasons[it,nseas] = True

    endseas = it
##    print 'Last season: ',nseas,' beginnig: ',begseas,' ending: ',endseas

# Checking integrity of last season (has to have 3 months)
##
    if integrity:
        fullseas=True
        Nmon=np.unique(tvals[begseas:endseas+1,1])
        for it in range(begseas,endseas): fullseas=fullseas*seasons[it,nseas]
        if len(Nmon) < 3 or not fullseas:
            seasons[begseas:endseas+1,nseas] = False

    return seasons

def realdatetime_CFcompilant(times, Srefdate, tunits):
    """ Function to transform a matrix with real time values ([year, month, day, hour, minute, second]) to a netCDF one
    times= matrix with times
    Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
    tunits= units of time respect to Srefdate
    >>> realdatetime_CFcompilant(np.array([ [1976, 2, 17, 8, 20, 0], [1976, 2, 18, 8, 20, 0]], dtype=int), '19491201000000', 'hours')
    [ 229784.33333333  229808.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)

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

    return cfdates

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 netCDFdatetime_realdatetime(units, tcalendar, times):
    """ Function to transfrom from netCDF CF-compilant times to real time
    [units]= CF time units [tunits] since [YYYY]-[MM]-[HH] [[HH]:[MI]:[SS]]
    [tcalendar]= time calendar
    [times]= CF time values
    """
    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(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    else:
        refdate = dateStr_date(Srefdate)

    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

class cls_yearly_means(object):
    """ Class to compute the yearly mean of a time series of values
      tvals=time values (year, month, day, hour, minute, second) format
      values=time-series of values
      dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
      tunits= time units
      self.dimt = Number of values
      self.yrv = yearly mean values
      self.yrt = years of the mean values
    >>> timesv = netCDFdatetime_realdatetime('days since 1949-12-01 00:00:00', 'standard', np.arange(1129))
    >>> valuesv = np.zeros((1129), dtype=np.float)
    >>> valuesv[0:31] = 0.
    >>> valuesv[31:394] = 1.
    >>> valuesv[395:761] = 2.
    >>> valuesv[761:1127] = 3.
    >>> valuesv[1127:1129] = 4.

    >>> yrlmeans = cls_yearly_means(timesv, valuesv, '19491201000000', 'days')
    >>> print yrlmeans.dimt,yrlmeans.yrv, yrlmeans.yrt 
    5 [ 0.  1.  2.  3.  4.] [ 1949.  1950.  1951.  1952.  1953.]
    """
    def __init__(self, tvals, values, dateref, tunits):

        import numpy.ma as ma
        fillVal=1.e20
        fname = 'cls_yearly_means'

        if tvals is None:
            self.dimt = None
            self.yrv = None
            self.yrt = None

        elif tvals == 'h':
            print fname + '_____________________________________________________________'
            print cls_yearly_means.__doc__
            quit()

        else:
            dt=len(values)
            if dt > 0:
                years = np.unique(tvals[:,0])
                Nyrs = len(years)
                yrmean = np.zeros((Nyrs), dtype=np.float)
                yrt = np.zeros((Nyrs), dtype=np.float)

                iiyr = tvals[0,0]
                yrmean[0] = values[0]
                yrt[0] = iiyr
                iit = 0
                for iyr in range(Nyrs):
                    nvals = 1
                    for it in range(iit+1,dt):
#                        print iyr, iiyr, it, tvals[it,0],':',values[it],'->',yrmean[iyr]
                        if tvals[it,0] == iiyr:
                            yrmean[iyr] = yrmean[iyr] + values[it]
                            nvals = nvals + 1
                        else:
                            yrmean[iyr] = yrmean[iyr] / (nvals*1.)
                            iiyr = tvals[it,0]
                            yrmean[iyr + 1] = values[it]
                            yrt[iyr+1] = iiyr
                            iit = it
                            break

                yrmean[Nyrs-1] = yrmean[Nyrs-1]/nvals
                self.dimt = Nyrs
                self.yrv = yrmean
                self.yrt = yrt
            else:
                print errormsg
                print '  ' + fname + ': No values passed!'
                print '    values:', values
                quit(-1)

        return

class cls_seasonal_means(object):
    """ Class to compute the seasonal mean of a time series of values
    tvals=time values (year, month, day, hour, minute, second) format
    values=time-series of values
    dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
    tunits= time units
    self.dimt = Number of values
    self.seasv = seasonal mean values
    self.seast = seasonal time mean values
    """
    def __init__(self, tvals, values, dateref, tunits):

        import numpy.ma as ma
        fillVal=1.e20
        fname = 'cls_seasonal_means'

        if tvals is None:
            self.dimt = None
            self.seasv = None
            self.seast = None

        else:
            tseasons=times_4seasons(tvals, True)

            dt=len(values)
            seasvals=np.ones((dt/12,4), dtype=np.float64)*fillVal
            seastimes=np.zeros((dt/12,4), dtype=np.float64)
            dates = realdatetime_CFcompilant(tvals, dateref, tunits)

            for iseas in range(4):
                for it in range(dt):
                    if tseasons[it,iseas]:
                        tfirstseas=int(it/11)
                        firstseas=True
##                        print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
                        break
                for itt in range(it, dt-1,12):
                    seasvals[int(itt/12),iseas]=(values[itt] + values[itt+1] + values[itt+2])/3.
                    seastimes[int(itt/12),iseas]=dates[itt] + (dates[itt+2] - dates[itt])/2.
##                    print itt, values[itt], values[itt+1], values[itt+2], '--->', seasvals[int(itt/12),iseas]

            self.dimt = dt/12
            self.seasv = ma.masked_equal(seasvals, fillVal)
            self.seast = ma.masked_equal(seastimes, fillVal)

        return

class cls_seasonal_accums(object):
    """ Class to compute the seasonal accumulations of a time series of values
    tvals=time values (year, month, day, hour, minute, second) format
    values=time-series of values
    dateref= reference date [YYYY][MM][DD][HH][MI][SS] format
    tunits= time units
    self.dimt = Number of values
    self.seasv = seasonal mean values
    self.seast = seasonal time mean values
    """
    def __init__(self, tvals, values, dateref, tunits):
        import numpy.ma as ma
        fillVal=1.e20

        if tvals is None:
            self.dimt = None
            self.seasv = None
            self.seast = None

        else:
            tseasons=times_4seasons(tvals, True)

            dt=len(values)
            seasvals=np.ones((dt/12,4), dtype=np.float64)*fillVal
            seastimes=np.zeros((dt/12,4), dtype=np.float64)
            dates = realdatetime_CFcompilant(tvals, dateref, tunits)

            for iseas in range(4):
                for it in range(dt):
                    if tseasons[it,iseas]:
                        tfirstseas=int(it/11)
                        firstseas=True
#                        print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
                        break
                for itt in range(it, dt-1,12):
                    seasvals[int(itt/12),iseas]=values[itt] + values[itt+1] + values[itt+2]
                    seastimes[int(itt/12),iseas]=dates[itt] + (dates[itt+2] - dates[itt])/2.

            self.dimt = dt/12
            self.seasv = ma.masked_equal(seasvals, fillVal)
            self.seast = ma.masked_equal(seastimes, fillVal)

        return

def seasonal_means(tvals, values):
    """ Function to compute the seasonal mean of a time series of values
    tvals=time values (year, month, day, hour, minute, second) format
    values=time-series of values
    """
    fillVal=1.e20

    tseasons=times_4seasons(tvals, True)

    dt=len(values)
    seasvals=np.ones((dt/4,4), dtype=np.float64)
    seasvals=seasvals*fillVal

    for iseas in range(4):
        for it in range(dt):
            if tseasons[it,iseas]:
                tfirstseas=int(it/11)
                firstseas=True
#                print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
                break
        for itt in range(it, dt-1,12):
            seasvals[int(itt/12),iseas]=(values[itt] + values[itt+1] + values[itt+2])/3.

    return seasvals

def seasonal_accum(tvals, values):
    """ Function to compute the seasonal accumulation of a time series of values
    tvals=time values (year, month, day, hour, minute, second) format
    values=time-series of values
    """
    fillVal=1.e20

    tseasons=times_4seasons(tvals, True)

    dt=len(values)
    seasvals=np.ones((dt/4,4), dtype=np.float64)
    seasvals=seasvals*fillVal

    for iseas in range(4):
        for it in range(dt):
            if tseasons[it,iseas]:
                tfirstseas=int(it/11)
                firstseas=True
#                print '  first time-step of season ', iseas, ' is: ',it,' position: ',tfirstseas
                break
        for itt in range(it, dt-1,12):
            seasvals[int(itt/11),iseas]=values[itt] + values[itt+1] + values[itt+2]

    return seasvals

def seasmean(timename, filename, varn):
    """ Function to compute the seasonal mean of a variable
    timename= name of the variable time in the file
    filename= netCDF file name
    varn= name of the variable
    """
    fillValue=1.e20

    ncf = NetCDFFile(filename,'a')
    ofile = 'seasmean_' + varn + '.nc'

    if not ncf.variables.has_key(varn):
        print errormsg
        print '    seasmean: File  "' + filename + '" does not have variable ' + varn + ' !!!!'
        print errormsg
        ncf.close()
        quit(-1)

    varobj = ncf.variables[varn]
    varinf = variable_inf(varobj)
    timeobj = ncf.variables[timename]
    timeinf = cls_time_information(ncf, timename)

    timevals=timeobj[:]
    netcdftimes = netCDFdatetime_realdatetime(timeinf.units + ' since ' + timeinf.Urefdate, timeinf.calendar, timevals)
    
    timeseasons=times_4seasons(netcdftimes, True)

    timeseasvals=seasonal_means(netcdftimes, timevals)
    
    ncfnew = NetCDFFile(ofile,'w')
    ncfnew.createDimension('seastime', timeinf.dimt/4)
    ncfnew.createDimension('seas', 4)

    newvarobj = ncfnew.createVariable('seastime', 'f4', ('seastime', 'seas',), fill_value=fillValue)
    newvarobj[:]=timeseasvals
    attr = newvarobj.setncattr('calendar', timeinf.calendar)
    attr = newvarobj.setncattr('units', timeinf.units)

    newvarobj = ncfnew.createVariable('seasn', str, ('seas'))
    attr = newvarobj.setncattr('standard_name', 'seasn')
    attr = newvarobj.setncattr('long_name', 'name of the season')
    attr = newvarobj.setncattr('units', '3-months')
    newvarobj[0]='DJF'
    newvarobj[1]='MAM'
    newvarobj[2]='JJA'
    newvarobj[3]='SON'

    newdims=[]
    idim = 0

    for ndim in varinf.dimns:
        if ndim == timename:
            newdims.append(unicode('seastime'))
            if not idim == 0:
                print errormsg
                print '  seasmean: File  "' + filename + '" does not have time dimension "' + timename + '" as the first one. It has it as # ', idim, ' !!!!'
                print errormsg
                ncf.close()
                quit(-1)
        else:
            newdims.append(ndim)

        if not ncfnew.dimensions.has_key(ndim):
            ncfnew.sync()
            ncfnew.close()
            ncf.close()            

            fdimadd(filename + ',' + ndim, ofile)
            ncf = NetCDFFile(filename,'r')
            ncfnew = NetCDFFile(ofile,'a')

    print ncfnew.dimensions
    namedims=tuple(newdims)
    print namedims

    varobj = ncf.variables[varn]
    varinf = variable_inf(varobj)

    newvarobj=ncfnew.createVariable(varn + '_seasmean', 'f4', namedims, fill_value=fillValue)
    attr = newvarobj.setncattr('standard_name', varn + '_seasmean')
    attr = newvarobj.setncattr('long_name', 'seasonal mean of ' + varn)
    attr = newvarobj.setncattr('units', varinf.units)

    newvar = np.ones(varinf.dims[0], dtype=varinf.dtype)
    if varinf.Ndims == 1:
        newvar = newvarobj[:]
        newvarobj[:] = seasonal_means(netcdftimes, newvar)
    if varinf.Ndims == 2:
        for i in range(varinf.dims[1]):
            newvar = newvarobj[:,i]
            newvarobj[:,i] = seasonal_means(netcdftimes, newvar)
    if varinf.Ndims == 3:
        for i in range(varinf.dims[1]):
            for j in range(varinf.dims[2]):
                newvar = newvarobj[:,i,j]
                newvarobj[:,i,j] = seasonal_means(netcdftimes, newvar)
    if varinf.Ndims == 4:
        for i in range(varinf.dims[1]):
            for j in range(varinf.dims[2]):
                for k in range(varinf.dims[3]):
                    newvar = newvarobj[:,i,j,k]
                    newvarobj[:,i,j,k] = seasonal_means(netcdftimes, newvar)
    if varinf.Ndims == 5:
        for i in range(varinf.dims[1]):
            for j in range(varinf.dims[2]):
                for k in range(varinf.dims[3]):
                    for l in range(varinf.dims[4]):
                        newvar = newvarobj[:,i,j,k,l]
                        newvarobj[:,i,j,k,l] = seasonal_means(netcdftimes, newvar)
    if varinf.Ndims == 6:
        for i in range(varinf.dims[1]):
            for j in range(varinf.dims[2]):
                for k in range(varinf.dims[3]):
                    for l in range(varinf.dims[4]):
                        for m in range(varinf.dims[5]):
                            newvar = newvarobj[:,i,j,k,l,m]
                            newvarobj[:,i,j,k,l,m] = seasonal_means(netcdftimes, newvar)
    else:
        print errormsg
        print '  seasmean: number of dimensions ', varinf.Ndims, ' not ready !!!!'
        print errormsg
        ncf.close()
        quit(-1)

    ncfnew.sync()
    ncfnew.close()

    print 'seasmean: Seasonal mean file "' + ofile + '" written!'

def load_variable_lastdims(var, prevdims, Nlastdims):
    """ Function to load the last [Nlastdims] dimensions of a variable [var] at the other dimensions at [prevdims]
    >>> load_variable_lastdims(np.array(range(5*5*5)).reshape(5,5,5), [1,2], 1)
    [35 36 37 38 39]
    """
    fname='load_variable_lastdims'

    dims=var.shape
    Ndims=len(dims)

    Nprevdims = len(prevdims)
    if not Nprevdims + Nlastdims == Ndims:
        print erromsg
        print '  ' + fname + ': number of dimensions previous (',Nprevdim,') and last (',Nlastdims, ') does not match with variable size: ',Ndims,' !!!!'
        print errormsg
        quit(-1)

    if Nlastdims > Ndims-1:
        print errormsg
        print '  ' + fname + ': number of last dimensions ', Nlastdims,' >= than the number of dimensions of the variable ', Ndims,' !!!!'
        print errormsg
        quit(-1)

    if Ndims == 1:
        loadvar = var[:]
    elif Ndims == 2:
        loadvar = var[prevdims[0], :]
    elif Ndims == 3:
        if Nlastdims == 1:
            loadvar = var[prevdims[0], prevdims[1], :]
        else:
            loadvar = var[prevdims[0], :, :]
    elif Ndims == 4:
        if Nlastdims == 1:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], :]
        elif Nlastdims == 2:
            loadvar = var[prevdims[0], prevdims[1], :, :]
        else:
            loadvar = var[prevdims[0], :, :, :]
    elif Ndims == 5:
        if Nlastdims == 1:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], :]
        elif Nlastdims == 2:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], :, :]
        elif Nlastdims == 3:
            loadvar = var[prevdims[0], prevdims[1], :, :, :]
        else:
            loadvar = var[prevdims[0], :, :, :, :]
    elif Ndims == 6:
        if Nlastdims == 1:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], prevdims[4], :]
        elif Nlastdims == 2:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], prevdims[3], :, :]
        elif Nlastdims == 3:
            loadvar = var[prevdims[0], prevdims[1], prevdims[2], :, :, :]
        elif Nlastdims == 4:
            loadvar = var[prevdims[0], prevdims[1], :, :, :, :]
        else:
            loadvar = var[prevdims[0], :, :, :, :]
    else:
        print errormsg
        print '  ' + fname + ' variable size ', Ndims, ' not ready!!!'
        print errormsg
        quit(-1)

    return loadvar

def indices_mask(mask):
    """ Function to retrieve the indices from a 2D mask
    mask= 2D mask matrix
    mat=np.array(range(100), dtype=float).reshape((10,10))
    >>> mask1= mat >= 45. 
    >>> mask2= mat <= 55.
    >>> mask = mask1 * mask2
    indices_mask(mask)
    [[4 5]
    [4 6]
    [4 7]
    [4 8]
    [4 9]
    [5 0]
    [5 1]
    [5 2]
    [5 3]
    [5 4]
    [5 5]]
    """
    dy=mask.shape[0]
    dx=mask.shape[1]

    Nptmask = np.sum(mask)
    maskindices = np.zeros((Nptmask,2), dtype=int)
    ij=0
    for j in range(dy):
        for i in range(dx):
            if mask[j,i] == True:
                maskindices[ij,0] = j
                maskindices[ij,1] = i
                ij = ij + 1

    return maskindices

def typemod(value, typeval):
    """ Function to give back a value in a given dtype
    >>> print(typemod(8.2223, 'np.float64'))
    <type 'numpy.float64'>
    >>> print(typemod(8.2223, 'tuple'))
    <type 'tuple'>
    """
    fname='typemod'

    if typeval == 'int':
        return int(value)
    elif typeval == 'long':
        return long(value)
    elif typeval == 'float':
        return float(value)
    elif typeval == 'complex':
        return complex(value)
    elif typeval == 'str':
        return str(value)
    elif typeval == 'bool':
        return bool(value)
    elif typeval == 'list':
        newval = []
        newval.append(value)
        return newval
    elif typeval == 'dic':
        newval = {}
        newval[value] = value
        return newval
    elif typeval == 'tuple':
        newv = []
        newv.append(value)
        newval = tuple(newv)
        return newval
    elif typeval == 'np.int8':
        return np.int8(value)
    elif typeval == 'np.int16':
        return np.int16(value)
    elif typeval == 'np.int32':
        return np.int32(value)
    elif typeval == 'np.int64':
        return np.int64(value)
    elif typeval == 'np.uint8':
        return np.uint8(value)
    elif typeval == 'np.uint16':
        return np.uint16(value)
    elif typeval == 'np.np.uint32':
        return np.uint32(value)
    elif typeval == 'np.uint64':
        return np.uint64(value)
    elif typeval == 'np.float':
        return np.float(value)
    elif typeval == 'np.float16':
        return np.float16(value)
    elif typeval == 'np.float32':
        return np.float32(value)
    elif typeval == 'float32':
        return np.float32(value)
    elif typeval == 'np.float64':
        return np.float64(value)
    elif typeval == 'float64':
        return np.float64(value)
    elif typeval == 'np.complex':
        return np.complex(value)
    elif typeval == 'np.complex64':
        return np.complex64(value)
    elif typeval == 'np.complex128':
        return np.complex128(value)
    else:
        print errormsg
        print fname + ':  data type "' + typeval + '" is not ready !'
        print errormsg
        quit(-1)

    return

def diffdate_units(diffdate, units):
    """ Function to transform a difference of dates to a certain units
    diffdate = difference of dates (as timedelta)
    units = units to transform: 'weeks', 'days', 'hours', 'minutes', 'seconds', 'miliseconds', 'nanoseconds'
    >>> diffdate_units(dt.datetime(1976, 2, 17, 8, 32, 5) - dt.datetime(1949, 12, 1, 0, 0, 0), 'seconds')
    827224325.0
    """
    fname = 'diffdate_units'

    if units == 'weeks':
        diffunits = diffdate.days/7. + diffdate.seconds/(3600.*24.*7.)
    elif units == 'days':
        diffunits = diffdate.days + diffdate.seconds/(3600.*24.)
    elif units == 'hours':
        diffunits = diffdate.days*24. + diffdate.seconds/(3600.)
    elif units == 'minutes':
        diffunits = diffdate.days*24.*60. + diffdate.seconds/(60.)
    elif units == 'seconds':
        diffunits = diffdate.days*24.*3600. + diffdate.seconds
    elif units == 'miliseconds':
        diffunits = diffdate.days*24.*3600.*1000. + diffdate.seconds*1000.
    elif units == 'nanoseconds':
        diffunits = diffdate.days*24.*3600.*1000000. + diffdate.seconds*1000000.
    else:
        print errormsg
        print fname + ': time units "' + units + '" not ready!!!'
        print errormsg
        quit(-1)

    return diffunits

def days_month(year,month):
    """ Function to give the number of days of a month of a given year
    >>> days_month(1976,2)
    29
    """
    import datetime as dt

    date1=dt.date(year,month,1)
    date2=dt.date(year,month+1,1)

    diffdate = date2-date1
    return diffdate.days

def mid_period(yr1,mon1,day1,hour1,min1,sec1,yr2,mon2,day2,hour2,min2,sec2):
    """ Function to give the mid poiint of a period
    >>> mid_period(1976,2,1,0,0,0,1976,3,1,0,0,0)
    [1976    2   15   12    0    0]
    """
    import datetime as dt

    date1=dt.datetime(yr1,mon1,day1,hour1,min1,sec1)
    date2=dt.datetime(yr2,mon2,day2,hour2,min2,sec2)

    diffdate = date2-date1
    diffseconds = diffdate.days*24*3600 + diffdate.seconds
    diff2 = dt.timedelta(seconds=diffseconds/2)
    datenew = date1 + diff2
    datenewvals = np.array([datenew.year, datenew.month, datenew.day, datenew.hour,  \
      datenew.minute, datenew.second])

    return datenewvals

def days_year(year):
    """ Function to give the number of days of a year
    >>> days_year(1976)
    366
    """
    import datetime as dt

    date1=dt.date(year,1,1)
    date2=dt.date(year+1,1,1)

    diffdate = date2-date1
    return diffdate.days

def days_period(yeari,yearf):
    """ Function to give the number of days for a period
    >>> days_period(1976,1980)
    1827
    """

    ndays=0    
    for iyr in range(yearf-yeari+1):
        ndays=ndays+days_year(yeari+iyr)

    return ndays

class ysubsetstats(object):
    """ Class to compute multi-year statistics of a given subset of values 
    values= values to use
    dates=dates of the values in matrix format ([year], [month], [day], [hour], [minute], [second])
    sec=section to use for the period='year', 'month', 'day', 'hour', 'minute', 'second'
    vali=initial value of the period
    valf=final value of the period
    Nq= number of quantiles (20 for 5% bins, it is going to produce 21 values)
    missval= missing value
      self.Nvalues = Number of non masked values
      self.min = minimum non masked value
      self.max = maximum non masked value
      self.mean = mean non masked value
      self.mean2 = quandratic mean non masked value
      self.stdev = standard deviation non masked value
      self.quantiles = quantiles non masked value
    """  
    def __init__(self, values, dates, sec, vali, valf, Nq, missVal):
        import numpy.ma as ma

        fname = 'ysubsetstats'
  
        if values is None:
            self.Nvalues = None
            self.min = None

            self.max = None
            self.mean = None
            self.mean2 = None
            self.stdev = None
            self.quantiles = None

        else:
            Nvalues = len(values)
            Ndates = len(dates)

            if Nvalues != Ndates:
                print errormsg
                print fname + ': number of values ', Nvalues,' does not coincide with the number of dates ',Ndates, '!!!!!!'
                print errormsg
                quit(-1)

            initialsub = np.zeros(Nvalues, dtype=bool)
            endsub = initialsub.copy()

            if missVal > 0.:
                valuesmask = ma.masked_greater_equal(values, missVal*0.9)
            else:
                valuesmask = ma.masked_less_equal(values, missVal*0.9)

            if sec == 'year':
                isec = 0
            elif sec == 'month':
                isec = 1
            elif sec == 'day':
                isec = 2
            elif sec == 'hour':
                isec = 3
            elif sec == 'minute':
                isec = 4
            elif sec == 'second':
                isec = 5
            else:
                print errormsg
                print fname + ': sction "' + sec + '" not ready!!!!'
                print errormsg
                quit(-1)
    
            if vali < valf:
                timesmaskperi = ma.masked_less(dates[:,isec], vali)
                timesmaskperf = ma.masked_greater(dates[:,isec], valf)
            else:
                timesmaskperi = ma.masked_less(dates[:,isec], valf)
                timesmaskperf = ma.masked_greater(dates[:,isec], vali)

            finalmask = valuesmask.mask+timesmaskperi.mask+timesmaskperf.mask
            finalvalues = ma.array(values, mask=finalmask)

            finalvalues2 = finalvalues*finalvalues

            self.Nvalues = finalvalues.count()
            self.min = finalvalues.min()
            self.max = finalvalues.max()
            self.mean = finalvalues.mean()
            self.mean2 = finalvalues2.mean()
            self.stdev = finalvalues.std()
            self.quantiles = mask_quantiles(finalvalues, Nq) 

        return

def give_fix1dim_values(mat,iddim,dimval):
    """ Function to return matrix values when one dimension has been removed at a 
      given value
    mat = matrix
    iddim = dimesion to remove (starting from 0 tacking the right most)
    dimval = value of the dimension

    >>> matrix = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
    give_fix1dim_values(matrix,2,2)
    [[[  0.   1.   2.]
      [  3.   4.   5.]
      [  6.   7.   8.]]

     [[  9.  10.  11.]
      [ 12.  13.  14.]
      [ 15.  16.  17.]]]
    """
    fname='give_fix1dim_values'

# Not working ma.squeeze
##    import numpy.ma as ma
##    print ma.squeeze(mat, axis = (idim,))
    
    matshape = mat.shape
    mattype = mat.dtype
    matmask = np.ones( tuple(matshape), dtype=bool)

    zerodims=np.ones( (6), dtype=int)
    Ndims = len(matshape)
    if Ndims > 6: 
        print errmsg
        print '  ' + fname + ' number of dimensions: ',Ndims,' not ready!!!!!'
        quit(-1)

    zerodims[5-Ndims+1:6] = matshape

    newdims = list(matshape)
    newdims[Ndims-1-iddim] = matshape[iddim]-1

    mask = np.ones(tuple(zerodims), dtype=bool)
    vals = np.zeros(tuple(zerodims), dtype=mattype)
    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],
      0:zerodims[5]] = mat

    for i5 in range(zerodims[5]):
        for i4 in range(zerodims[4]):
            for i3 in range(zerodims[3]):
                for i2 in range(zerodims[2]):
                    for i1 in range(zerodims[1]):
                        for i0 in range(zerodims[0]):
#                            print 'iddim:',iddim,'dimval:',dimval,'-->',i0,i1,i2,i3,i4,i5
                            if iddim == 5 and dimval == i0:
                                mask[i0,i1,i2,i3,i4,i5] = False
                            elif iddim == 4 and dimval == i1:
                                mask[i0,i1,i2,i3,i4,i5] = False
                            elif iddim == 3 and dimval == i2:
                                mask[i0,i1,i2,i3,i4,i5] = False
                            elif iddim == 2 and dimval == i3:
                                mask[i0,i1,i2,i3,i4,i5] = False
                            elif iddim == 1 and dimval == i4:
                                mask[i0,i1,i2,i3,i4,i5] = False
                            elif iddim == 0 and dimval == i5:
                                mask[i0,i1,i2,i3,i4,i5] = False

    newmat = vals[mask].reshape(tuple(newdims))
#    print 'newmat dim:',iddim,'val:',dimval,':',newdims,'______'
#    print newmat
   
    return newmat

def portion_fix1dim_values(mat,iddim,dimval):
    """ Function to return that portion of the matrix values when one dimension has
      been removed at a given value
    mat = matrix
    iddim = dimesion to remove (starting from 0 tacking the right most)
    dimval = value of the dimension

    >>> matrix = np.array( (range(27)), dtype=np.float).reshape(3,3,3)
    portion_fix1dim_values(matrix,2,2)
    [[ 18.  19.  20.]
     [ 21.  22.  23.]
     [ 24.  25.  26.]]
    """
    fname='portion_fix1dim_values'

# Not working ma.squeeze
##    import numpy.ma as ma
##    print ma.squeeze(mat, axis = (idim,))
    
    matshape = mat.shape
    mattype = mat.dtype
    matmask = np.zeros( tuple(matshape), dtype=bool)

    zerodims=np.ones( (6), dtype=int)
    Ndims = len(matshape)

    if Ndims > 6: 
        print errmsg
        print '  ' + fname + ' number of dimensions: ',Ndims,' not ready!!!!!'
        quit(-1)

    zerodims[5-Ndims+1:6] = matshape

    newdims = []
    for idim in range(Ndims):
        if idim != iddim:
            newdims.append(matshape[idim])

    newdims.reverse()

    mask = np.zeros(tuple(zerodims), dtype=bool)
    vals = np.zeros(tuple(zerodims), dtype=mattype)
    vals[0:zerodims[0],0:zerodims[1],0:zerodims[2],0:zerodims[3],0:zerodims[4],
      0:zerodims[5]] = mat

    for i5 in range(zerodims[5]):
        for i4 in range(zerodims[4]):
            for i3 in range(zerodims[3]):
                for i2 in range(zerodims[2]):
                    for i1 in range(zerodims[1]):
                        for i0 in range(zerodims[0]):
#                            print 'iddim:',iddim,'dimval:',dimval,'-->',i0,i1,i2,i3,i4,i5
                            if iddim == 5 and dimval == i0:
                                mask[i0,i1,i2,i3,i4,i5] = True
                            elif iddim == 4 and dimval == i1:
                                mask[i0,i1,i2,i3,i4,i5] = True
                            elif iddim == 3 and dimval == i2:
                                mask[i0,i1,i2,i3,i4,i5] = True
                            elif iddim == 2 and dimval == i3:
                                mask[i0,i1,i2,i3,i4,i5] = True
                            elif iddim == 1 and dimval == i4:
                                mask[i0,i1,i2,i3,i4,i5] = True
                            elif iddim == 0 and dimval == i5:
                                mask[i0,i1,i2,i3,i4,i5] = True

#    print 'newmat dim:',iddim,'val:',dimval,':',newdims,'______'
    newmat = vals[mask].reshape(tuple(newdims))
#    print newmat
   
    return newmat

def retype(val, vtype):
    """ Function to transform a value to a given type
    retype(val, vtype)
      [val]= value
      [vtype]= type to transform
    >>> retype(0, type(True))
    False
    """

    fname = 'retype'

    if val == 'h':
        print fname + '_____________________________________________________________'
        print retype.__doc__
        quit()

    if vtype == type(int(1)):
        newval = int(val)
    elif vtype == type(float(1)):
        newval = float(val)
#    elif vtype == type(float32(1)):
#        newval = float32(val)
    elif vtype == type(np.int(1)):
        typeinf = np.iinfo(np.int)
        if (abs(np.float64(val)) > typeinf.max):
            print warnmsg
            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
               vtype,' limits!!'
            print '  Changing value to kind largest allowed value:', typeinf.max
            newval = abs(np.float64(val))/np.float64(val) * np.int(typeinf.max)
        else:
            newval = np.int(val)
    elif vtype == type(np.int16(1)):
        typeinf = np.iinfo(np.int16)
        if (abs(np.float64(val)) > typeinf.max):
            print warnmsg
            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
               vtype,' limits!!'
            print '  Changing value to kind largest allowed value:', typeinf.max
            newval = abs(np.float64(val))/np.float64(val) * np.int16(typeinf.max)
        else:
            newval = np.int16(val)
    elif vtype == type(np.int32(1)):
        typeinf = np.iinfo(np.int32)
        if (np.abs(np.float64(val)) > typeinf.max):
            print warnmsg
            print '  ' + fname + ': value to transform ', val,' overpasses type:',   \
               vtype,' limits!!'
            print '  Changing value to kind largest allowed value:', typeinf.max
            newval = abs(np.float64(val))/np.float64(val) * np.int32(typeinf.max)
        else:
            newval = np.int32(np.float64(val))
    elif vtype == type(np.int64(1)):
        newval = np.int64(val)
    elif vtype == type(np.float(1)):
        newval = np.float(val)
#    elif vtype == type(np.float16(1)):
#        newval = np.float16(val)
    elif vtype == type(np.float32(1)):
        newval = np.float32(val)
    elif vtype == type(np.float64(1)):
        newval = np.float64(val)
    elif vtype == type(True):
        if val == 0:
            newval = False
        else:
            newval = True
    elif vtype == '|S1':
        newval = str(val)
    else:
        print errormsg
        print '  ' + fname + ': variable type "', vtype, '" not ready!!!'
        quit(-1)

    return newval

class stats_space2D(object):
    """spatial statistics for a 2D file:
      vals= variable ro check (assuming var[t,dy,dx])
        self.minv[t]: spatial minimum at each time-step
        self.maxv[t]: spatial maximum at each time-step
        self.meanv[t]: spatial mean at each time-step
        self.mean2v[t]: spatial quadratic mean at each time-step
        self.stdv[t]: spatial standard deviation at each time-step
        self.anomv[t]: spatial anomaly from the whole mean at each time-step
    """
  
    def __init__(self, vals):
        from scipy import stats as sts
        fname = 'stats_space2D'

        if vals1 == 'h':
            print fname + '_____________________________________________________________'
            print stats_space2Dvars.__doc__
            quit()

        if vals1 is None:
            self.minv = None
            self.maxv = None
            self.meanv = None
            self.mean2v = None
            self.stdv = None
        else:
            dimt = vals.shape[0]
            stats=np.zeros((dimt,5), dtype=np.float)
            absmean = np.mean(vals)

            for it in range(dimt):
                stats[it,0]=np.min(vals[it,:,:])
                stats[it,1]=np.max(vals[it,:,:])
                stats[it,2]=np.mean(vals[it,:,:])
                stats[it,3]=np.mean(vals[it,:,:]*vals[it,:,:])
                stats[it,4]=absmean - stats[it,2]

            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
            stats = np.where(stats is nan, fillValue, stats)
            stats = np.where(stats is np.inf, fillValue, stats)
            stats = np.where(stats is None, fillValue, stats)

            self.minv=stats[:,0]
            self.maxv=stats[:,1]
            self.meanv=stats[:,2]
            self.mean2v=stats[:,3]
            self.stdv=np.sqrt(stats[:,3]-stats[:,2]*stats[:,2])
            self.anomv=stats[:,4]

        return

class stats_space2Dvars(object):
    """spatial statistics beween 2 2D files:
    valsA = variable 1 (assuming var[t,dy,dx])
    valsB = variable 2 (assuming var[t,dy,dx])
    self.min[t,var], self.max[t,var], self.mean[t,var], self.mean2[t,var], 
      self.std[t,var]
      [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]

    self.mae=mean[t,abs(var1-var2)]
    self.correlation=correlation[var1,var2] (and p-value) 
    self.bias=[t,mean(var1)-mean(var2)]
    self.meancorrelation=correlation[mean(var1),mean(var2)] (and p-value) 
    """

    def __init__(self, vals1, vals2):
        from scipy import stats as sts
        fname = 'stats_space2Dvars'

        if vals1 == 'h':
            print fname + '_____________________________________________________________'
            print stats_space2Dvars.__doc__
            quit()

        if vals1 is None:
            self.minv1Av2 = None
            self.maxv1Av2 = None
            self.meanv1Av2 = None
            self.mean2v1Av2 = None
            self.stdv1Av2 = None
            self.minv1Sv2 = None
            self.maxv1Sv2 = None
            self.meanv1Sv2 = None
            self.mean2v1Sv2 = None
            self.stdv1Sv2 = None
            self.minv1Dv2 = None
            self.maxv1Dv2 = None
            self.meanv1Dv2 = None
            self.mean2v1Dv2 = None
            self.stdv1Dv2 = None
            self.minv1Pv2 = None
            self.maxv1Pv2 = None
            self.meanv1Pv2 = None
            self.mean2v1Pv2 = None
            self.stdv1Pv2 = None
            self.mae = None
            self.corr = None
        else:
            dimt = vals1.shape[0]
            stats=np.zeros((dimt,26), dtype=np.float)
            meanvals1 = np.zeros((dimt), dtype=np.float)
            meanvals2 = np.zeros((dimt), dtype=np.float)
    
            for it in range(dimt):
                v1 = vals1[it,:,:].flatten() 
                v2 = vals2[it,:,:].flatten() 

# add
                vs = v1 + v2
                stats[it,0] = np.min(vs)
                stats[it,1] = np.max(vs)
                stats[it,2] = np.mean(vs)
                stats[it,3] = np.mean(vs*vs)
                stats[it,4] = np.sqrt(stats[it,3] - stats[it,2]*stats[it,2])
# sub
                stats[it,20] = np.mean(np.abs(v1-v2))
                vs = v1 - v2
                stats[it,5] = np.min(vs)
                stats[it,6] = np.max(vs)
                stats[it,7] = np.mean(vs)
                stats[it,8] = np.mean(vs*vs)
                stats[it,9] = np.sqrt(stats[it,8] - stats[it,7]*stats[it,7])

# mul
                vs = v1 * v2
                stats[it,10] = np.min(vs)
                stats[it,11] = np.max(vs)
                stats[it,12] = np.mean(vs)
                stats[it,13] = np.mean(vs*vs)
                stats[it,14] = np.sqrt(stats[it,13] - stats[it,12]*stats[it,12])
# div
                vs = v1 / v2
                stats[it,15] = np.min(vs)
                stats[it,16] = np.max(vs)
                stats[it,17] = np.mean(vs)
                stats[it,18] = np.mean(vs*vs)
                stats[it,19] = np.sqrt(stats[it,18] - stats[it,17]*stats[it,17])
# corr
                stats[it,21], stats[it,22] = mask_pearsonr(v1, v2)

# Mean values
                meanvals1[it] = np.mean(v1)
                meanvals2[it] = np.mean(v2)
                stats[it,23] =  meanvals1[it] - meanvals2[it]
            
            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
            stats = np.where(stats is np.nan, fillValue, stats)
            stats = np.where(stats is np.inf, fillValue, stats)
            stats = np.where(stats is None, fillValue, stats)

            self.minv1Av2 = stats[:,0]
            self.maxv1Av2 = stats[:,1]
            self.meanv1Av2 = stats[:,2]
            self.meanv21Av2 = stats[:,3]
            self.stdv1Av2 = stats[:,4]

            self.minv1Sv2 = stats[:,5]
            self.maxv1Sv2 = stats[:,6]
            self.meanv1Sv2 = stats[:,7]
            self.meanv21Sv2 = stats[:,8]
            self.stdv1Sv2 = stats[:,9]

            self.minv1Pv2 = stats[:,10]
            self.maxv1Pv2 = stats[:,11]
            self.meanv1Pv2 = stats[:,12]
            self.meanv21Pv2 = stats[:,13]
            self.stdv1Pv2 = stats[:,14]

            self.minv1Dv2 = stats[:,15]
            self.maxv1Dv2 = stats[:,16]
            self.meanv1Dv2 = stats[:,17]
            self.meanv21Dv2 = stats[:,18]
            self.stdv1Dv2 = stats[:,19]

            self.mae = stats[:,20]
            self.corr = stats[:,21]
            self.p_value = stats[:,22]

            self.bias = stats[:,23]
            self.meancorr, self.meanp_value = sts.pearsonr(meanvals1, meanvals2)

        return

class stats_time2D(object):
    """temporal statistics for a 2D file:
      vals= variable ro check (assuming var[t,dy,dx])
        self.minv[t]: temporal minimum at each grid point
        self.maxv[t]: temporal maximum at each grid point
        self.meanv[t]: temporal mean at each grid point
        self.mean2v[t]: temporal quadratic mean at each grid point
        self.stdv[t]: temporal standard deviation at each grid point
        self.anomv[t]: temporal anomaly from the whole mean at each grid point
    """
  
    def __init__(self, vals):
        from scipy import stats as sts
        fname = 'stats_time2D'

        if vals1 == 'h':
            print fname + '_____________________________________________________________'
            print stats_time2Dvars.__doc__
            quit()

        if vals1 is None:
            self.minv = None
            self.maxv = None
            self.meanv = None
            self.mean2v = None
            self.stdv = None
        else:
            dimx = vals.shape[2]
            dimy = vals.shape[1]
            stats=np.zeros((dimy,dimx,5), dtype=np.float)
            absmean = np.mean(vals,axis=0)

            stats[:,:,0]=np.min(vals, axis=0)
            stats[:,:,1]=np.max(vals, axis=0)
            stats[:,:,2]=np.mean(vals, axis=0)
            stats[:,:,3]=np.mean(vals*vals, axis=0)
            stats[:,:,4]=absmean - stats[:,:,2]

            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
            stats = np.where(stats is np.nan, fillValue, stats)
            stats = np.where(stats is np.inf, fillValue, stats)
            stats = np.where(stats is None, fillValue, stats)

            self.minv=stats[:,:,0]
            self.maxv=stats[:,:,1]
            self.meanv=stats[:,:,2]
            self.mean2v=stats[:,:,3]
            self.stdv=np.sqrt(stats[:,:,3]-stats[:,:,2]*stats[:,:,2])
            self.anomv=stats[:,:,4]

        return

class stats_time2Dvars(object):
    """temporal statistics beween 2 2D files:
    valsA = variable 1 (assuming var[t,dy,dx])
    valsB = variable 2 (assuming var[t,dy,dx])
    self.min[dy,dx,var], self.max[dy,dx,var], self.mean[dy,dx,var], 
      self.mean2[dy,dx,var], self.std[dy,dx,var]
      [var] = var1+var2[v1Av2], var1-var2[v1Sv2], var1/var2[v1Dv2], var1*var2[v1Pv2]

    self.mae=mean[dy,dx,abs(var1-var2)]
    self.correlation=correlation[var1,var2] (and p-value) 
    self.bias=[dy,dx,mean(var1)-mean(var2)]
    self.meancorrelation=correlation[mean(var1),mean(var2)] (and p-value) 
    """

    def __init__(self, vals1, vals2):
        from scipy import stats as sts
        fname = 'stats_time2Dvars'

        if vals1 == 'h':
            print fname + '_____________________________________________________________'
            print stats_time2Dvars.__doc__
            quit()

        if vals1 is None:
            self.minv1Av2 = None
            self.maxv1Av2 = None
            self.meanv1Av2 = None
            self.mean2v1Av2 = None
            self.stdv1Av2 = None
            self.minv1Sv2 = None
            self.maxv1Sv2 = None
            self.meanv1Sv2 = None
            self.mean2v1Sv2 = None
            self.stdv1Sv2 = None
            self.minv1Dv2 = None
            self.maxv1Dv2 = None
            self.meanv1Dv2 = None
            self.mean2v1Dv2 = None
            self.stdv1Dv2 = None
            self.minv1Pv2 = None
            self.maxv1Pv2 = None
            self.meanv1Pv2 = None
            self.mean2v1Pv2 = None
            self.stdv1Pv2 = None
            self.mae = None
            self.corr = None
        else:
            dimx = vals1.shape[1]
            dimy = vals1.shape[1]
            stats=np.zeros((dimy,dimx,24), dtype=np.float)
            meanvals1 = np.zeros((dimy,dimx), dtype=np.float)
            meanvals2 = np.zeros((dimy,dimx), dtype=np.float)
# add
            vs = vals1 + vals2
            stats[:,:,0] = np.min(vs,axis=0)
            stats[:,:,1] = np.max(vs,axis=0)
            stats[:,:,2] = np.mean(vs,axis=0)
            stats[:,:,3] = np.mean(vs*vs,axis=0)
            stats[:,:,4] = np.sqrt(stats[:,:,3] - stats[:,:,2]*stats[:,:,2])
# sub
            stats[:,:,20] = np.mean(np.abs(vals1-vals2), axis=0)
            vs = vals1 - vals2
            stats[:,:,5] = np.min(vs,axis=0)
            stats[:,:,6] = np.max(vs,axis=0)
            stats[:,:,7] = np.mean(vs,axis=0)
            stats[:,:,8] = np.mean(vs*vs,axis=0)
            stats[:,:,9] = np.sqrt(stats[:,:,8] - stats[:,:,7]*stats[:,:,7])
# mul
            vs = vals1 * vals2
            stats[:,:,10] = np.min(vs,axis=0)
            stats[:,:,11] = np.max(vs,axis=0)
            stats[:,:,12] = np.mean(vs,axis=0)
            stats[:,:,13] = np.mean(vs*vs,axis=0)
            stats[:,:,14] = np.sqrt(stats[:,:,13] - stats[:,:,12]*stats[:,:,12])
# div
            vs = vals1 / vals2
            stats[:,:,15] = np.min(vs,axis=0)
            stats[:,:,16] = np.max(vs,axis=0)
            stats[:,:,17] = np.mean(vs,axis=0)
            stats[:,:,18] = np.mean(vs*vs,axis=0)
            stats[:,:,19] = np.sqrt(stats[:,:,18] - stats[:,:,17]*stats[:,:,17])


# Mean values
            meanvals1[:,:] = np.mean(vals1,axis=0)
            meanvals2[:,:] = np.mean(vals2,axis=0)
            stats[:,:,23] =  meanvals1[:,:] - meanvals2[:,:]
            
# corr
            self.meancorr, self.meanp_value = sts.pearsonr(meanvals1.flatten(), meanvals2.flatten())

            for j in range(dimy):
                for i in range(dimx):
                    stats[j,i,21], stats[j,i,22] = sts.pearsonr(vals1[:,j,i], vals2[:,j,i])

            stats = np.where(stats > 0.1*fillValue, fillValue, stats)
            stats = np.where(stats is np.nan, fillValue, stats)
            stats = np.where(stats is np.inf, fillValue, stats)
            stats = np.where(stats is None, fillValue, stats)

            self.minv1Av2 = stats[:,:,0]
            self.maxv1Av2 = stats[:,:,1]
            self.meanv1Av2 = stats[:,:,2]
            self.meanv21Av2 = stats[:,:,3]
            self.stdv1Av2 = stats[:,:,4]

            self.minv1Sv2 = stats[:,:,5]
            self.maxv1Sv2 = stats[:,:,6]
            self.meanv1Sv2 = stats[:,:,7]
            self.meanv21Sv2 = stats[:,:,8]
            self.stdv1Sv2 = stats[:,:,9]

            self.minv1Pv2 = stats[:,:,10]
            self.maxv1Pv2 = stats[:,:,11]
            self.meanv1Pv2 = stats[:,:,12]
            self.meanv21Pv2 = stats[:,:,13]
            self.stdv1Pv2 = stats[:,:,14]

            self.minv1Dv2 = stats[:,:,15]
            self.maxv1Dv2 = stats[:,:,16]
            self.meanv1Dv2 = stats[:,:,17]
            self.meanv21Dv2 = stats[:,:,18]
            self.stdv1Dv2 = stats[:,:,19]

            self.mae = stats[:,:,20]
            self.corr = stats[:,:,21]
            self.p_value = stats[:,:,22]

            self.bias = stats[:,:,23]

        return

#vals1 = np.arange(27).reshape(3,3,3)*1.
#vals2 = np.arange(1,28).reshape(3,3,3)*1.

#printing_class(stats_time2Dvars(vals1,vals2))

def file_nlines(filen,char):
    """ Function to provide the number of lines of a file
    filen= name of the file
    char= character as no line
    >>> 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: 
        if line[0:1] != char: nlines = nlines + 1

    fo.close()

    return nlines

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 table_tex(tablevals, colnames, rownames, of):
    """ Function to write into a LaTeX tabukar from a table of values
      tablevals = (ncol nrow) of values
      colnames = list with ncol labels for the columns
      rownames = list with nrow labels for the rows
      of= object ASCII file to write the table
    """
    errormsg = 'ERROR -- error -- ERROR -- error'

    fname = 'table_tex'

    Ncol = tablevals.shape[0]
    Nrow = tablevals.shape[1]

    if Ncol + 1 != len(colnames):
        print errormsg
        print '  ' + fname + ': wrong number of column names!!'
        print '    data has:', Ncol, 'and you provided:', len(colnames)
        print '    remember to provide one more for the label of the rows!'
        print '    you provide:',colnames
        quit(-1)

    if Nrow != len(rownames):
        print errormsg
        print '  ' + fname + ': wrong number of row names!!'
        print '    data has:', Nrow, 'and you provided:', len(rownames)
        print '    you provide:',rownames
        quit(-1)

    colcs = ''
    colns = ''
    for icol in colnames:
        colcs = colcs + 'c'
        if icol == colnames[0]:
            colns = ' {\\bfseries{' + icol + '}}'
        else:
            colns = colns + ' & {\\bfseries{' + icol + '}}'

    of.write('\n')
    of.write("%Automatically written file from function '" + fname + "'\n")
    of.write('\\begin{tabular}{l'+colcs+'}\n')
    of.write(colns + ' \\\\ \\hline\n')

    ir = 0
    for irow in rownames:
        rowns = '{\\bfseries{' + irow + '}}'
        for ic in range(Ncol):
            rowns = rowns + ' & ' + str(tablevals[ic,ir])
        rowns = rowns + ' \\\\'

        of.write(rowns + '\n')
        ir = ir + 1

    of.write('\\end{tabular}\n')

    return

def radius_dist(dx,dy,ptx,pty):
    """ Function to generate a matrix with the distance at a given point
    radius_dist(dx,dy,ptx,pty)
      [dx/y]: dimension of the matrix
      [ptx/y]: grid point coordinates of the point
    >>> radius_dist(3,5,2,2)
    [[ 1.41421356  1.          1.41421356  2.23606798  3.16227766]
     [ 1.          0.          1.          2.          3.        ]
     [ 1.41421356  1.          1.41421356  2.23606798  3.16227766]]
    """

    fname = 'radius_dist'

    if ptx < 0 or ptx > dx-1 or pty < 0 or pty > dy-1:
        print errormsg
        print '  ' + fname + ': wrong point coordinates:',dx,',',dy,'for matrix;',dx \
         ,'x',dy
        quit(-1)

    xdist =  np.zeros((dx,dy), dtype=np.float)
    ydist =  np.zeros((dx,dy), dtype=np.float)
    dist =  np.zeros((dx,dy), dtype=np.float)

    for ix in range(dx):
        xdist[ix,:] = np.float(ix-ptx)
    for iy in range(dy):
        ydist[:,iy] = np.float(iy-pty)

    dist = np.sqrt(xdist*xdist + ydist*ydist)

    return dist

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 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 vertical_interpolation2D(varb,vart,zorigvals,znewval,kinterp):
    """ Function to vertically integrate a 3D variable
    vertical_interpolation2D(varb,vart,zorigvals,znewval)
      varb= values at the base of the interval of interpolation (2D field)
      vart= values at the top of the interval of interpolation (2D field)
      zorigvals= pair of original values (2, 2D field)
      znewval= new value to interpolate
      Possible cases:
        zorigvals[0,:,:] <= znewval < zorigvals[1,:,:]
        znewval <= zorigvals[0,:,:] < zorigvals[1,:,:]
        zorigvals[0,:,:] < zorigvals[1,:,:] <= znewval
      kinterp: kind of interpolation
        'lin': linear
    >>> dx=5
    >>> dy=7
    >>> vals1 = np.ones((dy,dx), dtype=np.float)
    >>> vals2 = np.ones((dy,dx), dtype=np.float)*2.

    >>> zbase = np.zeros((2,dy,dx), dtype=np.float)
    >>> zbase[0,:,:] = 0.5
    >>> zbase[1,:,:] = 1.

    >>> vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')
    [[ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]
     [ 1.5  1.5  1.5  1.5  1.5]]
    """

    fname = 'vertical_interpolation2D'

    if varb == 'h':
        print fname + '_____________________________________________________________'
        print vertical_interpolation2D.__doc__
        quit()

    newvar = np.zeros((varb.shape), dtype=np.float)
    if kinterp == 'lin':
##        print '  ' + fname + ' Vertical linear interpolation at',znewval
# Standard linear interpolation (y = a + b*incz)
#   a = zorig[0], b = (vart-varb)/(zorig[1]-zorig[0])), incz = znewval - zorig[0]
        a = varb
        b = np.where(zorigvals[1,:,:] == zorigvals[0,:,:], 0.,                       \
          (vart - varb)/(zorigvals[1,:,:] - zorigvals[0,:,:]))
        incz = np.ones((varb.shape), dtype=np.float)*znewval - zorigvals[0,:,:]

        newvar = a + b*incz
# Too code for not be used... but maybe?
#        dx=varb.shape[1]
#        dy=varb.shape[0]
#        for j in range(dy):
#            for i in range(dx):
#                  if zorigvals[1,j,i] == zorigvals[0,j,i]: print 'equals!',          \
#                    zorigvals[1,j,i], zorigvals[0,j,i],':',newvar[j,i],'@',vart[j,i]
#                  if zorigvals[0,j,i] != zorigvals[0,j,i]: print '0 Nan',            \
#                    zorigvals[0,j,i],':',newvar[j,i],'@',vart[j,i]
#                  if zorigvals[1,j,i] != zorigvals[1,j,i]: print '1 Nan',            \
#                    zorigvals[1,j,i],':',newvar[j,i],'@',vart[j,i]
#                  if zorigvals[0,j,i] is None: print '0 None', zorigvals[0,j,i],':', \
#                    newvar[j,i],'@',vart[j,i]
#                  if zorigvals[1,j,i] is None: print '1 None', zorigvals[1,j,i],':', \
#                    newvar[j,i],'@',vart[j,i]
    else:
        print errormsg
        print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
        quit(-1)

    return newvar

#dx=5
#dy=7
#vals1 = np.ones((dy,dx), dtype=np.float)
#vals2 = np.ones((dy,dx), dtype=np.float)*2.

#zbase = np.zeros((2,dy,dx), dtype=np.float)
#zbase[0,:,:] = 0.5
#for i in range(dx):
#    for j in range(dy):
#        zbase[1,j,i] = np.sqrt((j-dy/2.)**2. + (i-dx/2.)**2.) /                      \
#          np.sqrt((dy/2.)**2.+(dy/2.)**2.) + 1.

#zbase[1,:,:] = 1.

#newz = 0.75

#print vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')

def vertical_interpolation(varb,vart,zorigvals,znewval,kinterp):
    """ Function to vertically integrate a 1D variable
    vertical_interpolation(varb,vart,zorigvals,znewval)
      varb= values at the base of the interval of interpolation
      vart= values at the top of the interval of interpolation
      zorigvals= pair of original values (2)
      znewval= new value to interpolate
      Possible cases:
        zorigvals[0,:] <= znewval < zorigvals[1,:]
        znewval <= zorigvals[0,:] < zorigvals[1,:]
        zorigvals[0,:] < zorigvals[1,:] <= znewval
      kinterp: kind of interpolation
        'lin': linear
    >>> dx=5
    >>> vals1 = np.ones((dx), dtype=np.float)
    >>> vals2 = np.ones((dx), dtype=np.float)*2.

    >>> zbase = np.zeros((2,dx), dtype=np.float)
    >>> zbase[0,:] = 0.5
    >>> for i in range(dx):
    >>>     zbase[1,i] = i + 1.

    >>> newz = 0.75
    >>> vertical_interpolation2D(vals1,vals2, zbase, newz,'lin')
    [ 1.5         1.16666667  1.1         1.07142857  1.05555556]
    """

    fname = 'vertical_interpolation'

    if varb == 'h':
        print fname + '_____________________________________________________________'
        print vertical_interpolation.__doc__
        quit()

    newvar = np.zeros((varb.shape), dtype=np.float)
    if kinterp == 'lin':
##        print '  ' + fname + ' Vertical linear interpolation at',znewval
# Standard linear interpolation (y = a + b*incz)
#   a = zorig[0], b = (vart-varb)/(zorig[1]-zorig[0])), incz = znewval - zorig[0]
        a = varb
        b = np.where(zorigvals[1,:] == zorigvals[0,:], 0.,                           \
          (vart - varb)/(zorigvals[1,:] - zorigvals[0,:]))
        incz = np.ones((varb.shape), dtype=np.float)*znewval - zorigvals[0,:]

        newvar = a + b*incz
    else:
        print errormsg
        print '  ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
        quit(-1)

    return newvar

def interpolate_3D(zorigvs, varv, zintvs, kint):
    """ Function to interpolate a 3D variable
    interpolate_3D(zintvs, varv, zorigvals)
      zorigvs= original 3D z values
      varv= 3D variable to interpolate
      zintvs= new series of z values to interpolate
      kint= kind of interpolation
        'lin': linear
    """
    fname = 'interpolate_3D'

    if zorigvs == 'h':
        print fname + '_____________________________________________________________'
        print interpolate_3D.__doc__
        quit()

    Ninv = len(zintvs)
    dimv = varv.shape

# Sense of the original z-variable
    Lzorig = zorigvs.shape[0]
    incz = zorigvs[Lzorig-1,0,0] - zorigvs[0,0,0]

    varin = np.zeros((Ninv,dimv[1],dimv[2]), dtype=np.float)

    for ilev in range(Ninv):
#        print ' vertical interpolate level: ',zintvs[ilev]
        valinf = np.zeros((dimv[1], dimv[2]), dtype=np.float)
        valsup = np.zeros((dimv[1], dimv[2]), dtype=np.float)
        zorigv = np.zeros((2,dimv[1], dimv[2]), dtype=np.float)
        for j in range(valinf.shape[0]):
            for i in range(valinf.shape[1]):
                
                if np.min(zorigvs[:,j,i]) > zintvs[ilev] and incz > 0.:
                    valinf[j,i] = varv[0,j,i]
                    valsup[j,i] = varv[1,j,i]
                    zorigv[0,j,i] = zorigvs[0,j,i]
                    zorigv[1,j,i] = zorigvs[1,j,i]
                elif np.max(zorigvs[:,j,i]) < zintvs[ilev] and incz < 0.:
                    valinf[j,i] = varv[0,j,i]
                    valsup[j,i] = varv[1,j,i]
                    zorigv[0,j,i] = zorigvs[0,j,i]
                    zorigv[1,j,i] = zorigvs[1,j,i]
                elif np.max(zorigvs[:,j,i]) < zintvs[ilev] and incz > 0.:
                    valinf[j,i] = varv[Lzorig-2,j,i]
                    valsup[j,i] = varv[Lzorig-1,j,i]
                    zorigv[0,j,i] = zorigvs[Lzorig-2,j,i]
                    zorigv[1,j,i] = zorigvs[Lzorig-1,j,i]
                elif np.min(zorigvs[:,j,i]) > zintvs[ilev] and incz < 0.:
                    valinf[j,i] = varv[Lzorig-2,j,i]
                    valsup[j,i] = varv[Lzorig-1,j,i]
                    zorigv[0,j,i] = zorigvs[Lzorig-2,j,i]
                    zorigv[1,j,i] = zorigvs[Lzorig-1,j,i]
#                    print 'top: ',i,j,':',zorigv[0,j,i], zintvs[ilev], zorigv[1,j,i]
                else:
                    for z in range(Lzorig-1):
                        if (zorigvs[z,j,i]-zintvs[ilev])*(zorigvs[z+1,j,i]-          \
                          zintvs[ilev]) <= 0.:
                            valinf[j,i] = varv[z,j,i]
                            valsup[j,i] = varv[z+1,j,i]
                            zorigv[0,j,i] = zorigvs[z,j,i]
                            zorigv[1,j,i] = zorigvs[z+1,j,i]
                            break

#        print 'zorigvs: ',zorigvs[:,0,0]
#        print '  valinf:', valinf[0,0]
#        print '  valsup:', valsup[0,0]
#        print '  zorigv 0:',zorigv[0,0,0]
#        print '  zorigv 1:',zorigv[1,0,0]

        varin[ilev,:,:] = vertical_interpolation2D(valinf,valsup,zorigv,             \
          zintvs[ilev], kint)
#        print '  varin:',varin[ilev,0,0]
#        quit()

#    quit()
    return varin

def interpolate_2D(zorigvs0, varv0, zdim, zintvs, kint):
    """ Function to interpolate a 2D variable
    interpolate_2D(zintvs, varv, zorigvals)
      zorigvs= original 2D z values
      varv= 2D variable to interpolate
      zdim= number of the dimension with the z-values
      zintvs= new series of z values to interpolate
      kint= kind of interpolation
        'lin': linear
    """
    fname = 'interpolate_2D'

    if zorigvs0 == 'h':
        print fname + '_____________________________________________________________'
        print interpolate_2D.__doc__
        quit()

    Ninv = len(zintvs)
    if zdim == 0:
        varv = varv0
        zorigvs = zorigvs0
    else:
        varv = np.transpose(varv0)
        zorigvs = np.transpose(zorigvs0)

    dimv = varv.shape

# Sense of the original z-variable
    Lzorig = zorigvs.shape[0]
    incz = zorigvs[Lzorig-1,0] - zorigvs[0,0]

    varin = np.zeros((Ninv,dimv[1]), dtype=np.float)

    for ilev in range(Ninv):
        valinf = np.zeros((dimv[1]), dtype=np.float)
        valsup = np.zeros((dimv[1]), dtype=np.float)
        zorigv = np.zeros((2,dimv[1]), dtype=np.float)
        for i in range(valinf.shape[0]):   
            if np.min(zorigvs[:,i]) > zintvs[ilev] and incz > 0.:
                valinf[i] = varv[0,i]
                valsup[i] = varv[1,i]
                zorigv[0,i] = zorigvs[0,i]
                zorigv[1,i] = zorigvs[1,i]
            elif np.max(zorigvs[:,i]) < zintvs[ilev] and incz < 0.:
                valinf[i] = varv[0,i]
                valsup[i] = varv[1,i]
                zorigv[0,i] = zorigvs[0,i]
                zorigv[1,i] = zorigvs[1,i]
            elif np.max(zorigvs[:,i]) < zintvs[ilev] and incz > 0.:
                valinf[i] = varv[Lzorig-2,i]
                valsup[i] = varv[Lzorig-1,i]
                zorigv[0,i] = zorigvs[Lzorig-2,i]
                zorigv[1,i] = zorigvs[Lzorig-1,i]
            elif np.min(zorigvs[:,i]) > zintvs[ilev] and incz < 0.:
                valinf[i] = varv[Lzorig-2,i]
                valsup[i] = varv[Lzorig-1,i]
                zorigv[0,i] = zorigvs[Lzorig-2,i]
                zorigv[1,i] = zorigvs[Lzorig-1,i]
            else:
                for z in range(Lzorig-1):
                    if (zorigvs[z,i]-zintvs[ilev])*(zorigvs[z+1,i]-                  \
                      zintvs[ilev]) <= 0.:
                        valinf[i] = varv[z,i]
                        valsup[i] = varv[z+1,i]
                        zorigv[0,i] = zorigvs[z,i]
                        zorigv[1,i] = zorigvs[z+1,i]
                        break

        varin[ilev,:] = vertical_interpolation(valinf,valsup,zorigv,                 \
          zintvs[ilev], kint)

    if zdim == 0:
        varin0 = varin
    else:
        varin0 = np.transpose(varin)

    return varin0

def list_toVec(listv, kind):
    """ Function to transform from a list of values to a vector
    list_toVec(listv, kind)
      listv= list with the values 
      kind= kind of values
      >>> list_toVec(['1', 2, '3', 4, 5, '6.9', 7, 9], 'npfloat')
      [ 1.   2.   3.   4.   5.   6.9  7.   9. ]
    """

    fname = 'list_toVec'

    if listv == 'h':
        print fname + '_____________________________________________________________'
        print list_toVec.__doc__
        quit()

    Nvals = len(listv)

    if kind == 'npfloat':
        vecv = np.zeros((Nvals), dtype=np.float)
        for iv in range(Nvals):
            vecv[iv] = np.float(listv[iv])
    else:
        print errormsg
        print '  ' + fname + ': type of value "' + kind + '" not ready!!'
        quit(-1)

    return vecv

def running_mean(values, run):
    """ Function to compute a running mean of a series of values
      running_mean(vals, int)
      [vals]: series of values
      [run]: interval to use to compute the running mean 
        NOTE: resultant size would be the same but with None except at [run/2,size(vals)-run/2]
      >>> running_mean(np.arange(100),10)
      [None None None None None 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5
       15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5 28.5 29.5
       30.5 31.5 32.5 33.5 34.5 35.5 36.5 37.5 38.5 39.5 40.5 41.5 42.5 43.5 44.5
       45.5 46.5 47.5 48.5 49.5 50.5 51.5 52.5 53.5 54.5 55.5 56.5 57.5 58.5 59.5
       60.5 61.5 62.5 63.5 64.5 65.5 66.5 67.5 68.5 69.5 70.5 71.5 72.5 73.5 74.5
       75.5 76.5 77.5 78.5 79.5 80.5 81.5 82.5 83.5 84.5 85.5 86.5 87.5 88.5 89.5
       90.5 91.5 92.5 93.5 None None None None None]
    """
    fname = 'running_mean'

    if values == 'h':
        print fname + '_____________________________________________________________'
        print running_mean.__doc__
        quit()

    runmean = np.ones((len(values)), dtype=np.float)*fillValue

    for ip in range(run/2,len(values)-run/2):
        runmean[ip] = np.mean(values[ip-run/2:ip+run/2])

    runmean = np.where(runmean == fillValue, None, runmean)

    return runmean

def files_folder(folder,headfile):
    """ Function to retrieve a list of files from a folder [fold] and head [headfile]
    >>> files_folder('/media/data2/etudes/WRF_LMDZ/WL_HyMeX_HighRes/medic950116/', 
      'wrfout_d03')
    ['wrfout_d03_1995-01-13_02:00:00', 'wrfout_d03_1995-01-13_04:00:00', 'wrfout_d03_1995-01-13_06:00:00', 
     'wrfout_d03_1995-01-13_06:15:00', 'wrfout_d03_1995-01-13_08:15:00', 'wrfout_d03_1995-01-13_10:15:00',
     'wrfout_d03_1995-01-13_12:15:00', 'wrfout_d03_1995-01-13_14:15:00', 'wrfout_d03_1995-01-13_16:15:00',     
     (...)
     'wrfout_d03_1995-01-17_18:15:00', 'wrfout_d03_1995-01-17_20:15:00', 'wrfout_d03_1995-01-17_22:15:00']
    """
    import subprocess as sub
    fname = 'files_folder'

    if folder == 'h':
        print fname + '_____________________________________________________________'
        print files_folder.__doc__
        quit()

#    ins = folder + "/" + headfile + "*"
    ins = folder + "/"
    print 'ins:',ins
    files = sub.Popen(["/bin/ls","-1",ins], stdout=sub.PIPE)
    fileslist = files.communicate()
    listfiles0 = str(fileslist).split('\\n')

# Filtering output
    Lheader=len(headfile)
    listfiles = []
    for ifile in listfiles0:
        if ifile[0:Lheader] == headfile:
            listfiles.append(ifile)

    return listfiles

def count_cond(var, val, con):
    """ Function to count values of a variable which attain a condition
      [var]= variable
      [val]= value
      [con]= statistics/condition
        'eq': equal than [val]
        'neq': not equal than [val]
        'lt': less than [val]
        'le': less and equal than [val]
        'gt': greater than [val]
        'ge': greater and equal than [val]
        'in': inside the range [val[0]] <= [var] <= [val[1]] (where [val]=[val1, val2]
        'out': inside the range [val[0]] > [var]; [val[1]] < [var] (where [val]=[val1, val2]
      >>> vals = np.arange(100)
      >>> count_cond(vals,50,'eq')
      1
      >>> count_cond(vals,50,'neq')
      99
      >>> count_cond(vals,50,'lt')
      50
      >>> count_cond(vals,50,'le')
      51
      >>> count_cond(vals,50,'gt')
      49
      >>> count_cond(vals,50,'ge')
      50
      >>> count_cond(vals,[25,75],'in')
      51
      >>> count_cond(vals,[25,75],'out')
      49
    """
    fname = 'count_cond'

    cons = ['eq', 'neq', 'lt', 'le', 'gt', 'ge', 'in', 'out']

    if con == 'eq':
        Nvals = np.sum(var == val)
    elif con == 'neq':
        Nvals = np.sum(var != val)
    elif con == 'lt':
        Nvals = np.sum(var < val)
    elif con == 'le':
        Nvals = np.sum(var <= val)
    elif con == 'gt':
        Nvals = np.sum(var > val)
    elif con == 'ge':
        Nvals = np.sum(var >= val)
    elif con == 'in':
        if len(val) != 2:
            print errormsg
            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
              " values!!" 
            print '  we have:',  val
            quit(-1)
        Nvals = np.sum((var >= val[0])*(var <= val[1]))
    elif con == 'out':
        if len(val) != 2:
            print errormsg
            print '  ' + fname + ": for condition '" + con + "' we must have two " + \
              " values!!" 
            print '  we have:',  val
            quit(-1)
        Nvals = np.sum((var < val[0])+(var > val[1]))
    else:
        print errormsg
        print '  ' + fname + ": condition '" + con + "' not ready!!"
        print '  only ready:', cons
        quit(-1)

    return Nvals

def get_namelist_vars(values, namelist):
    """ Function to get namelist-like  values ([varname] = [value]) 
    get_namelist_vars(namelist)
      values= [sectionname],[kout]
        [sectionname]: name of the section from which one want the values ('all', for all)
        [kout]: kind of output
          'tex3': printed output as LaTeX table of three columns \verb+namelist_name+ & value
          'column': printed output as namelist_name value
          'dict': as two python dictionary object (namelistname, value and namelistname, sectionname)
      namelist= namelist_like file to retrieve values
    >>> get_namelist_vars('geogrid,dic','/home/lluis/etudes/domains/medic950116/namelist.wps')
    {'e_sn': '32, 97,', 'stand_lon': '0.', 'opt_geogrid_tbl_path': "'./'", 'geog_data_path': "'/home/lluis/DATA/WRF/geog'", 'pole_lat': '90.0', 'ref_lat': '35.0,', 'map_proj': "'lat-lon',", 'parent_id': '1, 1,', 'geog_data_res': "'2m','2m',", 'e_we': '32, 112,', 'dx': '0.35,', 'dy': '0.35,', 'parent_grid_ratio': '1, 3,', 'pole_lon': '0.0', 'ref_lon': '20,', 'j_parent_start': '1, 17,', 'i_parent_start': '1, 31,'}
    """

    fname = 'get_namelist_vars'

# List of characters which split pairs name/value
    valuessep = ['=']
# List of characters which indicate a comment
    commentchars = ['#']

    if values == 'h':
        print fname + '_____________________________________________________________'
        print get_namelist_vars.__doc__
        quit()

    expectargs = '[sectionname],[kout]'
    check_arguments(fname,values,expectargs,',')

    secname = values.split(',')[0]
    kout = values.split(',')[1]

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

    ncml = open(namelist, 'r')

    sections = {}
    namelistvals = {}
    namelistsecs = {}
    sectionnames = []
    namessec = []
    allnames = []
    namelistvalssec = {}
    namelistsecssec = {}
    nmlname = ''
    sectionname = ''

    for line in ncml:
        linevals = reduce_spaces(line)
        Nvals = len(linevals)

        if Nvals >= 1 and linevals[0][0:1] == '&':
            if len(sectionnames) > 1:
                sections[sectionname] = namessec

            sectionname = linevals[0][1:len(linevals[0])+1]
#            print '    ' + fname + ": new section '" + sectionname + "' !!!"
            sectionnames.append(sectionname)
            namessec = []
            nmlname = ''
        elif Nvals >= 1 and not searchInlist(commentchars,linevals[0][0:1]):
            if Nvals >= 3 and searchInlist(valuessep,linevals[1]):
                nmlname = linevals[0]
                nmlval = numVector_String(linevals[2:Nvals],' ')
            elif Nvals == 1:
                for valsep in valuessep:
                    if linevals[0].find(valsep) != -1:
                        nmlname = linevals[0].split(valsep)[0]
                        nmlval = linevals[0].split(valsep)[1]
                        break
            elif Nvals == 2:
                for valsep in valuessep:
                    if linevals[0].find(valsep) != -1:
                        nmlname = linevals[0].split(valsep)[0]
                        nmlval = linevals[1]
                        break
                    elif linevals[1].find(valsep) != -1:
                        nmlname = linevals[0]
                        nmlval = linevals[1].split(valsep)[0]
                        break
            else:
                print warnmsg
                print '  ' + fname + ': wrong number of values', Nvals,              \
                  'in the namelist line!'
                print '    line: ',line
                print '    line values:',linevals
#                quit(-1)

            namelistvals[nmlname] = nmlval
            namelistsecs[nmlname] = sectionname

            namessec.append(nmlname)
            allnames.append(nmlname)

    if len(sectionname) > 1:
        sections[sectionname] = namessec

    if secname != 'all':
        if not searchInlist(sections.keys(),secname):
            print errormsg
            print '  ' + fname + ": section '" + values + "' does not exist !!"
            print '    only found:',sectionnames
            quit(-1)

        namestouse = []
        for nml in allnames:
            for nnml in sections[secname]:
                namelistvalssec[nnml] = namelistvals[nnml]
                namelistsecssec[nnml] = secname
                if nml == nnml: namestouse.append(nml)
    else:
        namestouse = allnames
        namelistvalssec = namelistvals
        namelistsecssec = namelistsecs

    if kout == 'tex3':
        ofile='get_namelist_vars_3col.tex'
        fobj = open(ofile, 'w')

        vals = namestouse
        Nvals = len(vals)
        Nvals3 = int(Nvals/3)
        latextab = '\\begin{center}\n\\begin{tabular}{lclclc}\n'
        for il in range(2):
            latextab = latextab + '{\\bfseries{name}} & {\\bfseries{value}} & '
        latextab= latextab+ '{\\bfseries{name}} & {\\bfseries{value}} \\\\ \\hline\n'

        if np.mod(Nvals,3) != 0:
            Nvals0 = Nvals - np.mod(Nvals,3)
        else:
            Nvals0 = Nvals

        for ival in range(0,Nvals0,3):
            line = ''
            print '  ival:',ival
            for il in range(2):
                line = line + '\\verb+' + vals[ival+il] + '+ & ' +                   \
                   namelistvalssec[vals[ival+il]].replace('_','\\_') +' & '
            line = line + '\\verb+' + vals[ival+2] + '+ & ' +                       \
               namelistvalssec[vals[ival+2]].replace('_','\\_') + ' \\\\\n'
            latextab = latextab + line

        latextab = latextab + '%not multiple of three!!!!\n'
        print 'mod:', np.mod(Nvals,3),Nvals0,Nvals
        if np.mod(Nvals,3) != 0:
            ival = Nvals0
            line = ''
            for il in range(np.mod(Nvals,3)):
                print 'ival:',ival + il
                line = line + '\\verb+' + vals[ival+il] + '+ & ' +                   \
                      namelistvalssec[vals[ival+il]].replace('_','\\_') + ' & '
            for il in range(2-np.mod(Nvals,3)):
                line = line + ' & & '
            latextab = latextab + line + ' & \\\\\n'
        latextab = latextab + '\\end{tabular}\n\\end{center}\n'

#        print latextab
        fobj.write(latextab)
 
        fobj.close()
        print fname + ": successful writen '" + ofile + "' LaTeX tale file !!"

        return 
    elif kout == 'column':
        for dictv in namestouse:
            print dictv + ' = ' + namelistvalssec[dictv]
        return
    elif kout == 'dict':
        return namelistvalssec, namelistsecssec
    else:
        print errormsg
        print '  ' + fname + ": kind of output '" + kout + "' not ready!!!"
        quit(-1)

    return

def WRF_d0Nref(values, geofold):
    """ Function for the generation of an extra WRF domain from a given one 
    field_stats(values, ncfile, varn)
      [values]= [namelist],[minLON],[minLAT],[maxLON],[maxLAT],[Ndomref],[factor],[geogres]
        [namelist]: namelist.wps to use
        [minLON],[minLAT]: minimum longitude and latitude
        [maxLON],[maxLAT]: maximum longitude and latitude
        [Ndomref]: domain reference number
        [factor]: resolution factor with respect d01
        [geogres]: which topography resolution to use
      [geofold]= folder where to found the geo_em.d0[N-1].nc file
    """

    fname='WRF_d0Nref'

    if values == 'h':
        print fname + '_____________________________________________________________'
        print WRF_d0Nref.__doc__
        quit()

    expectargs = '[namelist],[minLON],[minLAT],[maxLON],[maxLAT],[Ndomref],[factor],'
    expectargs = expectargs + '[geogres]'
 
    check_arguments(fname,values,expectargs,',')

    namelist = values.split(',')[0]
    minLON = np.float(values.split(',')[1])
    minLAT = np.float(values.split(',')[2])
    maxLON = np.float(values.split(',')[3])
    maxLAT = np.float(values.split(',')[4])
    Ndomref = int(values.split(',')[5])
    factor = int(values.split(',')[6])
    geogres = values.split(',')[7]

    onml = 'namelist.wps_d0' + str(Ndomref+1)

# Namelist variables to which new value is a copy of the last one
    repeatvalue = ['start_date' ,'end_date']
# Namelist variables to which value is increased by 1 (single value)
    add1 = ['max_dom']
# Namelist variables to which new value is an increase by 1 of the last one
    add1value = ['parent_id']

#    ncfile = geofold + '/geo_em.d' + zeros(Ndomref,2) + '.nc'
    ncfile = geofold + '/geo_em.d0' + str(Ndomref) + '.nc'

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

    ncobj = NetCDFFile(ncfile, 'r')

    nlon='XLONG_M'
    nlat='XLAT_M'

    lonobj = ncobj.variables[nlon]
    latobj = ncobj.variables[nlat]

    lon = lonobj[0,:,:]
    lat = latobj[0,:,:]

    ncobj.close()

    diffmin = np.sqrt((lon - minLON)**2. + (lat - minLAT)**2.)
    diffmax = np.sqrt((lon - maxLON)**2. + (lat - maxLAT)**2.)

    posdiffmin = index_mat(diffmin,np.min(diffmin))
    posdiffmax = index_mat(diffmax,np.min(diffmax))
#    print 'min vals. min:',np.min(diffmin),'max:',np.min(diffmax)

#    print 'min:', minLON, ',', minLAT, ':', posdiffmin, '.', lon[tuple(posdiffmin)], \
#      ',', lat[tuple(posdiffmin)]
#    print 'max:', maxLON, ',', maxLAT, ':', posdiffmax, '.', lon[tuple(posdiffmax)], \
#      ',', lat[tuple(posdiffmax)]

    print '  ' + fname +': coordinates of the SW corner:', posdiffmin, ' NE corner:',\
      posdiffmax

    nmlobj = open(namelist, 'r')
    onmlobj = open(onml, 'w')
    for line in nmlobj:
        listline = reduce_spaces(line)
        nvals = len(listline)
        if nvals > 0:
            if searchInlist(add1,listline[0]):
                nmlval = int(listline[nvals-1].replace(',','')) + 1
                onmlobj.write(' ' + listline[0] + ' = ' + str(nmlval) + ',\n')
            elif searchInlist(add1value,listline[0]):
                if Ndomref == 1:
                    nmlval = 1
                else:
                    nmlval = int(listline[nvals-1].replace(',','')) + 1   
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1], ' ') + ' ' +str(nmlval)+',\n')
            elif searchInlist(repeatvalue,listline[0]):
                nmlval = listline[nvals-1]
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+'\n')
            elif listline[0] == 'parent_grid_ratio':
                nmlval = factor
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            elif listline[0] == 'i_parent_start':
                nmlval = posdiffmin[1] - 10
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            elif listline[0] == 'j_parent_start':
                nmlval = posdiffmin[0] - 10
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            elif listline[0] == 'e_we':
                nmlval = (posdiffmax[1] - posdiffmin[1] + 20)*factor
                nmlval = nmlval + np.mod(nmlval,factor) + 1
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            elif listline[0] == 'e_sn':
                nmlval = (posdiffmax[0] - posdiffmin[0] + 20)*factor
                nmlval = nmlval + np.mod(nmlval,factor) + 1
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            elif listline[0] == 'geog_data_res':
                nmlval = "'" + geogres + "'"
                onmlobj.write(' ' + listline[0] + ' = ' +                            \
                  numVector_String(listline[2:nvals+1],' ') + ' ' + str(nmlval)+',\n')
            else:
                onmlobj.write(line)
        else:
            onmlobj.write('\n')

    nmlobj.close()
    onmlobj.close()

    print fname + ": successfully written new namelist '" + onml + "' !!"

    return

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 wdismean(pos,val4):
    """ Function to compute the mean value weighted to its 4 distances
      pos= x,y position within the 'square'
      val4= values at the four vertexs (2,2)
    >>>position = [0.005,0.005]
    >>>val4 = np.zeros((2,2), dtype=np.float)
    >>>val4 = np.arange(4).reshape(2,2)
    0.035707955234
    """
    fname = 'wdismean'

    dist = np.zeros((2,2), dtype=np.float)    
    dist[0,0] = np.sqrt(pos[0]*pos[0]+pos[1]*pos[1])
    dist[0,1] = np.sqrt(pos[0]*pos[0]+(1.-pos[1])*(1.-pos[1]))
    dist[1,0] = np.sqrt((1.-pos[0])*(1.-pos[0])+pos[1]*pos[1])
    dist[1,1] = np.sqrt((1.-pos[0])*(1.-pos[0])+(1.-pos[1])*(1.-pos[1]))

    if np.min(dist) == 0.:
        pos0 = index_mat(dist, 0.)
        dist[:,:] = 0.
        posmean = val4[pos0[0], pos0[1]]
    else:
        totdist = np.sum(1./dist)
        posmean = np.sum(val4.flatten()/dist.flatten())/totdist

    return posmean

def read_ASCIIlist(filen):
    """ Function to build a list from an ASCII lines file
      filen= name of the ASCII file
    """
    import os
    fname = 'read_ASCIIlist'

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

    of = open(filen, 'r')
    filevals = []
    for linef in of:
        filevals.append(linef.replace('\n',''))
    of.close()

    return filevals

def radial_points(angle,dist):
    """ Function to provide a number of grid point positions for a given angle
      angle= desired angle (rad)
      dist= number of grid points
    >>> radial_points(0.5*np.pi,5)
    [[  6.12323400e-17   1.00000000e+00]
     [  1.22464680e-16   2.00000000e+00]
     [  1.83697020e-16   3.00000000e+00]
     [  2.44929360e-16   4.00000000e+00]
     [  3.06161700e-16   5.00000000e+00]]
    """
    fname = 'radial_points'
    radpos = np.zeros((dist,2), dtype=np.float)

    for ip in range(dist):
        radpos[ip,0] = (ip+1.)*np.cos(angle)
        radpos[ip,1] = (ip+1.)*np.sin(angle)

    return radpos

def read_SQradii(filename='SQradii.dat', Npt=-1, maxradii=None):
    """ Function to read the outcome of the function `squared_radial'
      filename= Name of the outcome file of the function 
      Npt= Number of points on each grid
      maxradii= allowed maximum radii
      'SQradii.dat': file with the radii and paris of integers after ran squared_radial(1000)
    """
    fname = 'read_SQradii'

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

    infile = open(filename, 'r')

    gridradii = {}
    Nradii = {}
    radii = []
    for line in infile:
        values = values_line(line, ' ', ['[', ']'])
        rad = np.float(values[0])
# Limited by a given radii
        if maxradii is not None and rad > maxradii: break
        Nrad = int(values[2])
        couples = []
        Nradgood = 0
        if Nrad > 1:
           Ncouples = 0
           for ic in range(Nrad):
               ix = int(values[3 + ic*2])
               iy = int(values[3 + ic*2 + 1])
               if ix <= Npt and iy <= Npt:
                   couples.append(np.array([ix,iy]))
                   Nradgood = Nradgood + 1
               elif ix > Npt and iy > Npt: break
        else:
             ix = int(values[3])
             iy = int(values[4])
             if ix <= Npt and iy <= Npt:
                 couples.append(np.array([ix,iy]))
                 Nradgood = Nradgood + 1
             elif ix > Npt and iy > Npt: break

        if Nradgood > 0:
            Nradii[rad] = Nradgood
            finalcouples = np.zeros((Nradgood,2), dtype=int)
            for ic in range(Nradgood):
                finalcouples[ic,:] = couples[ic]
            gridradii[rad] = finalcouples

#    for rad in sorted(gridradii.keys()):
#        print rad, gridradii[rad]

    return gridradii

def grid_combinations(x,y):
    """ Function to provide all the possible grid points combination for a given pair of values
      x,y= pair of grid points
    >>>grid_combinations(1,2)
    [[1, 2], [-1, 2], [-1, -2], [1, -2], [2, 1], [-2, 1], [-2, -1], [2, -1]]
    """
    fname = 'grid_combinations'

    gridcomb = []
    gridcomb.append([x,y])
    if x != 0: gridcomb.append([-x,y])
    if x != 0 and y != 0: gridcomb.append([-x,-y])
    if y != 0: gridcomb.append([x,-y])

    if x != y:
        gridcomb.append([y,x])
        if y != 0: gridcomb.append([-y,x])
        if x != 0 and y != 0: gridcomb.append([-y,-x])
        if x != 0: gridcomb.append([y,-x])

    return gridcomb

def radii_points(xpos,ypos,Lrad,Nang,dx,dy):
    """ Function to provide the grid points for radial cross sections, by angle and radii and `squared' radii
      xpos= x position of the center
      ypos= y position of the center
      Lrad= length of the maximum radi in grid points
      Nang= number of equivalent angles
      dx= x size of the domain of values
      dy= y size of the domain of values
      >>> radpos, SQradpos = radii_points(12,12,1,8,25,25)
      radpos:
      [[[ 13.          12.        ]]
       [[ 12.70710678  12.70710678]]
       [[ 12.          13.        ]]
       [[ 11.29289322  12.70710678]]
       [[ 11.          12.        ]]
       [[ 11.29289322  11.29289322]]
       [[ 12.          11.        ]]
       [[ 12.70710678  11.29289322]]]
      SQradpos:
      {1.0: [[13, 12], [11, 12], [12, 13], [12, 11]], 1.41421356237: [[13, 13], [11, 13], [11, 11], [13, 11]]}
    """
    fname = 'radii_points'

    radiipos = np.zeros((Nang, Lrad, 2), dtype = np.float)
    SQradiipos = {}

# Getting the range of radii:
    pairsSQradii = read_SQradii(Npt=Lrad)
    SQradii = pairsSQradii.keys()
    Nradii = len(SQradii)

# Angle loop
    for ia in range(Nang):
        angle = 2.*np.pi*ia/Nang
        posangle = np.zeros((Lrad,2), dtype=np.float)

# Points at that given angle
        pangle = radial_points(angle,Lrad)
        posangle[:,0] = pangle[:,0] + xpos
        posangle[:,1] = pangle[:,1] + ypos

        for ip in range(Lrad):
            iipos = int(posangle[ip,0])
            jjpos = int(posangle[ip,1])

# Creation of a matrix with all the grid points at each time-step
            if iipos >= 0 and iipos < dx-1 and jjpos >= 0 and jjpos < dy-1:
                radiipos[ia,ip,0] = posangle[ip,0]
                radiipos[ia,ip,1] = posangle[ip,1]

# SQ radii values

# Radii loop (avoiding 0)
    pairswithin = {}
    for ir in range(1,Nradii):
        pairs = pairsSQradii[SQradii[ir]]

# pairs loop
        Npairs = len(pairs.flatten())/2
        pairsw = []
        for ip in range(Npairs):
            if Npairs == 0: 
                break
            else:
                allpairs = grid_combinations(pairs[ip,0],pairs[ip,1])                
            for iap in range(len(allpairs)):
                ppos = allpairs[iap]
                iipos = xpos + ppos[0]
                jjpos = ypos + ppos[1]

# Creation of a matrix with all the grid points at each time-step
                if iipos >= 0 and iipos < dx and jjpos >= 0 and jjpos < dy:
                    pairsw.append([iipos,jjpos])

        SQradiipos[SQradii[ir]] = pairsw

    return radiipos, SQradiipos

def hor_weight_int(SWval, SEval, NWval, NEval, xpos, ypos):
    """ Function to weighted by distance interpolate a horizontal value from it closest 4 neighbourgs 
      field= a 4 points representing the environment of a given point
      xpos, ypos: position to which one want to interpolate (relative to 0,0 at the corner SW)
    >>> hor_weight_int(1., 1., 1., 2., 0.75, 0.75)
    1.44888128193
    """
    fname = 'hor_weight_int'

    vals = np.array([SWval, SEval, NWval, NEval])
    print vals
    if xpos > 1.: 
        print errormsg
        print '  ' + fname + ': Wrong position to interpolate! (',xpos,',',ypos,     \
          ') must be within (1,1) !!'
        quit(-1)

    vertexs = np.zeros((4,2), dtype=np.float) 
    vertexs[0,0] = 0.
    vertexs[0,1] = 0.
    vertexs[1,0] = 1.
    vertexs[1,1] = 0.
    vertexs[2,0] = 0.
    vertexs[2,1] = 1.
    vertexs[3,0] = 1.
    vertexs[3,1] = 1.
    print vertexs

    dist = np.sqrt((vertexs[:,0]-xpos)**2 + (vertexs[:,1]-ypos)**2)
    print dist
    done = False
    for id in range(4):
        if dist[id] == 0.: 
            intval = vals[id]
            done = True

    if not done: intval = np.sum(vals/dist)/np.sum(1./dist)

    return intval

def rmNOnum(string):
    """ Removing from a string all that characters which are not numbers
    # From: http://stackoverflow.com/questions/4289331/python-extract-numbers-from-a-string
    """
    fname = 'rmNOnum'

    newstr = str(re.findall("[-+]?\d+[\.]?\d*", string)[0])

    return newstr

def values_fortran_fmt(lin,fFMT):
    """ Function to give back the values of an ASCII line according to its fortran printing format
      lin= ASCII line
      fFMT= list with the fortran FORMAT formats
    forline = 'Natchitoches (RGNL)        1 11 0011  ( 31.733, -93.100) (  28, 157) ( 31.761, -93.113)   41.2 meters'
    formats = ['A26', 'I2', 'I3', 'A6', 'A2', 'F7.3', 'A1', 'F8.3', 'A3', 'I4', 'A1', 'I4',
      'A3', 'F7.3', 'A1', 'F8.3', 'A2', 'F6.1', 'A7']
    >>> values_fortran_fmt(forline, formats)
    ['Natchitoches (RGNL)        ', 1, 11, ' 0011  ', ' ( ', 31.733, ', ', -93.1, ') ( ', 28, ', ', 157, ')
      ( ', 31.761, ', ', -93.113, ')  ', 41.2, ' meters']
    """
    fname = 'values_fortran_fmt'

    afmts = ['A', 'D', 'F', 'I']

    if lin == 'h':
        print fname + '_____________________________________________________________'
        print values_fortran_fmt.__doc__
        quit()

    fvalues = []
    ichar=0
    ival = 0
    for ft in fFMT:
        Lft = len(ft)
        if ft[0:1] == 'A' or ft[0:1] == 'a':
            Lfmt = int(ft[1:Lft+1])
            fvalues.append(lin[ichar:ichar+Lfmt+1])
            ichar = ichar + Lfmt
        elif ft[0:1] == 'F' or ft[0:1] == 'f':
            if ft.find('.') != -1:
                Lft = len(ft.split('.')[0])
            Lfmt = int(ft[1:Lft])
            fvalues.append(np.float(rmNOnum(lin[ichar:ichar+Lfmt+1])))
            ichar = ichar + Lfmt
        elif ft[0:1] == 'D' or ft[0:1] == 'd':
            if ft.find('.') != -1:
                Lft = len(ft.split('.')[0])
            Lfmt = int(ft[1:Lft])
            fvalues.append(np.float64(rmNOnum(lin[ichar:ichar+Lfmt+1])))
            ichar = ichar + Lfmt
        elif ft[0:1] == 'I' or ft[0:1] == 'i':
            Lfmt = int(ft[1:Lft+1])
            fvalues.append(int(rmNOnum(lin[ichar:ichar+Lfmt+1])))
            ichar = ichar + Lfmt
        elif ft.find('X') != -1 or ft.find('x') != -1:
            print '    ' + fname + ': skipping space'
            ichar = ichar + int(rmNOnum(ft))
        else:
            print errormsg
            print '  ' + fname + ": format '" + ft[0:1] + "' not ready!!"
            print '    Available formats:',afmts
            quit(-1)

        ival = ival + 1

    return fvalues

def reduce_last_spaces(string):
    """ Function to reduce the last right spaces from a string
      string= string to remove the spaces at the righ side of the string
    >>> reduce_last_spaces('Hola Lluis      ')
    'Hola Lluis'
    """
    fname = 'reduce_last_spaces'

    Lstring = len(string)

    firstspace = -1
    for istr in range(Lstring-1,0,-1):
        if string[istr:istr+1] == ' ':
            firstspace = istr
        else:
            break

    if firstspace == -1:
        print errormsg
        print '  ' + fname + ": string '" + string + "' is not ended by spaces!!"
        print '    char. number of last character:',ord(string[Lstring:Lstring+1])
        quit(-1)

    newstr = string[0:firstspace]

    return newstr

def squared_radial(Npt):
    """ Function to provide the series of radii as composite of pairs (x,y) of gid cells
      Npt= largest amount of grid points on x and y directions
    >>> squared_radial(5)
    [[0.0, 0, 0], [1.0, 1, 0], [1.4142135623730951, 1, 1], [2.0, 2, 0], [2.2360679774997898, 2, 1], [2.8284271247461903, 2, 2], [3.0, 3, 0], [3.1622776601683795, 3, 1], [3.6055512754639891, 3, 2], [4.0, 4, 0], [4.1231056256176606, 4, 1], [4.2426406871192848, 3, 3], [4.4721359549995796, 4, 2], [5.0, array([4, 3]), array([5, 0])], [5.0990195135927845, 5, 1], [5.3851648071345037, 5, 2], [5.6568542494923806, 4, 4], [5.8309518948453007, 5, 3], [6.4031242374328485, 5, 4], [7.0710678118654755, 5, 5]]
    """
    fname = 'squared_radial'

    radii = []
    gridradii = {}
    singradii = []
    Nradii = {}

# Computing all radii
##
    for i in range(Npt+1):
        if np.mod(i,100) == 0: print 'done: ',i
        for j in range(i+1):
#            ir = np.float(i)
#            jr = np.float(j)
            ir = np.float64(i)
            jr = np.float64(j)
            rad = np.sqrt(ir*ir + jr*jr)
            if not searchInlist(singradii, rad):
                singradii.append(rad)
                gridradii[rad] = np.array([i,j], dtype=int)
                Nradii[rad] = 1
            else:
#                print warnmsg
#                print '  ' + fname + ': repeated radii',rad,'!!'
#                print '    Previous values:', gridradii[rad]

                oldvalues = gridradii[rad]
                Nradii[rad] = Nradii[rad] + 1
                newvalues = np.zeros((Nradii[rad],2), dtype=int)
                if Nradii[rad] == 2:
                    newvalues[0,:] = oldvalues
                    Nstart = 1
                else:
                    Nstart = 0
                for Nr in range(Nstart,Nradii[rad]-1):
                    newvalues[Nr,:] = oldvalues[Nr,:]
                newvalues[Nradii[rad]-1,:] = np.array([i,j], dtype=int)
                gridradii[rad] = newvalues

# Sorting values
##
## Only required to downwrite the output file
#    radiif = open('SQradii.dat', 'w')
#    radii = []
#    prevrad = 0.
#    for rad in sorted(gridradii.keys()):
#        dradii = rad - prevrad
## Check fo repeated values
#        radii.append([rad, gridradii[rad][0], gridradii[rad][1]])
#        radiif.write(str(rad) + ' ' + str(dradii) + ' ' + str(Nradii[rad]) + ' [ '+  \
#          numVector_String(gridradii[rad].flatten(), ' ') + ' ] \n')
#        prevrad = rad.copy()

#    radiif.close()
    return gridradii

def fillvalue_kind(vartype, fillval0):
    """ Function to provide the fiven fill_Value for a kind of variable
      [vartype]= type of the variable
      [fillval0]= value to take as fill value, 'std' for the standard value
           Float = 1.e20
           Character = '-'
           Integer = -99999
           Integer16 = -99999
           Float64 = 1.e20
           Integer32 = -99999
    """
    fname = 'fillvalue_kind'

    if vartype == type('c'):
        if fillval0 == 'std':
            fillval = fillvalueC
        else:
            fillval = fillval0
    elif vartype == type(int(1)):
        if fillval0 == 'std':
            fillval = fillValueI
        else:
            fillval = int(fillval0)
    elif vartype == type(np.int16(1)):
        if fillval0 == 'std':
            fillval = np.int(9999999)
        else:
            fillval = np.int16(fillval0)
    elif vartype == type(np.int32(1)):
        if fillval0 == 'std':
            fillval = np.int32(fillValueI)
        else:
            fillval = np.int32(fillval0)
    elif vartype == type(np.float(1.)):
        if fillval0 == 'std':
            fillval = np.float(fillValueF)
        else:
            fillval = np.float(fillval0)
    elif vartype == type(np.float32(1.)):
        if fillval0 == 'std':
            fillval = np.float32(fillValueF)
        else:
            fillval = np.float32(fillval0)
    elif vartype == type(np.float64(1.)):
        if fillval0 == 'std':
            fillval = np.float64(fillValueF)
        else:
            fillval = np.float64(fillval0)
    else:
        print errormsg
        print '  ' + fname + ': variable type', vartype, 'not ready !!'
        quit(-1)

    return fillval

def varval_ind(var, ind):
    """ Function to provide a value from a variable providing a value for each index
      var= variable
      ind= list of indices for each dimension of the variable
    """
    fname = 'varval_ind'
    if len(var.shape) != len(ind):
        print errormsg
        print '  ' + fname + ': different shape:', var.shape,'and indices:', ind, '!!'
        quit(-1)

    varslice=[]
    for iv in ind:
        varslice.append(iv)

    val = varval_ind[slice(varslice)]

    return val

def CoarselonlatFind(ilon, ilat, lonv, latv, per):
    """ Function to search a given value from a coarser version of the data
      ilon, ilat= original 2D matrices with the longitudes and the latitudes
      lonv, latv= longitude and latitude to find
      per= period (as fraction over 1) of the fractions of the original grid to use to explore
      >>> npt=100
      >>> lonval0 = np.arange(0.,npt*1.,1.)
      >>> latval0 = np.arange(0.,npt*1.,1.)
      >>> lonval = np.zeros((npt,npt), dtype=np.float)
      >>> latval = np.zeros((npt,npt), dtype=np.float)

      >>> for i in range(npt):
      >>>     lonval[:,i] = lonval0[i]
      >>>     latval[i,:] = latval0[i]
      >>> CoarselonlatFind(lonval, latval, 5.25, 5.25, .1)
      (array([5, 5]), 0.35355339059327379)
    """
    fname = 'CoarselonlatFind'

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

    nlon = np.min(ilon)
    xlon = np.max(ilon)
    nlat = np.min(ilat)
    xlat = np.max(ilat)

    if lonv < nlon or lonv > xlon:
        print errormsg
        print '  ' + fname + ': longitude outside data range!!'
        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
        quit(-1) 
    if latv < nlat or latv > xlat:
        print errormsg
        print '  ' + fname + ': latitude outside data range!!'
        print '    given value:', latv,'outside [',nlat,',',xlat,']'
        quit(-1) 

    fracx = int(dx*per)
    fracy = int(dy*per)

    fraclon = ilon[::fracy,::fracx]
    fraclat = ilat[::fracy,::fracx]
    fracdx = fraclon.shape[1]
    fracdy = fraclon.shape[0]

#    print 'fraclon _______'
#    print fraclon

#    print 'fraclat _______'
#    print fraclat

# Fraction point
    difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
    mindifffracLl = np.min(difffraclonlat)
    ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)

#    print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
#    print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',             \
#      fraclat[ilatlonfrac[0],ilatlonfrac[1]]
#    print 'values lon, lat:', lonv, latv
      
# Providing fraction range
    if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                            \
      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
        if ilatlonfrac[0] > 0: 
            iybeg = (ilatlonfrac[0]-1)*fracy
            iyend = ilatlonfrac[0]*fracy+1
        else:
            iybeg = 0
            iyend = fracy+1
        if ilatlonfrac[1] > 0: 
            ixbeg = (ilatlonfrac[1]-1)*fracx
            ixend = ilatlonfrac[1]*fracx+1
        else:
            ixbeg = 0
            ixend = fracx+1
    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
      fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
#        ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
#        ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
        if ilatlonfrac[0] > 0: 
            iybeg = (ilatlonfrac[0]-1)*fracy
            iyend = ilatlonfrac[0]*fracy+1
        else:
            iybeg = 0
            iyend = fracy+1
        if ilatlonfrac[1] < fracdx:
            if ilatlonfrac[1] != 0:
                ixbeg = (ilatlonfrac[1]-1)*fracx
                ixend = ilatlonfrac[1]*fracx+1
            else:
                ixbeg = 0
                ixend = fracx+1
        else:
            ixbeg = fracdx*fracx
            ixend = dx+1
    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                           \
      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
        if ilatlonfrac[0] < fracdy: 
            if ilatlonfrac[0] != 0:
                iybeg = (ilatlonfrac[0]-1)*fracy
                iyend = ilatlonfrac[0]*fracy+1
            else:
                iybeg = 0
                iyend = fracy+1
        else:
            iybeg = fracdy*fracy
            iyend = dy+1
        if ilatlonfrac[1] < fracdx: 
            if ilatlonfrac[1] != 0:
                ixbeg = (ilatlonfrac[1]-1)*fracx
                ixend = ilatlonfrac[1]*fracx+1
            else:
                ixbeg = 0
                ixend = fracx+1
        else:
            ixbeg = fracdx*fracx
            ixend = dx+1
    elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                          \
      fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
#        ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
#        ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
        if ilatlonfrac[0] < fracdy: 
            if ilatlonfrac[0] != 0:
                iybeg = (ilatlonfrac[0]-1)*fracy
                iyend = ilatlonfrac[0]*fracy+1
            else:
                iybeg = 0
                iyend = fracy+1
        else:
            iybeg = fracdy*fracy
            iyend = dy+1
        if ilatlonfrac[1] > 0: 
            ixbeg = (ilatlonfrac[1]-1)*fracx
            ixend = ilatlonfrac[1]*fracx+1
        else:
            ixbeg = 0
            ixend = fracx+1

#    ibeg = [iybeg, ixbeg]
#    iend = [iyend, ixend]
#    print 'find window; beginning', ibeg, 'end:', iend
    lon = ilon[iybeg:iyend,ixbeg:ixend]
    lat = ilat[iybeg:iyend,ixbeg:ixend]

#    print 'lon _______'
#    print lon
#    print 'lat _______'
#    print lat

# Find point
    difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
    mindiffLl = np.min(difflonlat)
    ilatlon = index_mat(difflonlat, mindiffLl)

    ilatlon[0] = ilatlon[0] + iybeg
    ilatlon[1] = ilatlon[1] + ixbeg

#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
#    quit()

    return ilatlon, mindiffLl

def CoarselonlatFindAll(onc, ilon, ilat, longvs, latvs, per, frac, out):
    """ Function to search all values from a coarser version of the data
      onc= netCDF object from an existing file
      ilon, ilat= original 2D matrices with the longitudes and the latitudes
      lonv, latv= longitude and latitude to find
      per= period (as fraction over 1) of the fractions of the original grid to use to explore
      frac= frequency of points at which syncronize file with calculations
      out= True/False if outcome is required
      >>> npt=100
      >>> lonval0 = np.arange(0.,npt*1.,1.)
      >>> latval0 = np.arange(0.,npt*1.,1.)
      >>> lonval = np.zeros((npt,npt), dtype=np.float)
      >>> latval = np.zeros((npt,npt), dtype=np.float)
      >>> lonvalues = np.arange(0.25,10.,0.25)
      >>> latvalues = np.arange(0.25,10.,0.25)

      >>> for i in range(npt):
      >>>     lonval[:,i] = lonval0[i]
      >>>     latval[i,:] = latval0[i]
      >>> CoarselonlatFindAll(onewnc, lonval, latval, lovalues, latvalues, .1, 100)
      (array([5, 5]), 0.35355339059327379)
    """
    fname = 'CoarselonlatFindAll'

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

    nlon = np.min(ilon)
    xlon = np.max(ilon)
    nlat = np.min(ilat)
    xlat = np.max(ilat)

    Nallpts = len(longvs)

    fracx = int(dx*per)
    fracy = int(dy*per)

    fraclon = ilon[::fracy,::fracx]
    fraclat = ilat[::fracy,::fracx]
    fracdx = fraclon.shape[1]
    fracdy = fraclon.shape[0]

#    print 'fraclon _______'
#    print fraclon

#    print 'fraclat _______'
#    print fraclat

    if out:
        ilatlon = np.zeros((2,Nallpts), dtype=int)
        mindiffLl = np.zeros((Nallpts), dtype=np.float)

    iptpos=0
    for iv in range(Nallpts):
        lonv = longvs[iv]
        latv = latvs[iv]

        if lonv < nlon or lonv > xlon:
            print errormsg
            print '  ' + fname + ': longitude outside data range!!'
            print '    given value:', lonv,'outside [',nlon,',',xlon,']'
            quit(-1) 
        if latv < nlat or latv > xlat:
            print errormsg
            print '  ' + fname + ': latitude outside data range!!'
            print '    given value:', latv,'outside [',nlat,',',xlat,']'
            quit(-1) 

# Fraction point
        difffraclonlat = np.sqrt((fraclon-lonv)**2. + (fraclat-latv)**2.)
        mindifffracLl = np.min(difffraclonlat)
        ilatlonfrac = index_mat(difffraclonlat, mindifffracLl)

#        print 'values lon, lat:', lonv, latv
#        print 'mindifffracLl:', mindifffracLl, 'ilatlonfrac:', ilatlonfrac
#        print 'frac lon, lat:', fraclon[ilatlonfrac[0],ilatlonfrac[1]], ',',         \
#          fraclat[ilatlonfrac[0],ilatlonfrac[1]]
      
# Providing fraction range
        if fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                        \
          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]-1]
#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]]
            if ilatlonfrac[0] > 0: 
                iybeg = (ilatlonfrac[0]-1)*fracy
                iyend = ilatlonfrac[0]*fracy+1
            else:
                iybeg = 0
                iyend = fracy+1
            if ilatlonfrac[1] > 0: 
                ixbeg = (ilatlonfrac[1]-1)*fracx
                ixend = ilatlonfrac[1]*fracx+1
            else:
                ixbeg = 0
                ixend = fracx+1
        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
          fraclat[ilatlonfrac[0],ilatlonfrac[1]] >= latv:
#            ifracbeg = [ilatlonfrac[0]-1, ilatlonfrac[1]]
#            ifracend = [ilatlonfrac[0], ilatlonfrac[1]+1]
            if ilatlonfrac[0] > 0: 
                iybeg = (ilatlonfrac[0]-1)*fracy
                iyend = ilatlonfrac[0]*fracy+1
            else:
                iybeg = 0
                iyend = fracy+1
            if ilatlonfrac[1] < fracdx:
                if ilatlonfrac[1] != 0:
                    ixbeg = (ilatlonfrac[1]-1)*fracx
                    ixend = ilatlonfrac[1]*fracx+1
                else:
                    ixbeg = 0
                    ixend = fracx+1
            else:
                ixbeg = fracdx*fracx
                ixend = dx+1
        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] < lonv and                       \
          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]]
#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]+1]
            if ilatlonfrac[0] < fracdy: 
                if ilatlonfrac[0] != 0:
                    iybeg = (ilatlonfrac[0]-1)*fracy
                    iyend = ilatlonfrac[0]*fracy+1
                else:
                    iybeg = 0
                    iyend = fracy+1
            else:
                iybeg = fracdy*fracy
                iyend = dy+1
            if ilatlonfrac[1] < fracdx: 
                if ilatlonfrac[1] != 0:
                    ixbeg = (ilatlonfrac[1]-1)*fracx
                    ixend = ilatlonfrac[1]*fracx+1
                else:
                    ixbeg = 0
                    ixend = fracx+1
            else:
                ixbeg = fracdx*fracx
                ixend = dx+1
        elif fraclon[ilatlonfrac[0],ilatlonfrac[1]] >= lonv and                      \
          fraclat[ilatlonfrac[0],ilatlonfrac[1]] < latv:
#            ifracbeg = [ilatlonfrac[0], ilatlonfrac[1]-1]
#            ifracend = [ilatlonfrac[0]+1, ilatlonfrac[1]]
            if ilatlonfrac[0] < fracdy: 
                if ilatlonfrac[0] != 0:
                    iybeg = (ilatlonfrac[0]-1)*fracy
                    iyend = ilatlonfrac[0]*fracy+1
                else:
                    iybeg = 0
                    iyend = fracy+1
            else:
                iybeg = fracdy*fracy
                iyend = dy+1
            if ilatlonfrac[1] > 0: 
                ixbeg = (ilatlonfrac[1]-1)*fracx
                ixend = ilatlonfrac[1]*fracx+1
            else:
                ixbeg = 0
                ixend = fracx+1
    
#        ibeg = [iybeg, ixbeg]
#        iend = [iyend, ixend]
#        print 'find window; beginning', ibeg, 'end:', iend
        lon = ilon[iybeg:iyend,ixbeg:ixend]
        lat = ilat[iybeg:iyend,ixbeg:ixend]
#        print 'Looking within _______'
#        print 'lon:',lon[0,0],',',lon[-1,-1]
#        print 'lat:',lat[0,0],',',lat[-1,-1]
    
#        print 'lon _______'
#        print lon
#        print 'lat _______'
#        print lat
    
# Find point
        difflonlat = np.sqrt((lon-lonv)**2. + (lat-latv)**2.)
        if out:
            mindiffLl[iv] = np.min(difflonlat)
            ilatlon[:,iv] = index_mat(difflonlat, mindiffLl)
            ilatlon[0,iv] = ilatlon[0,iv] + iybeg
            ilatlon[1,iv] = ilatlon[1,iv] + ixbeg
        else:
            mindiffLl = np.min(difflonlat)
            ilatlon = index_mat(difflonlat, mindiffLl)
            ilatlon[0] = ilatlon[0] + iybeg
            ilatlon[1] = ilatlon[1] + ixbeg

#        print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
#        print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]
#        quit()

#        if mindiffLl == 0. and type(newvar[ilatlon[0],ilatlon[1]]) == type(amsk): 
#            percendone(iv,Ninpts,0.5,'done:')
        if mindiffLl == 0.: 
            iptpos = iptpos + 1
# Speeding it up!
#            if mindiffLl > mindiff:
#                print errormsg
#                print '  ' + fname + ': for point #', iv,'lon,lat in ' +             \
#                  'incomplet map:', lonvs[iv], ',', latvs[iv], 'there is ' +         \
#                  'not a set of lon,lat in the completed map closer than: ',         \
#                  mindiff, '!!'
#                print '    minimum difference:', mindiffLl
#                quit(-1)

# Speeding it up!
#            if ilatlon[0] >= 0 and ilatlon[1] >= 0:
#                newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
#                newvarin[ilatlon[0],ilatlon[1]] = iv
#                newvarinpt[iv] = 1
#                newvarindiff[iv] = mindiffLl
            newvar[ilatlon[0],ilatlon[1]] = ovar[iv]
            newvarin[ilatlon[0],ilatlon[1]] = iv
            newvarinpt[iv] = 1
            newvarindiff[iv] = mindiffLl

#                print 'Lluis iv:', newvarin[ilatlon[0],ilatlon[1]],                  \
#                  'localized:', newvarinpt[iv], 'values:',                           \
#                  newvar[ilatlon[0],ilatlon[1]], 'invalues:', ovar[iv],              \
#                  'mindist:', newvarindiff[iv], 'point:',ilatlon    
#            else:
#                print errormsg
#                print '  ' + fname + ': point iv:', iv, 'at', lonvs[iv], ',',        \
#                  latvs[iv],' not relocated !!'
#                print '    mindiffl:', mindiffLl, 'ilon:', ilatlon[1],               \
#                  'ilatlon:', ilatlon[1]
#                quit(-1)
            if np.mod(iptpos,fracs) == 0: 
                print 'Lluis syncronizing!'
                onc.sync()
    if out:
        return ilatlon, mindiffLl
    else:
        return

def lonlatFind(ilon, ilat, lonv, latv):
    """ Function to search a given value from a lon, lat 2D data
      ilon, ilat= original 2D matrices with the longitudes and the latitudes
      lonv, latv= longitude and latitude to find
      >>> npt=100
      >>> lonval0 = np.arange(0.,npt*1.,1.)
      >>> latval0 = np.arange(0.,npt*1.,1.)
      >>> lonval = np.zeros((npt,npt), dtype=np.float)
      >>> latval = np.zeros((npt,npt), dtype=np.float)

      >>> for i in range(npt):
      >>>     lonval[:,i] = lonval0[i]
      >>>     latval[i,:] = latval0[i]
      >>> lonlatFind(lonval, latval, 5.25, 5.25)
      (array([5, 5]), 0.35355339059327379)
    """
    fname = 'lonlatFind'

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

    nlon = np.min(ilon)
    xlon = np.max(ilon)
    nlat = np.min(ilat)
    xlat = np.max(ilat)

    if lonv < nlon or lonv > xlon:
        print errormsg
        print '  ' + fname + ': longitude outside data range!!'
        print '    given value:', lonv,'outside [',nlon,',',xlon,']'
        quit(-1) 
    if latv < nlat or latv > xlat:
        print errormsg
        print '  ' + fname + ': latitude outside data range!!'
        print '    given value:', latv,'outside [',nlat,',',xlat,']'
        quit(-1) 

# Find point
    difflonlat = np.sqrt((ilon-lonv)**2. + (ilat-latv)**2.)
    mindiffLl = np.min(difflonlat)
    ilatlon = index_mat(difflonlat, mindiffLl)

#    print 'mindiffLl:', mindiffLl, 'ilatlon:', ilatlon
#    print 'lon, lat:', lon[ilatlon[0],ilatlon[1]], ',', lat[ilatlon[0],ilatlon[1]]

    return ilatlon, mindiffLl
