# L. Fita, LMD-Jussieu. February 2015
## e.g. sfcEneAvigon # validation_sim.py -d X@west_east@None,Y@south_north@None,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,T@time@time -k single-station -l 4.878773,43.915876,12. -o /home/lluis/DATA/obs/HyMeX/IOP15/sfcEnergyBalance_Avignon/OBSnetcdf.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v HFX@H,LH@LE,GRDFLX@G
## e.g. AIREP # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@alti,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/AIREP/2012/10/AIREP_121018.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@t,WRFtd@td,WRFws@u,WRFwd@dd
## e.g. ATRCore # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/ATRCore/V3/ATR_1Hz-HYMEXBDD-SOP1-v3_20121018_as120051.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@air_temperature@subc@273.15,WRFp@air_pressure,WRFrh@relative_humidity,WRFrh@relative_humidity_Rosemount,WRFwd@wind_from_direction,WRFws@wind_speed
## e.g. BAMED # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/BAMED/BAMED_SOP1_B12_TOT5.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@tas_north,WRFp@pressure,WRFrh@hus,U@uas,V@vas

import numpy as np
import os
import re
from optparse import OptionParser
from netCDF4 import Dataset as NetCDFFile
from scipy import stats as sts
import numpy.ma as ma

main = 'validation_sim.py'
errormsg = 'ERROR -- errror -- ERROR -- error'
warnmsg = 'WARNING -- warning -- WARNING -- warning'

# version
version=1.0

# Filling values for floats, integer and string
fillValueF = 1.e20
fillValueI = -99999
fillValueS = '---'

StringLength = 50

# Number of grid points to take as 'environment' around the observed point
Ngrid = 1

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 set_attribute(ncvar, attrname, attrvalue):
    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
    ncvar = object netcdf variable
    attrname = name of the attribute
    attrvalue = value of the attribute
    """
    import numpy as np
    from netCDF4 import Dataset as NetCDFFile

    attvar = ncvar.ncattrs()
    if searchInlist(attvar, attrname):
        attr = ncvar.delncattr(attrname)

    attr = ncvar.setncattr(attrname, attrvalue)

    return ncvar

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

    return

def writing_str_nc(varo, values, Lchar):
    """ Function to write string values in a netCDF variable as a chain of 1char values
    varo= netCDF variable object
    values = list of values to introduce
    Lchar = length of the string in the netCDF file
    """

    Nvals = len(values)
    for iv in range(Nvals):    
        stringv=values[iv]  
        charvals = np.chararray(Lchar)
        Lstr = len(stringv)
        charvals[Lstr:Lchar] = ''

        for ich in range(Lstr):
            charvals[ich] = stringv[ich:ich+1]

        varo[iv,:] = charvals

    return

def index_3mat(matA,matB,matC,val):
    """ Function to provide the coordinates of a given value inside three matrix simultaneously
    index_mat(matA,matB,matC,val)
      matA= matrix with one set of values
      matB= matrix with the other set of values
      matB= matrix with the third set of values
      val= triplet of values to search
    >>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),np.arange(200,227).reshape(3,3,3),[22,122,222])
    [2 1 1]
    """
    fname = 'index_3mat'

    matAshape = matA.shape
    matBshape = matB.shape
    matCshape = matC.shape

    for idv in range(len(matAshape)):
        if matAshape[idv] != matBshape[idv]:
            print errormsg
            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
              'and B:',matBshape[idv],'does not coincide!!'
            quit(-1)
        if matAshape[idv] != matCshape[idv]:
            print errormsg
            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
              'and C:',matCshape[idv],'does not coincide!!'
            quit(-1)

    minA = np.min(matA)
    maxA = np.max(matA)
    minB = np.min(matB)
    maxB = np.max(matB)
    minC = np.min(matC)
    maxC = np.max(matC)

    if val[0] < minA or val[0] > maxA:
        print warnmsg
        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
          maxA,'!!'
    if val[1] < minB or val[1] > maxB:
        print warnmsg
        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
          maxB,'!!'
    if val[2] < minC or val[2] > maxC:
        print warnmsg
        print '  ' + fname + ': second value:',val[2],'outside matC range',minC,',',  \
          maxC,'!!'

    dist = np.zeros(tuple(matAshape), dtype=np.float)
    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2 +     \
      (matC - np.float(val[2]))**2)

    mindist = np.min(dist)
   
    matlist = list(dist.flatten())
    ifound = matlist.index(mindist)

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

    for dimid in range(Ndims):
        baseprevdims[dimid] = np.product(matAshape[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 index_2mat(matA,matB,val):
    """ Function to provide the coordinates of a given value inside two matrix simultaneously
    index_mat(matA,matB,val)
      matA= matrix with one set of values
      matB= matrix with the pother set of values
      val= couple of values to search
    >>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
    [2 1 1]
    """
    fname = 'index_2mat'

    matAshape = matA.shape
    matBshape = matB.shape

    for idv in range(len(matAshape)):
        if matAshape[idv] != matBshape[idv]:
            print errormsg
            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
              'and B:',matBshape[idv],'does not coincide!!'
            quit(-1)

    minA = np.min(matA)
    maxA = np.max(matA)
    minB = np.min(matB)
    maxB = np.max(matB)

    Ndims = len(matAshape)
#    valpos = np.ones((Ndims), dtype=int)*-1.
    valpos = np.zeros((Ndims), dtype=int)

    if val[0] < minA or val[0] > maxA:
        print warnmsg
        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
          maxA,'!!'
        return valpos
    if val[1] < minB or val[1] > maxB:
        print warnmsg
        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
          maxB,'!!'
        return valpos

    dist = np.zeros(tuple(matAshape), dtype=np.float)
    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2)

    mindist = np.min(dist)
   
    if mindist != mindist:
        print '  ' + fname + ': wrong minimal distance',mindist,'!!'
        return valpos
    else:
        matlist = list(dist.flatten())
        ifound = matlist.index(mindist)

    baseprevdims = np.zeros((Ndims), dtype=int)
    for dimid in range(Ndims):
        baseprevdims[dimid] = np.product(matAshape[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 index_mat(matA,val):
    """ Function to provide the coordinates of a given value inside a matrix
    index_mat(matA,val)
      matA= matrix with one set of values
      val= couple of values to search
    >>> index_mat(np.arange(27),22.3)
    22
    """
    fname = 'index_mat'

    matAshape = matA.shape

    minA = np.min(matA)
    maxA = np.max(matA)

    Ndims = len(matAshape)
#    valpos = np.ones((Ndims), dtype=int)*-1.
    valpos = np.zeros((Ndims), dtype=int)

    if val < minA or val > maxA:
        print warnmsg
        print '  ' + fname + ': first value:',val,'outside matA range',minA,',',     \
          maxA,'!!'
        return valpos

    dist = np.zeros(tuple(matAshape), dtype=np.float)
    dist = (matA - np.float(val))**2

    mindist = np.min(dist)
    if mindist != mindist:
        print '  ' + fname + ': wrong minimal distance',mindist,'!!'
        return valpos
   
    matlist = list(dist.flatten())
    valpos = matlist.index(mindist)

    return valpos

def index_mat_exact(mat,val):
    """ Function to provide the coordinates of a given exact 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 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 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 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 slice_variable(varobj, dimslice):
    """ Function to return a slice of a given variable according to values to its 
      dimensions
    slice_variable(varobj, dimslice)
      varobj= object wit the variable
      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension 
        [value]: 
          * [integer]: which value of the dimension
          * -1: all along the dimension
          * -9: last value of the dimension
          * [beg]@[end] slice from [beg] to [end]
    """
    fname = 'slice_variable'

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

    vardims = varobj.dimensions
    Ndimvar = len(vardims)

    Ndimcut = len(dimslice.split('|'))
    dimsl = dimslice.split('|')

    varvalsdim = []
    dimnslice = []

    for idd in range(Ndimvar):
        for idc in range(Ndimcut):
            dimcutn = dimsl[idc].split(':')[0]
            dimcutv = dimsl[idc].split(':')[1]
            if vardims[idd] == dimcutn: 
                posfrac = dimcutv.find('@')
                if posfrac != -1:
                    inifrac = int(dimcutv.split('@')[0])
                    endfrac = int(dimcutv.split('@')[1])
                    varvalsdim.append(slice(inifrac,endfrac))
                    dimnslice.append(vardims[idd])
                else:
                    if int(dimcutv) == -1:
                        varvalsdim.append(slice(0,varobj.shape[idd]))
                        dimnslice.append(vardims[idd])
                    elif int(dimcutv) == -9:
                        varvalsdim.append(int(varobj.shape[idd])-1)
                    else:
                        varvalsdim.append(int(dimcutv))
                break

    varvalues = varobj[tuple(varvalsdim)]

    return varvalues, dimnslice

