Changeset 337 in lmdz_wrf for trunk/tools/validation_sim.py

Timestamp:

Feb 28, 2015, 4:22:10 PM (10 years ago)

Author:

lfita

Message:

Working version for a 3D trajectory

File:

• trunk/tools/validation_sim.py (modified) (17 diffs)

trunk/tools/validation_sim.py

@@ -1 +1 @@
 
 # L. Fita, LMD-Jussieu. February 2015
-## e.g. # 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. 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
+## 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
+## 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
+
 import numpy as np
 import os
@@ -21 +25 @@
 
 StringLength = 50
+
+# Number of grid points to take as 'environment' around the observed point
+Ngrid = 1
 
 def searchInlist(listname, nameFind):
@@ -162 +169 @@
       matB= matrix with the pother set of values
       val= couple of values to search
-    >>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
+    >>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
     [2 1 1]
     """
@@ -182 +189 @@
     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)
@@ -196 +209 @@
     mindist = np.min(dist)
    
-    matlist = list(dist.flatten())
-    ifound = matlist.index(mindist)
-
-    Ndims = len(matAshape)
-    valpos = np.zeros((Ndims), dtype=int)
+    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])
@@ -213 +227 @@
     return valpos
 
-def index_mat(mat,val):
+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
@@ -433 +485 @@
     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
+        '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 == '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
+        'WRFt': temperature 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'][:] + 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
+                shape = ncobj.variables['P'].shape
+
+
+            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 == '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 erromsg
+                print '  ' + fname + ": variable '" + varn + "' nor ready !!"
+                quit(-1)
+
+            print type(ncobj.variables)
+    
+            self.dimensions = dimensions
+            self.shape = shape
+            self.__values = varNOcheckv
+
+    def __getitem__(self,elem):
+        return self.__values[elem]
+
 
 ####### ###### ##### #### ### ## #
@@ -443 +683 @@
   "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
   "'trajectory': following a trajectory"
-simopers =['sumc','subc','mulc','divc']
-#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
+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', 'WRFt', 'WRFz']
+varNOcheckinf = "'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; " +    \
+ "'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 operations: " + opersinf +          \
+  "(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
+statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
 
 parser = OptionParser()
-parser.add_option("-d", "--dimensions", dest="dims",
-  help="[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from each source ([DIM]='X','Y','Z','T'; None, no value)",
-  metavar="VALUES")
+parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
 parser.add_option("-D", "--vardimensions", dest="vardims",
-  help="[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, no value)", metavar="VALUES")
+  help=vardimshelp, metavar="VALUES")
 parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
   help=strkObs, metavar="FILE")
@@ -464 +717 @@
 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="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to validate and if necessary operation and value", metavar="VALUES")
+  help=varshelp, metavar="VALUES")
 
 (opts, args) = parser.parse_args()
@@ -612 +868 @@
           dims[dn][0] + "' !!"
         quit(-1)
-    if not osim.variables.has_key(vardims[dn][0]):
+    if not osim.variables.has_key(vardims[dn][0]) and not                            \
+      searchInlist(varNOcheck,vardims[dn][0]):
         print errormsg
         print '  ' + main + ": simulation file does not have varibale dimension '" + \
           vardims[dn][0] + "' !!"
         quit(-1)
-    valdimsim[dn] = osim.variables[vardims[dn][0]][:]
+    if searchInlist(varNOcheck,vardims[dn][0]):
+        valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
+    else:
+        valdimsim[dn] = osim.variables[vardims[dn][0]][:]
 
 # General characteristics
@@ -625 +885 @@
 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 dimtobs:
+    for it in range(dimtobs):
         trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
           [valdimobs['X'][it],valdimobss['Y'][it]])
@@ -648 +910 @@
 
 elif obskind == 'trajectory':
-    if dims.has_key('Z'):
-        trajpos = np.zeros((3,dimt),dtype=int)
-        for it in dimtobs:
-            trajpos[0:1,it] = index_2mat(valdimsim['Y'],valdimsim['X'],              \
-              [valdimobs['Y'][it],valdimobss['X'][it]])
-            trajpos[2,it] = index_mat(valdimsim['Z'],valdimobs['Z'][it])
+    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:
-        trajpos = np.zeros((2,dimt),dtype=int)
-        for it in dimtobs:
-            trajpos[:,it] = index_2mat(valdimsim['Y'],valdimsim['X'],                \
-              [valdimobs['Y'][it],valdimobss['X'][it]])
+        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 and 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
@@ -670 +973 @@
 # Concident times
 ##
-coindtvalues = []
+coindtvalues0 = []
 for it in range(dimtsim):    
     ot = 0
@@ -678 +981 @@
             ot = ito
             tdist = simobstimes[it] - valdimobs['T'][ito]
-            coindtvalues.append([it,ito,simobstimes[it],valdimobs['T'][ito],tdist])
-
-Ncoindt = len(coindtvalues)
+            coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],      \
+              tdist])
+
+coindtvalues = np.array(coindtvalues0, dtype=np.float)
+
+Ncoindt = len(coindtvalues[:,0])
 print main + ': found',Ncoindt,'coincident times between simulation and observations'
+
+if Ncoindt == 0:
+    print warnmsg
+    print main + ': no coincident times found !!'
+    print '  stopping it'
+    quit(-1)
 
 # Validating
@@ -690 +1002 @@
 # Dimensions
 newdim = onewnc.createDimension('time',None)
+newdim = onewnc.createDimension('obstime',None)
 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('stats',5)
 
 # Variable dimensions
 ##
-newvar = onewnc.createVariable('simtime','f8',('time'))
-basicvardef(newvar, 'simtime', 'time simulation', obstunits )
-set_attribute(newvar, 'calendar', 'standard')
-newvar[:] = coindtvalues[:][2]
-
 newvar = onewnc.createVariable('obstime','f8',('time'))
 basicvardef(newvar, 'obstime', 'time observations', obstunits )
 set_attribute(newvar, 'calendar', 'standard')
-newvar[:] = coindtvalues[:][3]
+newvar[:] = coindtvalues[:,3]
 
 newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
@@ -709 +1020 @@
 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)
+
+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)
+    trjsim = np.zeros((4,Ncoindt), dtype=int)
+    trjsimval = np.zeros((4,Ncoindt), dtype=np.float)
+else:
+    newdim = onewnc.createDimension('simtrj',3)
+    trjsim = np.zeros((3,Ncoindt), dtype=int)
+    trjsimval = np.zeros((3,Ncoindt), dtype=np.float)
+
 Nvars = len(valvars)
 for ivar in range(Nvars):
     simobsvalues = []
+# Values spatially around the point (+/- [Ngrid] points)
+    simobsSvalues = []
 
     varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
@@ -720 +1056 @@
           "' !!"
         quit(-1)
+
     if not osim.variables.has_key(valvars[ivar][0]):
-        print errormsg
-        print '  ' + main + ": simulation file has not '" + valvars[ivar][0] + "' !!"
-        quit(-1)
+        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]]
 
     ovobs = oobs.variables[valvars[ivar][1]]
-    ovsim = osim.variables[valvars[ivar][0]]
+    if dims.has_key('Z'):
+        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1, Ngrid*2+1),         \
+          dtype = np.float)
+    else:
+        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1), dtype = np.float)
 
     if obskind == 'multi-points':
         for it in range(Ncoindt):
-            slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                     \
-              dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                            \
+            slicev = dims['X'][0] + ':' + str(trajpos[2,it]) + '|' +                 \
+              dims['Y'][0]+ ':' + str(trajpos[1,it]) + '|' +                         \
+              dims['T'][0]+ ':' + str(coindtvalues[it][0])
+            slicevar, dimslice = slice_variable(ovsim, slicev)
+            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
+            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
+              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
+              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
               dims['T'][0]+':'+str(coindtvalues[it][0])
             slicevar, dimslice = slice_variable(ovsim, slicev)
-            simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
+            simobsSvalues[it,:,:] = slicevar
+
     elif obskind == 'single-station':
         for it in range(Ncoindt):
-            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                       \
-              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                              \
+            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                      \
+              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                             \
               dims['T'][0]+':'+str(coindtvalues[it][0])
             slicevar, dimslice = slice_variable(ovsim, slicev)
-            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
+            simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
+            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
+              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
+              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
+              dims['T'][0] + ':' + str(coindtvalues[it,0])
+            slicevar, dimslice = slice_variable(ovsim, slicev)
+            simobsSvalues[it,:,:] = slicevar
     elif obskind == 'trajectory':
         if dims.has_key('Z'):
             for it in range(Ncoindt):
-                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
-                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
-                  dims['Z'][0]+':'+str(trajpos[0,it]) + '|' +                        \
-                  dims['T'][0]+':'+str(coindtvalues[it][0])
-                slicevar, dimslice = slice_variable(ovsim, slicev)
-                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
-                print simobsvalues[varsimobs][:][it]
+                if notfound[it] == 0:
+                    ito = int(coindtvalues[it,1])
+                    trajpos[2,ito] = index_mat(valdimsim['Z'][it,:,trajpos[1,ito],   \
+                     trajpos[0,ito]], valdimobs['Z'][ito])
+                    slicev = dims['X'][0]+':'+str(trajpos[0,ito]) + '|' +            \
+                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
+                      dims['Z'][0]+':'+str(trajpos[2,ito]) + '|' +                   \
+                      dims['T'][0]+':'+str(int(coindtvalues[it,0]))
+                    slicevar, dimslice = slice_variable(ovsim, slicev)
+                    simobsvalues.append([ slicevar, ovobs[int(ito)]])
+                    slicev = dims['X'][0] + ':' + str(trajpos[0,ito]-Ngrid) + '@' +  \
+                      str(trajpos[0,ito]+Ngrid+1) + '|' + dims['Y'][0] + ':' +       \
+                      str(trajpos[1,ito]-Ngrid) + '@' + str(trajpos[1,ito]+Ngrid+1)+ \
+                      '|' + dims['Z'][0] + ':' + str(trajpos[2,ito]-Ngrid) + '@' +   \
+                      str(trajpos[2,ito]+Ngrid+1) + '|' + dims['T'][0] + ':' +       \
+                      str(int(coindtvalues[it,0]))
+                    slicevar, dimslice = slice_variable(ovsim, slicev)
+                    simobsSvalues[it,:,:,:] = slicevar
+                    if ivar == 0:
+                        trjsim[0,it] = trajpos[0,ito]
+                        trjsim[1,it] = trajpos[1,ito]
+                        trjsim[2,it] = trajpos[2,ito]
+                        trjsim[3,it] = coindtvalues[it,0]
+                else:
+                    simobsvalues.append([fillValueF, fillValueF])
+                    simobsSvalues[it,:,:,:]= np.ones((Ngrid*2+1,Ngrid*2+1,Ngrid*2+1),\
+                      dtype = np.float)*fillValueF
         else:
             for it in range(Ncoindt):
-                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
-                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
-                  dims['T'][0]+':'+str(coindtvalues[it][0])
-                slicevar, dimslice = slice_variable(ovsim, slicev)
-                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
+                if notfound[it] == 0:
+                    ito = coindtvalues[it,1]
+                    slicev = dims['X'][0]+':'+str(trajpos[2,ito]) + '|' +            \
+                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
+                      dims['T'][0]+':'+str(coindtvalues[ito,0])
+                    slicevar, dimslice = slice_variable(ovsim, slicev)
+                    simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
+                    slicev = dims['X'][0] + ':' + str(trajpos[0,it]-Ngrid) + '@' +   \
+                      str(trajpos[0,it]+Ngrid) + '|' + dims['Y'][0] + ':' +          \
+                      str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) +    \
+                      '|' + dims['T'][0] + ':' + str(coindtvalues[it,0])
+                    slicevar, dimslice = slice_variable(ovsim, slicev)
+                    simobsSvalues[it,:,:] = slicevar
+                else:
+                    simobsvalues.append([fillValue, fillValue])
+                    simobsSvalues[it,:,:] = np.ones((Ngrid*2+1,Ngrid*2+1),           \
+                      dtype = np.float)*fillValueF
                 print simobsvalues[varsimobs][:][it]
 
-    newvar = onewnc.createVariable(varsimobs, 'f', ('time','couple'))
-    descvar = 'couples of simulated: ' + valvars[ivar][0] + ' and observed ' +        \
+# Statistics around values
+    aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
+    for it in range(Ncoindt):
+        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])
+
+# Values to netCDF
+    newvar = onewnc.createVariable(varsimobs, 'f', ('time','couple'),                \
+      fill_value=fillValueF)
+    descvar = 'couples of simulated: ' + valvars[ivar][0] + ' and observed ' +       \
       valvars[ivar][1]
     basicvardef(newvar, varsimobs, descvar, ovobs.getncattr('units'))
- 
     newvar[:] = np.array(simobsvalues)
 
+# Around values
+    if dims.has_key('Z'):
+        newdim = onewnc.createDimension('zaround',Ngrid*2+1)
+        newvar = onewnc.createVariable( valvars[ivar][0] + 'around', 'f',            \
+          ('time','zaround','yaround','xaround'), fill_value=fillValueF)
+    else:
+        newvar = onewnc.createVariable( valvars[ivar][0] + 'around', 'f',            \
+          ('time','yaround','xaround'), fill_value=fillValueF)
+
+    descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
+    basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
+    newvar[:] = simobsSvalues
+
+# Statistics
+    newvar = onewnc.createVariable(varsimobs + 'staround', 'f', ('time','stats'),   \
+      fill_value=fillValueF)
+    descvar = 'around simulated statisitcs: ' + valvars[ivar][0]
+    basicvardef(newvar, varsimobs + 'staround', descvar, ovobs.getncattr('units'))
+    newvar[:] = aroundstats.transpose()
+
     onewnc.sync()
+
+newvar = onewnc.createVariable('simtrj','i',('time','simtrj'))
+basicvardef(newvar, 'simtrj', 'coordinates [X,Y,Z,T] of the coincident trajectory ' +\
+  'in sim', obstunits)
+newvar[:] = trjsim.transpose()
 
 # Global attributes