Index: /trunk/tools/get_stations.py
===================================================================
--- /trunk/tools/get_stations.py	(revision 1679)
+++ /trunk/tools/get_stations.py	(revision 1679)
@@ -0,0 +1,733 @@
+## Getting statistics from a series of surface stations and soundings from any given netCDF file
+# L. Fita, CIMA, November 2017
+#
+import numpy as np
+import matplotlib as mpl
+mpl.use('Agg')
+from matplotlib.pylab import *
+import matplotlib.pyplot as plt
+from mpl_toolkits.basemap import Basemap
+import os
+from netCDF4 import Dataset as NetCDFFile
+import nc_var_tools as ncvar
+import generic_tools as gen
+import drawing_tools as drw
+import diag_tools as diag
+
+#sfcstatfile='ubicacion_datos_sup.csv'
+sfcstatfile='ubicacion_sfc_test.csv'
+
+# Surface stations columns for data
+sfclatcol=0
+sfcloncol=1
+sfclabcol=3
+sfchgtcol=2
+
+sndstatfile='ubicacion_sondeos.csv'
+
+# Sounding stations columns for data
+sndlatcol=1
+sndloncol=2
+sndlabcol=0
+sndhgtcol=3
+
+# comment char
+comchar = '#'
+
+# separation char
+sepchar = ','
+
+# missing value
+missvalS = ['0', 'NaN']
+
+# simulation file
+simfilen = '/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00'
+
+# axes dimensions
+axesdims = {'X': ['west_east', 'west_east_stag'],                                    \
+  'Y': ['south_north', 'south_north_stag'], 'Z': ['bottom_top', 'bottom_top_stag'],  \
+  'T': ['Time']}
+
+# axes variables
+axesvars = {'X': ['XLONG', 'XLONG_U'], 'Y': ['XLAT', 'XLAT_V'], 'Z': ['ZNU', 'ZNW'], \
+  'T': ['WRFtime']}
+
+# dim Variables
+dimvariables = {'west_east': 'XLONG', 'south_north': 'XLAT', 'bottom_top': 'ZNU',    \
+  'Time': 'WRFtime', 'west_east_stag': 'XLONG_U', 'south_north_stag': 'XLAT_U',      \
+  'bottom_top_stag': 'ZNW'}
+
+# SFC Variables
+sfcvariables = {'tas':'T2', 'uas':'WRFuas', 'vas':'WRFvas', 'huss':'Q2', 'ps':'PSFC',\
+  'zhgt':'HGT'}
+sfcrefvar = 'tas'
+
+# Vertical Variables
+sndvariables= {'ta':'WRFta', 'ua':'WRFua', 'va':'WRFva', 'hus':'QVAPOR',             \
+  'zg':'WRFzg', 'plev':'WRFp', 'tda': 'WRFtd'}
+sndrefvar = 'ta'
+
+# Variable pressure
+UnitsPress = 'Pa'
+Pressdimref = 'bottom_top'
+
+# Does pressure variable change with time?
+presstime = True
+
+# Map range
+#nlon = -80.
+#xlon = -50.
+#nlat = -60.
+#xlat = -20.
+
+# Map range RELAMPAGO
+nlon = -70.
+xlon = -57.
+nlat = -37.
+xlat = -27.
+
+# CF Temporal values
+ReferenceDate = '19491201000000'
+UnitsTime = 'minutes'
+
+####### ###### ##### #### ### ## #
+
+# Variables not to check their existence inside file
+NONcheckingvars = ['WRFp', 'WRFta', 'WRFtd', 'WRFtime', 'WRFua', 'WRFuas', 'WRFva',  \
+  'WRFvas', 'WRFzg']
+
+def creation_sfcstation_file(filen, lv, Lv, hv, lab, tunits, sfcvars, dimt, ifilen):
+    """ Function to create the structure of the surface station file
+      filen: name of the file
+      lv: value of the longitude of the station
+      Lv: value of the latitude of the station
+      hv: value of the height of the station
+      lab: label of the station
+      tunits: uints of time
+      sfcvars: variables to be included
+      dimt: quantity of time-steps
+      ifilen: name of the file from which data is retrieved
+    """
+    fname = 'creation_sfcstation_file'
+
+    Lstring = 256
+
+    onewnc = NetCDFFile(filen, 'w')
+    # Dimensions
+    newdim = onewnc.createDimension('time', None)
+    newdim = onewnc.createDimension('Lstring', Lstring)
+
+    # Variable-dimensions
+    newvar = onewnc.createVariable('time', 'f8', ('time'))
+    ncvar.basicvardef(newvar, 'time', 'Time', tunits)
+    ncvar.set_attribute(newvar, 'calendar', 'standard')
+    newvar.setncattr('axis', 'T')
+    newvar.setncattr('_CoordinateAxisType', 'Time')
+    onewnc.sync()
+
+    # station Variables
+    Llab = len(lab)
+    newvar = onewnc.createVariable('station', 'c', ('Lstring'))
+    newvar[0:Llab] = lab[:]
+    ncvar.basicvardef(newvar, 'station', 'Name of the station', '-')
+
+    newvar = onewnc.createVariable('lon', 'f8')
+    newvar[:] = lv
+    ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_west')
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+    
+    newvar = onewnc.createVariable('lat', 'f8')
+    newvar[:] = Lv
+    ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')
+    
+    newvar = onewnc.createVariable('height', 'f8')
+    newvar[:] = hv
+    ncvar.basicvardef(newvar, 'height', 'Height', 'm')
+    newvar.setncattr('axis', 'Z')
+    newvar.setncattr('_CoordinateAxisType', 'Height')
+    onewnc.sync()
+
+    newvar = onewnc.createVariable('flon', 'f8')
+    newvar[:] = lv
+    ncvar.basicvardef(newvar, 'file_lon', 'Longitude closest grid-point from file',  \
+      'degrees_west')
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+    
+    newvar = onewnc.createVariable('flat', 'f8')
+    newvar[:] = Lv
+    ncvar.basicvardef(newvar, 'file_lat', 'Latitude closest grid-point from file',   \
+      'degrees_north')
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')
+    
+    newvar = onewnc.createVariable('fheight', 'f8')
+    newvar[:] = hv
+    ncvar.basicvardef(newvar, 'file_height', 'Height closest grid-point from file',  \
+      'm')
+    newvar.setncattr('axis', 'Z')
+    newvar.setncattr('_CoordinateAxisType', 'Height')
+
+    newvar = onewnc.createVariable('ipoint', 'i')
+    newvar[:] = 0
+    ncvar.basicvardef(newvar, 'file_i', 'x-axis closest grid-point from file', '-')
+
+    newvar = onewnc.createVariable('jpoint', 'i')
+    newvar[:] = 0
+    ncvar.basicvardef(newvar, 'file_j', 'y-axis closest grid-point from file', '-')
+
+    onewnc.sync()
+
+    # Variables
+    NOvarvals = np.ones((dimt), dtype=np.float)*gen.fillValueF
+    for vn in sfcvars:
+        CFvalues = gen.variables_values(vn)
+        newvar = onewnc.createVariable(vn, 'f4', ('time'), fill_value=gen.fillValueF)
+        newvar[:] = NOvarvals[:]
+        ncvar.basicvardef(newvar, CFvalues[0], CFvalues[4].replace('|',' '),         \
+          CFvalues[5])
+
+    # Global attributes
+    ncvar.add_global_PyNCplot(onewnc, 'get_stations.py', fname, '0.1')
+    ncvar.set_attribute(newvar, 'data_origin', 'SMN')
+
+    onewnc.sync()
+    onewnc.close()
+
+    return
+
+def creation_sndstation_file(filen, lv, Lv, hv, lab, tunits, punits, ptime, sndvars, \
+  dimt, dimz, ifilen):
+    """ Function to create the structure of the surface station file
+      filen: name of the file
+      lv: value of the longitude of the station
+      Lv: value of the latitude of the station
+      hv: value of the height of the station
+      lab: label of the station
+      tunits: uints of time
+      punits: uints of pressure
+      ptime: does pressure evolves in time? (True/False)
+      sndvars: variables to be included
+      dimt: quantity of time-steps
+      dimz: quantity of vertical levels
+      ifilen: name of the file from which data is retrieved
+    """
+    fname = 'creation_sndstation_file'
+
+    Lstring = 256
+
+    onewnc = NetCDFFile(filen, 'w')
+    # Dimensions
+    newdim = onewnc.createDimension('time', None)
+    newdim = onewnc.createDimension('plev', dimz)
+    newdim = onewnc.createDimension('Lstring', Lstring)
+
+    # Variable-dimensions
+    newvar = onewnc.createVariable('time', 'f8', ('time'))
+    ncvar.basicvardef(newvar, 'time', 'Time', tunits)
+    ncvar.set_attribute(newvar, 'calendar', 'standard')
+    newvar.setncattr('axis', 'T')
+    newvar.setncattr('_CoordinateAxisType', 'Time')
+
+    if ptime:
+        newvar = onewnc.createVariable('plev', 'f8', ('time','plev'))
+        ncvar.basicvardef(newvar, 'plev', 'air pressure', punits)
+        ncvar.set_attribute(newvar, 'down', 'up')
+        newvar.setncattr('axis', 'Z')
+        newvar.setncattr('_CoordinateAxisType', 'pressure')
+    else:
+        newvar = onewnc.createVariable('plev', 'f8', ('plev'))
+        ncvar.basicvardef(newvar, 'plev', 'air pressure', punits)
+        ncvar.set_attribute(newvar, 'down', 'up')
+        newvar.setncattr('axis', 'Z')
+        newvar.setncattr('_CoordinateAxisType', 'pressure')
+    onewnc.sync()
+
+    # station Variables
+    Llab = len(lab)
+    newvar = onewnc.createVariable('station', 'c', ('Lstring'))
+    newvar[0:Llab] = lab[:]
+    ncvar.basicvardef(newvar, 'station', 'Name of the station', '-')
+
+    newvar = onewnc.createVariable('lon', 'f8')
+    newvar[:] = lv
+    ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_west')
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+    
+    newvar = onewnc.createVariable('lat', 'f8')
+    newvar[:] = Lv
+    ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')
+    
+    newvar = onewnc.createVariable('height', 'f8')
+    newvar[:] = hv
+    ncvar.basicvardef(newvar, 'height', 'Height', 'm')
+    newvar.setncattr('axis', 'Z')
+    newvar.setncattr('_CoordinateAxisType', 'Height')
+    onewnc.sync()
+
+    newvar = onewnc.createVariable('flon', 'f8')
+    newvar[:] = lv
+    ncvar.basicvardef(newvar, 'file_lon', 'Longitude closest grid-point from file',  \
+      'degrees_west')
+    newvar.setncattr('axis', 'X')
+    newvar.setncattr('_CoordinateAxisType', 'Lon')
+    
+    newvar = onewnc.createVariable('flat', 'f8')
+    newvar[:] = Lv
+    ncvar.basicvardef(newvar, 'file_lat', 'Latitude closest grid-point from file',   \
+      'degrees_north')
+    newvar.setncattr('axis', 'Y')
+    newvar.setncattr('_CoordinateAxisType', 'Lat')
+    
+    newvar = onewnc.createVariable('fheight', 'f8')
+    newvar[:] = hv
+    ncvar.basicvardef(newvar, 'file_height', 'Height closest grid-point from file',  \
+      'm')
+    newvar.setncattr('axis', 'Z')
+    newvar.setncattr('_CoordinateAxisType', 'Height')
+
+    newvar = onewnc.createVariable('ipoint', 'i')
+    newvar[:] = 0
+    ncvar.basicvardef(newvar, 'file_i', 'x-axis closest grid-point from file', '-')
+
+    newvar = onewnc.createVariable('jpoint', 'i')
+    newvar[:] = 0
+    ncvar.basicvardef(newvar, 'file_j', 'y-axis closest grid-point from file', '-')
+
+    onewnc.sync()
+
+    # Variables
+    NOvarvals = np.ones((dimt,dimz), dtype=np.float)*gen.fillValueF
+    for vn in sndvars:
+        if not onewnc.variables.has_key(vn):
+            CFvalues = gen.variables_values(vn)
+            newvar = onewnc.createVariable(vn, 'f4', ('time', 'plev'),               \
+              fill_value=gen.fillValueF)
+            newvar[:] = NOvarvals[:]
+            ncvar.basicvardef(newvar, CFvalues[0], CFvalues[4].replace('|',' '),     \
+              CFvalues[5])
+
+    # Global attributes
+    ncvar.add_global_PyNCplot(onewnc, 'get_stations.py', fname, '0.1')
+    ncvar.set_attribute(newvar, 'data_origin', 'SMN')
+
+    onewnc.sync()
+    onewnc.close()
+
+    return
+
+
+#######    #######
+## MAIN
+    #######
+
+yrref = ReferenceDate[0:4]
+monref = ReferenceDate[4:6]
+dayref = ReferenceDate[6:8]
+horref = ReferenceDate[8:10]
+minref = ReferenceDate[10:12]
+secref = ReferenceDate[12:14]
+
+utime = UnitsTime + ' since ' + yrref + '-' + monref + '-' + dayref + ' ' + horref + \
+  ':' + minref + ':' + secref
+
+# surface stations 
+osfcstatf = open(sfcstatfile, 'r')
+
+sfclonvals = []
+sfclatvals = []
+sfclabvals = []
+sfchgtvals = []
+
+for line in osfcstatf:
+    if line[0:1] != comchar and len(line) > 1:
+        vals = line.replace('\n','').split(sepchar)
+        if not gen.searchInlist(missvalS, vals[sfcloncol]) and                       \
+          not gen.searchInlist(missvalS, vals[sfclatcol]):
+            stlon = np.float(vals[sfcloncol])
+            stlat = np.float(vals[sfclatcol])
+            if stlon >= nlon and stlon <= xlon and stlat >= nlat and stlat <= xlat:
+                sfclonvals.append(stlon)
+                sfclatvals.append(stlat)
+                sfclabvals.append(vals[sfclabcol])
+                if gen.searchInlist(missvalS, vals[sfchgtcol]):
+                    if vals[sfchgtcol] == '0': 
+                        sfchgtvals.append(np.float(vals[sfchgtcol]))
+                    else: 
+                        sfchgtvals.append(gen.fillValueF)
+                else:
+                    sfchgtvals.append(np.float(vals[sfchgtcol]))
+
+osfcstatf.close()
+
+Nsfcstats = len(sfclonvals)
+print '  Number of surface stations: ', Nsfcstats
+
+# sounding stations 
+osndstatf = open(sndstatfile, 'r')
+
+sndlonvals = []
+sndlatvals = []
+sndlabvals = []
+sndhgtvals = []
+
+for line in osndstatf:
+    if line[0:1] != comchar and len(line) > 1:
+        vals = line.replace('\n','').split(sepchar)
+        if vals[sndloncol] != missvalS and vals[sndlatcol] != missvalS:
+            stlon = np.float(vals[sndloncol])
+            stlat = np.float(vals[sndlatcol])
+            if stlon >= nlon and stlon <= xlon and stlat >= nlat and stlat <= xlat:
+                sndlonvals.append(stlon)
+                sndlatvals.append(stlat)
+                sndlabvals.append(vals[sndlabcol])
+                sndhgtvals.append(np.float(vals[sndhgtcol]))
+
+osndstatf.close()
+Nsndstats = len(sndlonvals)
+print '  Number of sounding stations: ', Nsndstats
+
+if not os.path.isfile(simfilen):
+    print gen.errormsg
+    print "  file '" + simfilen + "' does not exist !!"
+    quit(-1)
+
+onc = NetCDFFile(simfilen, 'r')
+olistv = onc.variables.keys()
+olistv.sort()
+
+# Getting basic values for each dimension
+dimvarvalues = {}
+for dimn in dimvariables.keys():
+    varn = dimvariables[dimn]
+    if not gen.searchInlist(NONcheckingvars, varn) and not gen.searchInlist(olistv, varn):
+        print gen.errormsg
+        print "  file '" + simfilen + "' does not have variable '" + varn + "' !!"
+        print '    available ones: ', olistv
+        quit(-1)
+
+    # Except for temporal variables, we don't want dimensions with time-axis 
+    #  (assuming fixed)
+    if varn == 'WRFtime':
+        varvalues, CFtu= ncvar.compute_WRFtime(onc.variables['Times'], ReferenceDate,\
+          UnitsTime)
+        dt = varvalues.shape[0]
+        dimvarvalues[varn] = varvalues
+    else:
+        ovar = onc.variables[varn]
+        slicevar = {}
+        if not gen.searchInlist(axesvars['T'], varn):
+            for dn in ovar.dimensions:
+                if not gen.searchInlist(axesdims['T'], dn):
+                    slicevar[dn] = -1
+                else:
+                    slicevar[dn] = 0
+        else:
+            slicevar[dn] = -1
+            dt = len(onc.dimensions[dn])
+        slicevar, vardims = ncvar.SliceVarDict(ovar, slicevar)
+        dimvarvalues[varn] = ovar[tuple(slicevar)]
+
+    if dimn == Pressdimref: dz = len(onc.dimensions[dimn])
+
+# Retrieving surface data
+##
+
+# Had already uas & vas been computed?
+computeduasvas = False
+
+# Diagnosted variables
+diagvars = {}
+
+# Coinident surface station
+sfccoinc = {}
+
+for ist in range(Nsfcstats):
+    jumpstation = False
+
+    lonv = sfclonvals[ist]
+    latv = sfclatvals[ist]
+    labelv = sfclabvals[ist]
+    heightv = sfchgtvals[ist]
+
+    stfilen = 'sfc_station_' + labelv + '.nc'
+
+    creation_sfcstation_file(stfilen, lonv, latv, heightv, labelv, utime,            \
+      sfcvariables.keys(), dt, simfilen)
+
+    onewnc = NetCDFFile(stfilen, 'a')
+
+    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
+    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
+
+    for sfcv in sfcvariables.keys():
+        fvn = sfcvariables[sfcv]
+        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
+          not onewnc.variables.has_key(sfcv):
+            print gen.errormsg
+            print "  Newfile for the station '" + stfilen + "' does not content " +  \
+              " variable '" + sfcv + "' !!"
+            print '    variables:', onewnc.variables.keys()
+            quit(-1) 
+
+        if not gen.searchInlist(NONcheckingvars,fvn):
+            ogetvar = onc.variables[fvn]
+            getdims = ogetvar.dimensions
+        else:
+            getdims = ['Time', 'south_north', 'west_east']
+            getvdims = ['WRFtime', 'XLAT', 'XLONG']
+
+        # Looking for the X,Y axis for location of station within file
+        for dn in getdims:
+            axis = gen.dictionary_key_list(axesdims, dn)
+            if not dimvariables.has_key(dn):
+                print gen.errormsg
+                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
+                print '    add it to proceed'
+                quit(-1)
+            if axis == 'X':
+                lvals = dimvarvalues[dimvariables[dn]]
+                xdiff = (lvals-lonv)**2
+            elif axis == 'Y':
+                Lvals = dimvarvalues[dimvariables[dn]]
+                ydiff = (Lvals-latv)**2
+
+        # Looking for the closest point
+        diff = np.sqrt(xdiff + ydiff)
+        mindiff = np.min(diff)
+        minji = gen.index_mat(diff, mindiff)
+        
+        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
+        if np.any(coincstats == 0):
+            print '  Surface station coincidence by closest grid point from file !!'
+            iist = gen.index_vec(coincstats, 0)
+            print '  ' + labelv, '&', sfclabvals[iist]
+            onewnc.close()
+            sub.call('rm ' + stfilen, shell=True)
+            sub.call('cp ' + 'sfc_station_' + sfclabvals[iist] + '.nc ' + stfilen,   \
+              shell=True)
+            onewnc = NetCDFFile(stfilen, 'a')
+            ostlab = onewnc.variables['station']
+            Llab = len(labelv)
+            ostlab[0:Llab] = labelv
+            onewnc.sync()
+            onewnc.close()
+            if sfccoinc.has_key(sfclabvals[iist]):
+                lvals = sfccoinc[sfclabvals[iist]]
+                lvals.append(labelv)
+                sfccoinc[sfclabvals[iist]] = lvals
+            else:
+                sfccoinc[sfclabvals[iist]] = [labelv]
+            continue
+        
+        # Writting information to file
+        if sfcv == sfcrefvar:
+            ojvar = onewnc.variables['jpoint']
+            ojvar[:] = minji[0]
+            ojvar = onewnc.variables['ipoint']
+            ojvar[:] = minji[1]
+            onewnc.sync()
+
+        # Computing a single time the diagnosted variables
+        if gen.searchInlist(NONcheckingvars,fvn):
+            if ist == 0:
+                if fvn == 'WRFuas' or fvn == 'WRFvas':
+                    if not computeduasvas:
+                        computeduasvas = True
+                        ua, va, dmn, dmvn = diag.compute_WRFuasvas(                  \
+                          onc.variables['U10'][:], onc.variables['V10'][:],          \
+                          onc.variables['SINALPHA'][:], onc.variables['COSALPHA'][:],\
+                          getdims, getvdims)
+                        diagvars['WRFuas'] = ua
+                        diagvars['WRFvas'] = va
+                    ogetvar = diagvars[fvn]
+            else:
+                ogetvar = diagvars[fvn]
+                
+        slicevar = []
+        for dn in getdims:
+            axis = gen.dictionary_key_list(axesdims, dn)
+            if axis == 'X': slicevar.append(minji[1])
+            elif axis == 'Y': slicevar.append(minji[0])
+            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
+
+        # Writting data and slicing
+        onewvar = onewnc.variables[sfcv]
+        onewvar[:] = ogetvar[tuple(slicevar)]    
+        onewnc.sync()
+
+    if jumpstation: continue
+
+    onewnc.close()
+
+if len(sfccoinc.keys()) > 0:
+    print 'Coincident surface stations _______'
+    gen.printing_dictionary(sfccoinc)
+
+# Retrieving sounding data
+
+# Had already diagnostics been computed?
+computedta = False
+computedtd = False
+computeduava = False
+computedp = False
+computedzg = False
+
+# Diagnosted variables
+diagvars = {}
+
+# Coinident sounding station
+sndcoinc = {}
+
+for ist in range(Nsndstats):
+    jumpstation = False
+
+    lonv = sndlonvals[ist]
+    latv = sndlatvals[ist]
+    labelv = sndlabvals[ist]
+    heightv = sndhgtvals[ist]
+
+    stfilen = 'snd_station_' + labelv + '.nc'
+
+    creation_sndstation_file(stfilen, lonv, latv, heightv, labelv, utime, UnitsPress,\
+      presstime, sndvariables.keys(), dt, dz, simfilen)
+
+    onewnc = NetCDFFile(stfilen, 'a')
+
+    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
+    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
+
+    for sndv in sndvariables.keys():
+        fvn = sndvariables[sndv]
+        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
+          not onewnc.variables.has_key(sndv):
+            print gen.errormsg
+            print "  Newfile for the station '" + stfilen + "' does not content " +  \
+              " variable '" + sndv + "' !!"
+            print '    variables:', onewnc.variables.keys()
+            quit(-1) 
+
+        if not gen.searchInlist(NONcheckingvars,fvn):
+            ogetvar = onc.variables[fvn]
+            getdims = ogetvar.dimensions
+        else:
+            getdims = ['Time', 'bottom_top', 'south_north', 'west_east']
+            getvdims = ['WRFtime', 'WRFp', 'XLAT', 'XLONG']
+
+        # Looking for the X,Y axis for location of station within file
+        for dn in getdims:
+            axis = gen.dictionary_key_list(axesdims, dn)
+            if not dimvariables.has_key(dn):
+                print gen.errormsg
+                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
+                print '    add it to proceed'
+                quit(-1)
+            if axis == 'X':
+                lvals = dimvarvalues[dimvariables[dn]]
+                xdiff = (lvals-lonv)**2
+            elif axis == 'Y':
+                Lvals = dimvarvalues[dimvariables[dn]]
+                ydiff = (Lvals-latv)**2
+
+        # Looking for the closest point
+        diff = np.sqrt(xdiff + ydiff)
+        mindiff = np.min(diff)
+        minji = gen.index_mat(diff, mindiff)
+        
+        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
+        if np.any(coincstats == 0):
+            print '  Sounding station coincidence by closest grid point from file !!'
+            iist = gen.index_vec(coincstats, 0)
+            print '  ' + labelv, '&', sndlabvals[iist]
+            onewnc.close()
+            sub.call('rm ' + stfilen, shell=True)
+            sub.call('cp ' + 'snd_station_' + sndlabvals[iist] + '.nc ' + stfilen,   \
+              shell=True)
+            onewnc = NetCDFFile(stfilen, 'a')
+            ostlab = onewnc.variables['station']
+            Llab = len(labelv)
+            ostlab[0:Llab] = labelv
+            onewnc.sync()
+            onewnc.close()
+            if sndcoinc.has_key(sndlabvals[iist]):
+                lvals = sfccoinc[sndlabvals[iist]]
+                lvals.append(labelv)
+                sndcoinc[sndlabvals[iist]] = lvals
+            else:
+                sndcoinc[sndlabvals[iist]] = [labelv]
+            continue
+        
+        # Writting information to file
+        if sfcv == sfcrefvar:
+            ojvar = onewnc.variables['jpoint']
+            ojvar[:] = minji[0]
+            ojvar = onewnc.variables['ipoint']
+            ojvar[:] = minji[1]
+            onewnc.sync()
+
+        # Computing a single time the diagnosted variables
+        if gen.searchInlist(NONcheckingvars,fvn):
+            if ist == 0:
+                if fvn == 'WRFp':
+                    if not computedp:
+                        computedp = True
+                        diagvars['WRFp'] = onc.variables['P'][:] +                   \
+                          onc.variables['PB'][:]
+                elif fvn == 'WRFta':
+                    if not computedta:
+                        computedta = True
+                        ta, dmn, dmvn = diag.compute_WRFta(onc.variables['T'][:],    \
+                          onc.variables['P'][:]+onc.variables['PB'][:],              \
+                          getdims, getvdims)
+                        diagvars['WRFta'] = ta
+                elif fvn == 'WRFtd':
+                    if not computedtd:
+                        computedtd = True
+                        td, dmn, dmvn = diag.compute_WRFtd(onc.variables['T'][:],    \
+                          onc.variables['P'][:]+onc.variables['PB'][:],              \
+                          onc.variables['QVAPOR'][:], getdims, getvdims)
+                        diagvars['WRFtd'] = td
+                elif fvn == 'WRFua' or fvn == 'WRFva':
+                    if not computeduava:
+                        computeduasvas = True
+                        ua, va, dmn, dmvn = diag.compute_WRFuava(                    \
+                          onc.variables['U'][:], onc.variables['V'][:],              \
+                          onc.variables['SINALPHA'][:], onc.variables['COSALPHA'][:],\
+                          getdims, getvdims)
+                        diagvars['WRFua'] = ua
+                        diagvars['WRFva'] = va
+                elif fvn == 'WRFzg':
+                    if not computedzg:
+                        computedzg = True
+                        diagvars['WRFzg'] = onc.variables['PH'][:] +                 \
+                          onc.variables['PHB'][:]
+                ogetvar = diagvars[fvn]
+            else:
+                ogetvar = diagvars[fvn]
+                
+        slicevar = []
+        for dn in getdims:
+            axis = gen.dictionary_key_list(axesdims, dn)
+            if axis == 'X': slicevar.append(minji[1])
+            elif axis == 'Y': slicevar.append(minji[0])
+            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
+
+        # Writting data and slicing
+        onewvar = onewnc.variables[sndv]
+        onewvar[:] = ogetvar[tuple(slicevar)]    
+        onewnc.sync()
+
+    if jumpstation: continue
+
+    onewnc.close()
+
+if len(sndcoinc.keys()) > 0:
+    print 'Coincident sounding stations _______'
+    gen.printing_dictionary(sndcoinc)
+