def func_compute_varNOcheck(ncobj, varn):
    """ Function to compute variables which are not originary in the file
      ncobj= netCDF object file
      varn = variable to compute: 
        'WRFdens': air density from WRF variables
        'WRFght': geopotential height from WRF variables
        'WRFp': pressure from WRF variables
        'WRFrh': relative humidty fom WRF variables
        'WRFt': temperature from WRF variables
        'WRFwds': surface wind direction from WRF variables
        'WRFwss': surface wind speed from WRF variables
        'WRFz': height from WRF variables
    """
    fname = 'compute_varNOcheck'

    if varn == 'WRFdens':
#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
#          'DNW)/g ...'
        grav = 9.81

# Just we need in in absolute values: Size of the central grid cell
##    dxval = ncobj.getncattr('DX')
##    dyval = ncobj.getncattr('DY')
##    mapfac = ncobj.variables['MAPFAC_M'][:]
##    area = dxval*dyval*mapfac
        dimensions = ncobj.variables['MU'].dimensions

        mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
        dnw = ncobj.variables['DNW'][:]

        varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
          dtype=np.float)
        levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)

        for it in range(mu.shape[0]):
            for iz in range(dnw.shape[1]):
                levval.fill(np.abs(dnw[it,iz]))
                varNOcheck[it,iz,:,:] = levval
                varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav

    elif varn == 'WRFght':
#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
        varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
        dimensions = ncobj.variables['PH'].dimensions

    elif varn == 'WRFp':
#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
        varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
        dimensions = ncobj.variables['P'].dimensions

    elif varn == 'WRFrh':
#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
#         ' equation (T,P) ...'
        p0=100000.
        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
        tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
        qv = ncobj.variables['QVAPOR'][:]

        data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
        data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

        varNOcheckv = qv/data2
        dimensions = ncobj.variables['P'].dimensions

    elif varn == 'WRFt':
#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
        p0=100000.
        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]

        varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
        dimensions = ncobj.variables['T'].dimensions

    elif varn == 'WRFwds':
#        print '    ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
        varNOcheckv = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
        dimensions = ncobj.variables['V10'].dimensions

    elif varn == 'WRFwss':
#        print '    ' + main + ': computing surface wind speed from WRF as SQRT(U**2 + V**2) ...'
        varNOcheckv = np.sqrt(ncobj.variables['U10'][:]*ncobj.variables['U10'][:] +  \
          ncobj.variables['V10'][:]*ncobj.variables['V10'][:])
        dimensions = ncobj.variables['U10'].dimensions

    elif varn == 'WRFz':
        grav = 9.81
#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
        varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
        dimensions = ncobj.variables['PH'].dimensions

    else:
        print erromsg
        print '  ' + fname + ": variable '" + varn + "' nor ready !!"
        quit(-1)

    return varNOcheck

class compute_varNOcheck(object):
    """ Class to compute variables which are not originary in the file
      ncobj= netCDF object file
      varn = variable to compute: 
        'WRFdens': air density from WRF variables
        'WRFght': geopotential height from WRF variables
        'WRFp': pressure from WRF variables
        'WRFrh': relative humidty fom WRF variables
        'TSrhs': surface relative humidty fom TS variables
        'WRFrhs': surface relative humidty fom WRF variables
        'WRFT': CF-time from WRF variables
        'WRFt': temperature from WRF variables
        'WRFtd': dew-point temperature from WRF variables
        'WRFwd': wind direction from WRF variables
        'TSwds': surface wind direction from TS variables
        'WRFwds': surface wind direction from WRF variables
        'WRFws': wind speed from WRF variables
        'TSwss': surface wind speed from TS variables
        'WRFwss': surface wind speed from WRF variables
        'WRFz': height from WRF variables
    """
    fname = 'compute_varNOcheck'

    def __init__(self, ncobj, varn):

        if ncobj is None:
            self = None
            self.dimensions = None
            self.shape = None
            self.__values = None
        else:
            if varn == 'WRFdens':
#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
#          'DNW)/g ...'
                grav = 9.81

# Just we need in in absolute values: Size of the central grid cell
##    dxval = ncobj.getncattr('DX')
##    dyval = ncobj.getncattr('DY')
##    mapfac = ncobj.variables['MAPFAC_M'][:]
##    area = dxval*dyval*mapfac
                dimensions = ncobj.variables['MU'].dimensions
                shape = ncobj.variables['MU'].shape

                mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
                dnw = ncobj.variables['DNW'][:]

                varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
                  dtype=np.float)
                levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)

                for it in range(mu.shape[0]):
                    for iz in range(dnw.shape[1]):
                        levval.fill(np.abs(dnw[it,iz]))
                        varNOcheck[it,iz,:,:] = levval
                        varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav

            elif varn == 'WRFght':
#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
                varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
                dimensions = ncobj.variables['PH'].dimensions
                shape = ncobj.variables['PH'].shape

            elif varn == 'WRFp':
#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
                varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
                dimensions = ncobj.variables['P'].dimensions
                shape = ncobj.variables['P'].shape

            elif varn == 'WRFrh':
#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
#         ' equation (T,P) ...'
                p0=100000.
                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
                tk = (ncobj.variables['T'][:])*(p/p0)**(2./7.)
                qv = ncobj.variables['QVAPOR'][:]

                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

                varNOcheckv = qv/data2
                dimensions = ncobj.variables['P'].dimensions
                shape = ncobj.variables['P'].shape

            elif varn == 'TSrhs':
#        print '    ' + main + ": computing surface relative humidity from TSs as 'Tetens'" +\
#         ' equation (T,P) ...'
                p0=100000.
                p=ncobj.variables['psfc'][:]
                tk = (ncobj.variables['t'][:])*(p/p0)**(2./7.)
                qv = ncobj.variables['q'][:]

                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

                varNOcheckv = qv/data2
                dimensions = ncobj.variables['psfc'].dimensions
                shape = ncobj.variables['psfc'].shape

            elif varn == 'WRFrhs':
#        print '    ' + main + ": computing surface relative humidity from WRF as 'Tetens'" +\
#         ' equation (T,P) ...'
                p0=100000.
                p=ncobj.variables['PSFC'][:]
                tk = (ncobj.variables['T2'][:] + 300.)*(p/p0)**(2./7.)
                qv = ncobj.variables['Q2'][:]

                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

                varNOcheckv = qv/data2
                dimensions = ncobj.variables['PSFC'].dimensions
                shape = ncobj.variables['PSFC'].shape

            elif varn == 'WRFT':
# To compute CF-times from WRF kind
#
                import datetime as dt

                times = ncobj.variables['Times']
                dimt = times.shape[0]
                varNOcheckv = np.zeros((dimt), dtype=np.float64)
                self.unitsval = 'seconds since 1949-12-01 00:00:00'
                refdate = datetimeStr_datetime('1949-12-01_00:00:00')

                dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
                shape = tuple([dimt])

                for it in range(dimt):
                    datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime',    \
                      'YmdHMS')
                    dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
                    difft = dateval - refdate
                    varNOcheckv[it] = difft.days*3600.*24. + difft.seconds +        \
                          np.float(int(difft.microseconds/10.e6))

            elif varn == 'WRFt':
#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
                p0=100000.
                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]

                varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
                dimensions = ncobj.variables['T'].dimensions
                shape = ncobj.variables['P'].shape

            elif varn == 'WRFtd':
#        print '    ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
# tacking from: http://en.wikipedia.org/wiki/Dew_point
                p0=100000.
                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]

                temp = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)

                qv = ncobj.variables['QVAPOR'][:]

                tk = temp - 273.15
                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)

                rh = qv/data2
                
                pa = rh * data1/100.
                varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))

                dimensions = ncobj.variables['T'].dimensions
                shape = ncobj.variables['P'].shape

            elif varn == 'WRFwd':
#        print '    ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...' 
                uwind = ncobj.variables['U'][:]
                vwind = ncobj.variables['V'][:]
                dx = uwind.shape[3]
                dy = vwind.shape[2]

# de-staggering
                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])

                theta = np.arctan2(ua, va)
                theta = np.where(theta < 0., theta + 2.*np.pi, theta)
                varNOcheckv = 360.*theta/(2.*np.pi)

                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
                shape = ua.shape

            elif varn == 'WRFwds':
#        print '    ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
                theta = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
                theta = np.where(theta < 0., theta + 2.*np.pi, theta)
                
                varNOcheckv = 360.*theta/(2.*np.pi)
                dimensions = ncobj.variables['V10'].dimensions
                shape = ncobj.variables['V10'].shape

            elif varn == 'TSwds':
#        print '    ' + main + ': computing surface wind direction from TSs as ATAN2(v,u) ...'
                theta = np.arctan2(ncobj.variables['v'][:], ncobj.variables['u'][:])
                theta = np.where(theta < 0., theta + 2.*np.pi, theta)

                varNOcheckv = 360.*theta/(2.*np.pi)
                dimensions = ncobj.variables['v'].dimensions
                shape = ncobj.variables['v'].shape

            elif varn == 'WRFws':
#        print '    ' + main + ': computing wind speed from WRF as SQRT(U**2 + V**2) ...' 
                uwind = ncobj.variables['U'][:]
                vwind = ncobj.variables['V'][:]
                dx = uwind.shape[3]
                dy = vwind.shape[2]
                 
# de-staggering
                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])

                varNOcheckv = np.sqrt(ua*ua + va*va)
                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
                shape = ua.shape

            elif varn == 'TSwss':
#        print '    ' + main + ': computing surface wind speed from TSs as SQRT(u**2 + v**2) ...'
                varNOcheckv = np.sqrt(ncobj.variables['u'][:]*                       \
                  ncobj.variables['u'][:] + ncobj.variables['v'][:]*                 \
                  ncobj.variables['v'][:]) 
                dimensions = ncobj.variables['u'].dimensions
                shape = ncobj.variables['u'].shape

            elif varn == 'WRFwss':
#        print '    ' + main + ': computing surface wind speed from WRF as SQRT(U**2 + V**2) ...'
                varNOcheckv = np.sqrt(ncobj.variables['U10'][:]*                     \
                  ncobj.variables['U10'][:] + ncobj.variables['V10'][:]*             \
                  ncobj.variables['V10'][:]) 
                dimensions = ncobj.variables['U10'].dimensions
                shape = ncobj.variables['U10'].shape

            elif varn == 'WRFz':
                grav = 9.81
#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...' 
                varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
                dimensions = ncobj.variables['PH'].dimensions
                shape = ncobj.variables['PH'].shape

            else:
                print errormsg
                print '  ' + fname + ": variable '" + varn + "' nor ready !!"
                quit(-1)

            self.dimensions = dimensions
            self.shape = shape
            self.__values = varNOcheckv

    def __getitem__(self,elem):
        return self.__values[elem]

def adding_station_desc(onc,stdesc):
    """ Function to add a station description in a netCDF file
      onc= netCDF object
      stdesc= station description lon, lat, height
    """
    fname = 'adding_station_desc'

    newvar = onc.createVariable( 'station', 'c', ('StrLength'))
    newvar[0:len(stdesc[0])] = stdesc[0]

    newdim = onc.createDimension('nst',1)

    if onc.variables.has_key('lon'):
        print warnmsg
        print '  ' + fname + ": variable 'lon' already exist !!"
        print "    renaming it as 'lonst'"
        lonname = 'lonst'
    else:
        lonname = 'lon'

    newvar = onc.createVariable( lonname, 'f4', ('nst'))
    basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
    newvar[:] = stdesc[1]

    if onc.variables.has_key('lat'):
        print warnmsg
        print '  ' + fname + ": variable 'lat' already exist !!"
        print "    renaming it as 'latst'"
        latname = 'latst'
    else:
        latname = 'lat'

    newvar = onc.createVariable( latname, 'f4', ('nst'))
    basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
    newvar[:] = stdesc[2]

    if onc.variables.has_key('height'):
        print warnmsg
        print '  ' + fname + ": variable 'height' already exist !!"
        print "    renaming it as 'heightst'"
        heightname = 'heightst'
    else:
        heightname = 'height'

    newvar = onc.createVariable( heightname, 'f4', ('nst'))
    basicvardef(newvar, heightname, 'height above sea level', 'm' )
    newvar[:] = stdesc[3]

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


def getting_ValidationValues(okind, dt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng, kvals):
    """ Function to get the values to validate accroding to the type of observation
      okind= observational kind
      dt= initial number of values to retrieve
      ds= dictionary with the names of the dimensions (sim, obs)
      trjpos= positions of the multi-stations (t, Y, X) or trajectory ([Z], Y, X)
      ovs= object with the values of the simulation
      ovs= object with the values of the observations
      tvalues= temporal values (sim. time step, obs. time step, sim t value, obs t value, t diff)
      oFill= Fill Value for the observations
      Ng= number of grid points around the observation
      kvals= kind of values
        'instantaneous':  values are taken as instantaneous values
        'tbackwardSmean':  simulated values are taken as time averages from back to the point
        'tbackwardOmean':  observed values are taken as time averages from back to the point
    return:
      sovalues= simulated values at the observation point and time
      soSvalues= values Ngrid points around the simulated point
      soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
      trjs= trajectory on the simulation space
    """
    fname = 'getting_ValidationValues'

    sovalues = []

    if kvals == 'instantaneous':
        dtf = dt
    elif kvals == 'tbackwardSmean':
        print '  ' + fname + ':',kvals,'!!'
        uniqt = np.unique(tvalues[:,3])
        dtf = len(uniqt)
        print '    initially we got',dt,'values which will become',dtf
        postf = {}
        for it in range(dtf):
            if it == 0:
                postf[uniqt[it]] = [0,0]
            elif it == 1:
                posprev = postf[uniqt[it-1]][1]
                posit = list(tvalues[:,3]).index(uniqt[it])
                postf[uniqt[it]] = [posprev, posit+1]
            else:
                posprev = postf[uniqt[it-1]][1]
                posit = list(tvalues[:,3]).index(uniqt[it])
                postf[uniqt[it]] = [posprev+1, posit+1]
    elif kvals == 'tbackwardOmean':
        print '  ' + fname + ':',kvals,'!!'
        uniqt = np.unique(tvalues[:,2])
        dtf = len(uniqt)
        print '     initially we got',dt,'values which will become',dtf
        print '     ==== NOT READY === '
        quit(-1)
    else:
        print errormsg
        print '  ' + fname + ": kind of values '" + kvals + "' not ready!!"
        quit(-1)

# Simulated values spatially around
    if ds.has_key('Z'):
        soSvalues = np.zeros((dt, Ng*2+1, Ng*2+1, Ng*2+1), dtype = np.float)
        if okind == 'trajectory':
            trjs = np.zeros((4,dt), dtype=int)
        else:
            trjs = None
    else:
        soSvalues = np.zeros((dt, Ng*2+1, Ng*2+1), dtype = np.float)
        if okind == 'trajectory':
            trjs = np.zeros((3,dt), dtype=int)
        else:
            trjs = None

    if okind == 'single-station':
        soTtvalues = np.zeros((dt,2), dtype=np.float)
    else:
        None

    if okind == 'multi-points':
        for it in range(dt):
            slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '|' +  ds['Y'][0] +      \
              ':' + str(trjpos[1,it]) + '|' + ds['T'][0]+ ':' + str(tvalues[it][0])
            slicevar, dimslice = slice_variable(ovs, slicev)
            sovalues.append([ slicevar, ovo[tvalues[it][1]]])
            slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' +                 \
              str(trjpos[2,it]+Ng) + '|' + ds['Y'][0] + ':' +                        \
              str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '|' +              \
              ds['T'][0]+':'+str(tvalues[it][0])
            slicevar, dimslice = slice_variable(ovs, slicev)
            soSvalues[it,:,:] = slicevar

    elif okind == 'single-station':
        for it in range(dt):
            ito = int(tvalues[it,1])
            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
                slicev = ds['X'][0] + ':' + str(stationpos[1]) + '|' +               \
                  ds['Y'][0] + ':' + str(stationpos[0]) + '|' +                      \
                  ds['T'][0] + ':' + str(int(tvalues[it][0]))
            else:
                slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
            slicevar, dimslice = slice_variable(ovs, slicev)
            if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
                sovalues.append([ slicevar, fillValueF])
#            elif ovo[int(ito)] != ovo[int(ito)]: 
#                sovalues.append([ slicevar, fillValueF])
            else:
                sovalues.append([ slicevar, ovo[int(ito)]])
            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
                slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' +            \
                  str(stationpos[1]+Ng+1) + '|' + ds['Y'][0] + ':' +                 \
                  str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '|' +      \
                  ds['T'][0] + ':' + str(int(tvalues[it,0]))
            else:
                slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
            slicevar, dimslice = slice_variable(ovs, slicev)
            soSvalues[it,:,:] = slicevar

            if it == 0:
                itoi = 0
                itof = int(tvalues[it,1]) / 2
            elif it == dt-1:
                itoi = int( (ito + int(tvalues[it-1,1])) / 2)
                itof = int(tvalues[it,1])
            else:
                itod = int( (ito - int(tvalues[it-1,1])) / 2 ) 
                itoi = ito - itod
                itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
                itof = ito + itod

            soTtvalues[it,0] = valdimobs['T'][itoi]
            soTtvalues[it,1] = valdimobs['T'][itof]

    elif okind == 'trajectory':
        if ds.has_key('Z'):
            for it in range(dt):
                ito = int(tvalues[it,1])
                if notfound[ito] == 0:
                    trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:,        \
                      trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
                    slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '|' +               \
                      ds['Y'][0]+':'+str(trjpos[1,ito]) + '|' +                      \
                      ds['Z'][0]+':'+str(trjpos[2,ito]) + '|' +                      \
                      ds['T'][0]+':'+str(int(tvalues[it,0]))
                    slicevar, dimslice = slice_variable(ovs, slicev)
                    sovalues.append([ slicevar, ovo[int(ito)]])
                    minx = np.max([trjpos[0,ito]-Ng,0])
                    maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
                    miny = np.max([trjpos[1,ito]-Ng,0])
                    maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
                    minz = np.max([trjpos[2,ito]-Ng,0])
                    maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])

                    slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '|' +  \
                      ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '|' +         \
                      ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '|' +         \
                      ds['T'][0] + ':' + str(int(tvalues[it,0]))
                    slicevar, dimslice = slice_variable(ovs, slicev)

                    sliceS = []
                    sliceS.append(it)
                    sliceS.append(slice(0,maxz-minz))
                    sliceS.append(slice(0,maxy-miny))
                    sliceS.append(slice(0,maxx-minx))

                    soSvalues[tuple(sliceS)] = slicevar
                    if ivar == 0:
                        trjs[0,it] = trjpos[0,ito]
                        trjs[1,it] = trjpos[1,ito]
                        trjs[2,it] = trjpos[2,ito]
                        trjs[3,it] = tvalues[it,0]
                else:
                    sovalues.append([fillValueF, fillValueF])
                    soSvalues[it,:,:,:]= np.ones((Ng*2+1,Ng*2+1,Ng*2+1),             \
                      dtype = np.float)*fillValueF
# 2D trajectory
        else:
            for it in range(dt):
                if notfound[it] == 0:
                    ito = tvalues[it,1]
                    slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '|' +               \
                      ds['Y'][0]+':'+str(trjpos[1,ito]) + '|' +                      \
                      ds['T'][0]+':'+str(tvalues[ito,0])
                    slicevar, dimslice = slice_variable(ovs, slicev)
                    sovalues.append([ slicevar, ovo[tvalues[it,1]]])
                    slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' +         \
                      str(trjpos[0,it]+Ng) + '|' + ds['Y'][0] + ':' +                \
                      str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) +            \
                      '|' + ds['T'][0] + ':' + str(tvalues[it,0])
                    slicevar, dimslice = slice_variable(ovs, slicev)
                    soSvalues[it,:,:] = slicevar
                else:
                    sovalues.append([fillValue, fillValue])
                    soSvalues[it,:,:] = np.ones((Ng*2+1,Ng*2+1),                     \
                      dtype = np.float)*fillValueF
                print sovalues[varsimobs][:][it]
    else:
        print errormsg
        print '  ' + fname + ": observatino kind '" + okind + "' not ready!!"
        quit(-1)

# Re-arranging final values 
##
    if kvals == 'instantaneous':
        return sovalues, soSvalues, soTtvalues, trjs

    elif kvals == 'tbackwardSmean':
        fsovalues = []
        if ds.has_key('Z'):
            fsoSvalues = np.zeros((dtf, Ng*2+1, Ng*2+1, Ng*2+1), dtype = np.float)
            if okind == 'trajectory':
                ftrjs = np.zeros((4,dtf), dtype=int)
            else:
                ftrjs = None
        else:
            fsoSvalues = np.zeros((dtf, Ng*2+1, Ng*2+1), dtype = np.float)
            if okind == 'trajectory':
                ftrjs = np.zeros((3,dtf), dtype=int)
            else:
                ftrjs = None

        if okind == 'single-station':
            fsoTtvalues = np.ones((dtf,2), dtype=np.float)*fillValueF
        else:
            None

        for it in range(1,dtf):
            tv = uniqt[it]
            intv = postf[tv]

# Temporal statistics
            sovs = np.array(sovalues[intv[0]:intv[1]])[:,0]
            minv = np.min(sovs)
            maxv = np.max(sovs)
            meanv = np.mean(sovs)
            stdv = np.std(sovs)

            fsovalues.append([meanv, np.array(sovalues[intv[0]:intv[1]])[0,1], minv, \
              maxv, stdv])
            if ds.has_key('Z'):
                if okind == 'trajectory':
                    for ip in range(4):
                        ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
                for iz in range(2*Ng+1):
                    for iy in range(2*Ng+1):
                        for ix in range(2*Ng+1):
                            fsoSvalues[it,iz,iy,ix] = np.mean(soSvalues[intv[0]:     \
                              intv[1],iz,iy,ix])
            else:
                if okind == 'trajectory':
                    for ip in range(3):
                        ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
                for iy in range(2*Ng+1):
                    for ix in range(2*Ng+1):
                        fsoSvalues[it,iy,ix] = np.mean(soSvalues[intv[0]:intv[1],    \
                          iy,ix])
            fsoTtvalues[it,0] = soTtvalues[intv[0],0]
            fsoTtvalues[it,1] = soTtvalues[intv[1],0]

        return fsovalues, fsoSvalues, fsoTtvalues, ftrjs

    elif kvals == 'tbackwardOmean':
        print '  ' + fname + ':',kvals,'!!'
        uniqt = np.unique(tvalues[:,2])
        dtf = len(uniqt)
        print '     initially we got',dt,'values which will become',dtf

    return 


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

strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"

kindobs=['multi-points', 'single-station', 'trajectory']
strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
  "'trajectory': following a trajectory"
simopers = ['sumc','subc','mulc','divc']
opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
  'substract [constant] to variables values; mulc,[constant]: multipy by ' +         \
  '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
  'variables values'
varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'TSrhs', 'WRFrhs', 'WRFT',       \
  'WRFt', 'WRFtd', 'WRFwd', 'TSwds', 'WRFwds', 'WRFws', 'TSwss', 'WRFwss', 'WRFz'] 
varNOcheckinf = "'WRFdens': air density from WRF variables; " +                      \
  "'WRFght': geopotentiali height from WRF variables; " +                            \
  "'WRFp': pressure from WRF variables; " +                                          \
  "'WRFrh': relative humidty fom WRF variables; " +                                  \
  "'WRFrhs': surface relative humidity from WRF variables; " +                       \
  "'WRFT': CF-time from WRF variables; " +                                           \
  "'WRFt': temperature from WRF variables; " +                                       \
  "'WRFtd': dew-point temperature from WRF variables; " +                            \
  "'WRFwd': wind direction from WRF variables; " +                                   \
  "'WRFwds': surface wind direction from WRF variables; " +                          \
  "'WRFws': wind speed from WRF variables; " +                                       \
  "'WRFwss': surface wind speed from WRF variables; " +                              \
  "'WRFz': height from WRF variables"

dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
  "each source ([DIM]='X','Y','Z','T'; None, no value)"
vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " +    \
  "names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " +   \
  "no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " +   \
  "validate and if necessary operation and value (sim. values) available " +         \
  "operations: " + opersinf + " (WRFdiagnosted variables also available: " +         \
  varNOcheckinf + ")"
statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae',  \
  'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR',            \
  'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2',                \
  'standard deviation']

parser = OptionParser()
parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
parser.add_option("-D", "--vardimensions", dest="vardims",
  help=vardimshelp, metavar="VALUES")
parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
  help=strkObs, metavar="FILE")
parser.add_option("-l", "--stationLocation", dest="stloc",  
  help="name (| for spaces), longitude, latitude and height of the station (only for 'single-station')", 
  metavar="FILE")
parser.add_option("-o", "--observation", dest="fobs",
  help="observations file to validate", metavar="FILE")
parser.add_option("-s", "--simulation", dest="fsim",
  help="simulation file to validate", metavar="FILE")
parser.add_option("-t", "--trajectoryfile", dest="trajf",
  help="file with grid points of the trajectory in the simulation grid ('simtrj')", 
  metavar="FILE")
parser.add_option("-v", "--variables", dest="vars",
  help=varshelp, metavar="VALUES")

(opts, args) = parser.parse_args()

####### ###### ##### #### ### ## #
# Number of different statistics according to the temporal coincidence
#  0: Exact time
#  1: Simulation values closest to observed times
#  2: Simulation values between consecutive observed times
Nstsim = 3

stdescsim = ['E', 'C', 'B']
prestdescsim = ['exact', 'closest', 'between']
Lstdescsim = ['exact time', 'cloest time', 'between observational time-steps']

#######    #######
## MAIN
    #######

ofile='validation_sim.nc'

if opts.dims is None:
    print errormsg
    print '  ' + main + ': No list of dimensions are provided!!'
    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
      ' is needed'
    quit(-1)
else:
    simdims = {}
    obsdims = {}
    print main +': couple of dimensions _______'
    dims = {}
    ds = opts.dims.split(',')
    for d in ds:
        dsecs = d.split('@')
        if len(dsecs) != 3:
            print errormsg
            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
            print '    [DIM]@[dimnsim]@[dimnobs]'
            quit(-1)
        if dsecs[1] != 'None':
            dims[dsecs[0]] = [dsecs[1], dsecs[2]]
            simdims[dsecs[0]] = dsecs[1]
            obsdims[dsecs[0]] = dsecs[2]

            print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
        
if opts.vardims is None:
    print errormsg
    print '  ' + main + ': No list of variables with dimension values are provided!!'
    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
      '[vardimtsim]@[vardimtobs] is needed'
    quit(-1)
else:
    print main +': couple of variable dimensions _______'
    vardims = {}
    ds = opts.vardims.split(',')
    for d in ds:
        dsecs = d.split('@')
        if len(dsecs) != 3:
            print errormsg
            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
            print '    [DIM]@[vardimnsim]@[vardimnobs]'
            quit(-1)
        if dsecs[1] != 'None':
            vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
            print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]

if opts.obskind is None:
    print errormsg
    print '  ' + main + ': No kind of observations provided !!'
    quit(-1)
else:
    obskind = opts.obskind
    if obskind == 'single-station':
        if opts.stloc is None:
            print errormsg
            print '  ' + main + ': No station location provided !!'
            quit(-1)
        else:
            stationdesc = [opts.stloc.split(',')[0].replace('|',' '),                \
              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
              np.float(opts.stloc.split(',')[3])]

if opts.fobs is None:
    print errormsg
    print '  ' + main + ': No observations file is provided!!'
    quit(-1)
else:
    if not os.path.isfile(opts.fobs):
        print errormsg
        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
        quit(-1)

if opts.fsim is None:
    print errormsg
    print '  ' + main + ': No simulation file is provided!!'
    quit(-1)
else:
    if not os.path.isfile(opts.fsim):
        print errormsg
        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
        quit(-1)

if opts.vars is None:
    print errormsg
    print '  ' + main + ': No list of couples of variables is provided!!'
    print '    a ',' list of values [varsim]@[varobs],... is needed'
    quit(-1)
else:
    valvars = []
    vs = opts.vars.split(',')
    for v in vs:
        vsecs = v.split('@')
        if len(vsecs) < 2:
            print errormsg
            print '  ' + main + ': wrong number of values in:',vsecs,                \
              ' at least 2 are needed !!'
            print '    [varsim]@[varobs]@[[oper][val]]'
            quit(-1)
        if len(vsecs) > 2:
            if not searchInlist(simopers,vsecs[2]): 
                print errormsg
                print main + ": operation on simulation values '" + vsecs[2] +       \
                  "' not ready !!"
                quit(-1)

        valvars.append(vsecs)

# Openning observations trajectory
##
oobs = NetCDFFile(opts.fobs, 'r')

valdimobs = {}
for dn in dims:
    print dn,':',dims[dn]
    if dims[dn][1] != 'None':
        if not oobs.dimensions.has_key(dims[dn][1]):
            print errormsg
            print '  ' + main + ": observations file does not have dimension '" +    \
              dims[dn][1] + "' !!"
            quit(-1)
        if vardims[dn][1] != 'None':
            if not oobs.variables.has_key(vardims[dn][1]):
                print errormsg
                print '  ' + main + ": observations file does not have varibale " +  \
                  "dimension '" + vardims[dn][1] + "' !!"
                quit(-1)
            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
    else:
        if dn == 'X':
            valdimobs[dn] = stationdesc[1]
        elif dn == 'Y':
            valdimobs[dn] = stationdesc[2]
        elif dn == 'Z':
            valdimobs[dn] = stationdesc[3]

osim = NetCDFFile(opts.fsim, 'r')

valdimsim = {}
for dn in dims:
    if dims[dn][0] != 'None':
        if not osim.dimensions.has_key(dims[dn][0]):
            print errormsg
            print '  ' + main + ": simulation file '" + opts.fsim +                  \
              "' does not have dimension '" + dims[dn][0] + "' !!"
            print '    it has: ',osim.dimensions
            quit(-1)

        if not osim.variables.has_key(vardims[dn][0]) and                            \
          not searchInlist(varNOcheck,vardims[dn][0]):
            print errormsg
            print '  ' + main + ": simulation file '" + opts.fsim +                  \
              "' does not have varibale dimension '" + vardims[dn][0] + "' !!"
            print '    it has variables:',osim.variables
            quit(-1)
        if searchInlist(varNOcheck,vardims[dn][0]):
            valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
        else:
            valdimsim[dn] = osim.variables[vardims[dn][0]][:]

# General characteristics
dimtobs = valdimobs['T'].shape[0]
dimtsim = valdimsim['T'].shape[0]

print main +': observational time-steps:',dimtobs,'simulation:',dimtsim

notfound = np.zeros((dimtobs), dtype=int)

if obskind == 'multi-points':
    trajpos = np.zeros((2,dimt),dtype=int)
    for it in range(dimtobs):
        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
          [valdimobs['X'][it],valdimobss['Y'][it]])
elif obskind == 'single-station':
    trajpos = None
    stationpos = np.zeros((2), dtype=int)
    if valdimsim.has_key('X') and valdimsim.has_key('Y'):
        stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],         \
          valdimobs['X']])
        iid = 0
        for idn in osim.variables[vardims['X'][0]].dimensions:
            if idn == dims['X'][0]:
                stationpos[1] = stsimpos[iid]
            elif idn == dims['Y'][0]:
                stationpos[0] = stsimpos[iid]

            iid = iid + 1
        print main + ': station point in simulation:', stationpos
        print '    station position:',valdimobs['X'],',',valdimobs['Y']
        print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',          \
          valdimsim['Y'][tuple(stsimpos)]
    else:
        print main + ': validation with two time-series !!'

elif obskind == 'trajectory':
    if opts.trajf is not None:
        if not os.path.isfile(opts.fsim):
            print errormsg
            print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
            quit(-1)
        else:
            otrjf = NetCDFFile(opts.fsim, 'r')
            trajpos[0,:] = otrjf.variables['obssimtrj'][0]
            trajpos[1,:] = otrjf.variables['obssimtrj'][1]
            otrjf.close()
    else:
        if dims.has_key('Z'):
            trajpos = np.zeros((3,dimtobs),dtype=int)
            for it in range(dimtobs):
                if np.mod(it*100./dimtobs,10.) == 0.:
                    print '    trajectory done: ',it*100./dimtobs,'%'
                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
                  [valdimobs['Y'][it],valdimobs['X'][it]])
                stationpos = np.zeros((2), dtype=int)
                iid = 0
                for idn in osim.variables[vardims['X'][0]].dimensions:
                    if idn == dims['X'][0]:
                        stationpos[1] = stsimpos[iid]
                    elif idn == dims['Y'][0]:
                        stationpos[0] = stsimpos[iid]
                    iid = iid + 1
                if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
              
                trajpos[0,it] = stationpos[0]
                trajpos[1,it] = stationpos[1]
# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
#                trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0],         \
#                  stationpos[1]], valdimobs['Z'][it])
        else:
            trajpos = np.zeros((2,dimtobs),dtype=int)
            for it in range(dimtobs):
                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
                  [valdimobs['Y'][it],valdimobss['X'][it]])
                stationpos = np.zeros((2), dtype=int)
                iid = 0
                for idn in osim.variables[vardims['X'][0]].dimensions:
                    if idn == dims['X'][0]:
                        stationpos[1] = stsimpos[iid]
                    elif idn == dims['Y'][0]:
                        stationpos[0] = stsimpos[iid]
                    iid = iid + 1
                if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1

                trajpos[0,it] = stationspos[0]
                trajpos[1,it] = stationspos[1]

        print main + ': not found',np.sum(notfound),'points of the trajectory'

# Getting times
tobj = oobs.variables[vardims['T'][1]]
obstunits = tobj.getncattr('units')
if vardims['T'][0] == 'WRFT':
    tsim = valdimsim['T'][:]
    simtunits = 'seconds since 1949-12-01 00:00:00'
else:
    tsim = osim.variables[vardims['T'][0]][:]
    otsim = osim.variables[vardims['T'][0]]
    simtunits = otsim.getncattr('units')

simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)

#
## Looking for exact/near times
#

# Exact Coincident times
##
exacttvalues0 = []
for it in range(dimtsim):    
    ot = 0
    for ito in range(ot,dimtobs-1):
        if valdimobs['T'][ito] == simobstimes[it]:
            ot = ito
            exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])

exacttvalues = np.array(exacttvalues0, dtype=np.float)
Nexactt = len(exacttvalues[:,0])

print main + ': found',Nexactt,'Temporal exact values in simulation and observations'

# Sim Closest times
##
Nsimt = 0
closesttvalues0 = []
closesttvalues0st = []
tsiminit = 0
tsimend = 0

dtsim = simobstimes[1] - simobstimes[0]

for it in range(dimtsim):
    ot = 0
    for ito in range(ot,dimtobs-1):
        if np.abs(valdimobs['T'][ito] - simobstimes[it]) <= dtsim/2.:
            ot = ito
            tdist = simobstimes[it] - valdimobs['T'][ito]
            closesttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],     \
              tdist])
            Nsimt = Nsimt + 1
            if tsiminit == 0: tsiminit=simobstimes[it]
            tsimend = simobstimes[it]

closesttvalues = np.array(closesttvalues0, dtype=np.float)

Nclosest = len(closesttvalues[:,0])
print main + ': found',Nclosest,'Simulation time-values closest to observations'

if Nclosest == 0:
    print warnmsg
    print main + ': no cclosest times found !!'
    print '  stopping it'
    quit(-1)

# Sim Coincident times
##
Nsimt = 0
coindtvalues0 = []
coindtvalues0st = []
tsiminit = 0
tsimend = 0

for it in range(dimtsim):    
    ot = 0
    for ito in range(ot,dimtobs-1):
        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
          simobstimes[it]:
            ot = ito
            tdist = simobstimes[it] - valdimobs['T'][ito]
            coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],     \
              tdist])
            Nsimt = Nsimt + 1
            if tsiminit == 0: tsiminit=simobstimes[it]
            tsimend = simobstimes[it]
        elif simobstimes[it] > valdimobs['T'][ito+1]:
            coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])

coindtvalues = np.array(coindtvalues0, dtype=np.float)
coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)

Ncoindt = len(coindtvalues[:,0])
print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' +     \
  'observed times) coincident times between simulation and observations'

if Ncoindt == 0:
    print warnmsg
    print main + ': no coincident times found !!'
    print '  stopping it'
    quit(-1)

# Validating
##

onewnc = NetCDFFile(ofile, 'w')

# Dimensions
for kst in range(Nstsim):
    newdim = onewnc.createDimension(prestdescsim[kst] + 'time',None)
    if stdescsim[kst] != 'E':
        newdim = onewnc.createDimension(prestdescsim[kst] + 'obstime',None)

newdim = onewnc.createDimension('bnds',2)
newdim = onewnc.createDimension('couple',2)
newdim = onewnc.createDimension('StrLength',StringLength)
newdim = onewnc.createDimension('xaround',Ngrid*2+1)
newdim = onewnc.createDimension('yaround',Ngrid*2+1)
newdim = onewnc.createDimension('gstats',13)
newdim = onewnc.createDimension('stats',5)
newdim = onewnc.createDimension('tstats',6)
newdim = onewnc.createDimension('Nstsim', 3)

# Variable dimensions
##
newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
basicvardef(newvar, 'couple', 'couples of values', '-')
writing_str_nc(newvar, ['sim','obs'], StringLength)

newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
basicvardef(newvar, 'statistics', 'statitics from values', '-')
writing_str_nc(newvar, statsn, StringLength)

newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
writing_str_nc(newvar, gstatsn, StringLength)

newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
writing_str_nc(newvar, ostatsn, StringLength)

newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
writing_str_nc(newvar, Lstdescsim, StringLength)

if obskind == 'trajectory':
    if dims.has_key('Z'):
        newdim = onewnc.createDimension('trj',3)
    else:
        newdim = onewnc.createDimension('trj',2)

    newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
    basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
    newvar[:] = trajpos.transpose()

if dims.has_key('Z'):
    newdim = onewnc.createDimension('simtrj',4)
else:
    newdim = onewnc.createDimension('simtrj',3)

Nvars = len(valvars)
for ivar in range(Nvars):
    simobsvalues = []

    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
    print '  ' + varsimobs + '... .. .'

    if not oobs.variables.has_key(valvars[ivar][1]):
        print errormsg
        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
          "' !!"
        quit(-1)

    if not osim.variables.has_key(valvars[ivar][0]):
        if not searchInlist(varNOcheck, valvars[ivar][0]):
            print errormsg
            print '  ' + main + ": simulation file has not '" + valvars[ivar][0] +   \
              "' !!"
            quit(-1)
        else:
            ovsim = compute_varNOcheck(osim, valvars[ivar][0])
    else:
        ovsim = osim.variables[valvars[ivar][0]]

    for idn in ovsim.dimensions:
        if not searchInlist(simdims.values(),idn):
            print errormsg
            print '  ' + main + ": dimension '" + idn + "' of variable '" +          \
              valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
            quit(-1)

    ovobs = oobs.variables[valvars[ivar][1]]
    if searchInlist(ovobs.ncattrs(),'_FillValue'): 
        oFillValue = ovobs.getncattr('_FillValue')
    else:
        oFillValue = None

    for kst in range(Nstsim):
        simobsvalues = [] 
       
        timedn = prestdescsim[kst] + 'time'
        timeobsdn = prestdescsim[kst] + 'obstime'
        print '    ' + prestdescsim[kst] + ' ...'

        if stdescsim[kst] == 'E':
# Observed and simualted exact times
            simobsvalues, simobsSvalues, simobsTtvalues, trjsim =                    \
              getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim,       \
              ovobs, exacttvalues, oFillValue, Ngrid, 'instantaneous')

            if ivar == 0:
                vname = prestdescsim[kst] + 'time'
                newvar = onewnc.createVariable(vname,'f8', (timedn))
                basicvardef(newvar, vname, 'exact coincident time observations and '+\
                  'simulation', obstunits)
                set_attribute(newvar, 'calendar', 'standard')
                newvar[:] = exacttvalues[:,3]

            dimt = Nexactt

        elif stdescsim[kst] == 'C':
# Simualted closest to Observed times
            simobsvalues, simobsSvalues, simobsTtvalues, trjsim =                    \
              getting_ValidationValues(obskind, Nclosest, dims, trajpos, ovsim,      \
              ovobs, closesttvalues, oFillValue, Ngrid, 'instantaneous')
            dimt = Nclosest


            if ivar == 0:
                vname = prestdescsim[kst] + 'time'
                newvar = onewnc.createVariable(vname,'f8', (timedn))
                basicvardef(newvar, vname, 'time simulations closest to observed ' + \
                  'values', obstunits )
                set_attribute(newvar, 'calendar', 'standard')
                newvar[:] = closesttvalues[:,2]

                vname  = prestdescsim[kst] + 'obstime'
                newvar = onewnc.createVariable(vname,'f8', (vname))
                basicvardef(newvar, vname, 'closest time observations', obstunits)
                set_attribute(newvar, 'calendar', 'standard')
                newvar[:] = closesttvalues[:,3]

        elif stdescsim[kst] == 'B':
# Observed values temporally around coincident times
            simobsvalues, simobsSvalues, simobsTtvalues, trjsim =                    \
              getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,\
              coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
            dimt = simobsSvalues.shape[0]

            if ivar == 0:
                vname = prestdescsim[kst] + 'time'
                newvar = onewnc.createVariable(vname,'f8', (timedn),                 \
                  fill_value = fillValueF)
                basicvardef(newvar, vname, 'simulation time between observations',   \
                  obstunits)
                set_attribute(newvar, 'calendar', 'standard')
                set_attribute(newvar, 'bounds', 'time_bnds')
                newvar[:] = simobsTtvalues[:,1]

                vname = prestdescsim[kst] + 'obstime'
                newvar = onewnc.createVariable(vname,'f8', (vname))
                basicvardef(newvar, vname, 'observed between time', obstunits )
                set_attribute(newvar, 'calendar', 'standard')
                newvar[:] = np.unique(coindtvalues[:,3])

# Re-arranging values...
        arrayvals = np.array(simobsvalues)
        if len(valvars[ivar]) > 2:
            const=np.float(valvars[ivar][3])
            if valvars[ivar][2] == 'sumc':
                simobsSvalues = simobsSvalues + const
                arrayvals[:,0] = arrayvals[:,0] + const
            elif valvars[ivar][2] == 'subc':
                simobsSvalues = simobsSvalues - const
                arrayvals[:,0] = arrayvals[:,0] - const
            elif valvars[ivar][2] == 'mulc':
                simobsSvalues = simobsSvalues * const
                arrayvals[:,0] = arrayvals[:,0] * const
            elif valvars[ivar][2] == 'divc':
                simobsSvalues = simobsSvalues / const
                arrayvals[:,0] = arrayvals[:,0] / const
            else:
                print errormsg
                print '  ' + fname + ": operation '"+valvars[ivar][2]+"' not ready!!"
                quit(-1)

        if kst == 0:
            simstats = np.zeros((Nstsim,5), dtype=np.float)
            obsstats = np.zeros((Nstsim,5), dtype=np.float)
            simobsstats = np.zeros((Nstsim,13), dtype=np.float)

# statisics sim
        simstats[kst,0] = np.min(arrayvals[:,0])
        simstats[kst,1] = np.max(arrayvals[:,0])
        simstats[kst,2] = np.mean(arrayvals[:,0])
        simstats[kst,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
        simstats[kst,4] = np.sqrt(simstats[kst,3] - simstats[kst,2]*simstats[kst,2])

# statisics obs
# Masking 'nan'
        obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF,            \
          arrayvals[:,1])

        obsmask = ma.masked_equal(obsmask0, fillValueF)
        obsmask2 = obsmask*obsmask

        obsstats[kst,0] = obsmask.min()
        obsstats[kst,1] = obsmask.max()
        obsstats[kst,2] = obsmask.mean()
        obsstats[kst,3] = obsmask2.mean()
        obsstats[kst,4] = np.sqrt(obsstats[kst,3] - obsstats[kst,2]*obsstats[kst,2])

# Statistics sim-obs
        diffvals = np.zeros((dimt), dtype=np.float)

        diffvals = arrayvals[:,0] - obsmask

        diff2vals = diffvals*diffvals
        sumdiff = diffvals.sum()
        sumdiff2 = diff2vals.sum()

        simobsstats[kst,0] = simstats[kst,0] - obsstats[kst,0]
        simobsstats[kst,1] = np.mean(arrayvals[:,0]*obsmask)
        simobsstats[kst,2] = diffvals.min()
        simobsstats[kst,3] = diffvals.max()
        simobsstats[kst,4] = diffvals.mean()
        simobsstats[kst,5] = np.abs(diffvals).mean()
        simobsstats[kst,6] = np.sqrt(diff2vals.mean())
        simobsstats[kst,7], simobsstats[kst,8] = sts.mstats.pearsonr(arrayvals[:,0], \
          obsmask)
# From:
#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
#Willmott, C. J. (1984). 'On the evaluation of model performance in physical 
#  geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
#  443-460
#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R. 
#  Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
#  Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
#   Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
#   35(1), 233-241 
#
# Deviation of residuals (SDR)
        simobsstats[kst,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[kst,0])*   \
          (diffvals-obsstats[kst,0]))))
# Index of Efficiency (IOE)
        obsbias2series = (obsmask - obsstats[kst,0])*(obsmask - obsstats[kst,0])
        sumobsbias2series = obsbias2series.sum()

        simobsstats[kst,10] = 1. - sumdiff2/(sumobsbias2series)
# Index of Agreement (IOA)
        simbias2series = arrayvals[:,0] - obsstats[kst,0]
        obsbias2series = obsmask - obsstats[kst,0]

        abssimbias2series = np.abs(simbias2series)
        absobsbias2series = np.where(obsbias2series<0, -obsbias2series,              \
          obsbias2series)

        abssimobsbias2series = (abssimbias2series+absobsbias2series)*(               \
          abssimbias2series + absobsbias2series)

        simobsstats[kst,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
# Fractional Mean Bias (FMB)
        simobsstats[kst,12]=(simstats[kst,0]-obsstats[kst,0])/(0.5*(simstats[kst,0] +\
          obsstats[kst,0]))

# Statistics around sim values
        aroundstats = np.zeros((5,dimt), dtype=np.float)
        for it in range(dimt):
            aroundstats[0,it] = np.min(simobsSvalues[it,])
            aroundstats[1,it] = np.max(simobsSvalues[it,])
            aroundstats[2,it] = np.mean(simobsSvalues[it,])
            aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
            aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]*       \
              aroundstats[2,it])

# sim Values to netCDF
        newvar = onewnc.createVariable(valvars[ivar][0] + '_' + stdescsim[kst] +     \
          'sim', 'f', (timedn), fill_value=fillValueF)
        descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
        basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
        newvar[:] = arrayvals[:,0]

# obs Values to netCDF
        if stdescsim[kst] != 'E':
            newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
          'obs', 'f', (timeobsdn), fill_value=fillValueF)
        else:
            newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
          'obs', 'f', (timedn), fill_value=fillValueF)
       
        descvar = prestdescsim[kst] + ' time observed: ' + valvars[ivar][1]
        basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
        newvar[:] = arrayvals[:,1]

# Around values
        if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
            vname = prestdescsim[kst] + valvars[ivar][0] + 'around'
        else:
            vname = prestdescsim[kst] + valvars[ivar][0] + 'Around'

        if dims.has_key('Z'):
            if not onewnc.dimensions.has_key('zaround'):
                newdim = onewnc.createDimension('zaround',Ngrid*2+1)
                newvar = onewnc.createVariable(vname, 'f', (timedn,'zaround',        \
                  'yaround','xaround'), fill_value=fillValueF)
        else:
            newvar = onewnc.createVariable(vname, 'f', (timedn,'yaround','xaround'), \
              fill_value=fillValueF)

        descvar = prestdescsim[kst] + 'around simulated values +/- grid values: ' +  \
          valvars[ivar][0]
        basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
        newvar[:] = simobsSvalues


# around sim Statistics
        if not searchInlist(onewnc.variables,prestdescsim[kst] + valvars[ivar][0] +  \
          'staround'):
            vname = prestdescsim[kst] + valvars[ivar][0] + 'staround'
        else:
            vname = prestdescsim[kst] + valvars[ivar][0] + 'Staround'

        newvar = onewnc.createVariable(vname, 'f', (timedn,'stats'),                 \
          fill_value=fillValueF)
        descvar = prestdescsim[kst] + ' around (' +  str(Ngrid) + ', ' + str(Ngrid) +\
          ') simulated statisitcs: ' + valvars[ivar][0]
        basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
        set_attribute(newvar, 'cell_methods', 'time_bnds')
        newvar[:] = aroundstats.transpose()

        if stdescsim[kst] == 'B':
            if not searchInlist(onewnc.variables, 'time_bnds'):
                newvar = onewnc.createVariable('time_bnds','f8',(timedn,'bnds'))
                basicvardef(newvar, 'time_bnds', timedn, obstunits )
                set_attribute(newvar, 'calendar', 'standard')
                newvar[:] = simobsTtvalues

# sim Statistics
    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
        vname = valvars[ivar][0] + 'stsim'
    else:
        vname = valvars[ivar][0] + 'stSim'

    newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'),                  \
      fill_value=fillValueF)
    descvar = 'simulated statisitcs: ' + valvars[ivar][0]
    basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
    newvar[:] = simstats

# obs Statistics
    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
        vname = valvars[ivar][1] + 'stobs'
    else:
        vname = valvars[ivar][1] + 'stObs'

    newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'),                  \
      fill_value=fillValueF)
    descvar = 'observed statisitcs: ' + valvars[ivar][1]
    basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
    newvar[:] = obsstats

# sim-obs Statistics
    if not searchInlist(onewnc.variables,varsimobs + 'st'):
        vname = varsimobs + 'st'
    else:
        vname = varSimobs + 'st'

    newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'gstats'),                 \
      fill_value=fillValueF)
    descvar = 'simulated-observed tatisitcs: ' + varsimobs
    basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
    newvar[:] = simobsstats

    onewnc.sync()

if trjsim is not None:
    newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
    basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' +         \
      'trajectory in sim', obstunits)
    newvar[:] = trjsim.transpose()

# Adding three variables with the station location, longitude, latitude and height
if obskind == 'single-station':
    adding_station_desc(onewnc,stationdesc)

# Global attributes
##
set_attribute(onewnc,'author_nc','Lluis Fita')
set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
  'LMD-Jussieu, UPMC, Paris')
set_attribute(onewnc,'country_nc','France')
set_attribute(onewnc,'script_nc',main)
set_attribute(onewnc,'version_script',version)
set_attribute(onewnc,'information',                                                 \
  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
set_attribute(onewnc,'simfile',opts.fsim)
set_attribute(onewnc,'obsfile',opts.fobs)

onewnc.sync()
onewnc.close()

print main + ": successfull writting of '" + ofile + "' !!"