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

# Getting configuration
from config_get_stations import *

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

# Gneric variables prepared to be computed
Cvars = ['C_hur', 'C_hur_Uhus', 'C_td', 'C_td_Uhus', 'C_wd', 'C_ws']

# WRF specific diagnostics
WRFvars = ['WRFp', 'WRFta', 'WRFtd', 'WRFtime', 'WRFua', 'WRFuas', 'WRFva', 'WRFvas',\
  'WRFzg']

# Variables not to check their existence inside file
NONcheckingvars = Cvars + WRFvars + ['None']

def creation_sfcstation_file(filen, lv, Lv, hv, tv, 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'))
    newvar[:] = tv
    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, tv, 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'))
    newvar[:] = tv
    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
    #######

Wcomputed = {}
for vn in diag.Wavailablediags:
    Wcomputed[vn] = False

Ccomputed = {}
for vn in diag.Cavailablediags:
    Ccomputed[vn] = False

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 
if sfcstatfile != 'None':
    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
else:
    Nsfcstats = 0
    print '  No surface stations !!'

# sounding stations 
if sndstatfile != 'None':
    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
else:
    Nsndstats = 0
    print '  No sounding stations !!'


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)

    if not gen.searchInlist(NONcheckingvars, varn):
        # 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])

if sfcstatfile != 'None':
  dt = onc.variables[sfcrefvar].shape[0]
else:
  dt = onc.variables[sndrefvar].shape[0]

# Looking for 2D 'lon', 'lat' related variables
for xn in axesdims['X']:
    Xvarvals = dimvarvalues[dimvariables[xn]]
    Yvarvals = dimvarvalues[dimvariables[CFdims['lat']]]
    if len(Xvarvals.shape) == 1:
        newXvarvals, newYvarvals = np.meshgrid(Xvarvals, Yvarvals) 
        dimvarvalues[dimvariables[xn]] = newXvarvals
        if len(Yvarvals.shape) == 1:
            dimvarvalues[dimvariables[CFdims['lat']]] = newYvarvals

for yn in axesdims['Y']:
    Xvarvals = dimvarvalues[dimvariables[CFdims['lon']]]
    Yvarvals = dimvarvalues[dimvariables[yn]]
    if len(Yvarvals.shape) == 1:
        newXvarvals, newYvarvals = np.meshgrid(Xvarvals, Yvarvals) 
        dimvarvalues[dimvariables[yn]] = newYvarvals

tvals = dimvarvalues[dimvariables[CFdims['time']]]

# Retrieving surface data
##

# 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, tvals, 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 = list(onc.variables[sfcrefvar].dimensions)
            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
              dimvariables[getdims[2]]]

        # 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]
            ohvar = onewnc.variables['fheight']
            ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
            onewnc.sync()

        # Computing a single time the diagnosted variables
        if gen.searchInlist(NONcheckingvars,fvn):
            if ist == 0:
                if fvn[0:3] == 'WRF':
                    fvn_diag = fvn[3:len(fvn)]
                    if not Wcomputed[fvn_diag]:
                        # Using a given series of prepared diagnostics
                        Wcomputed[fvn_diag] = True
                        vardiag = diag.W_diagnostic(fvn_diag, onc,                  \
                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
                        diagvars[fvn] = vardiag.values
                elif fvn[0:2] == 'C_':
                    fvn_diag = fvn.split('_')[1]
                    if not Ccomputed[fvn_diag]:
                        # Using a given series of prepared diagnostics
                        Ccomputed[fvn_diag] = True
                        vardiag = C_diagnostic(fvn_diag, onc, sfcvariables,          \
                          sndvariables, CFdims)
                        diagvars[fvn] = vardiag.values
                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

# 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, tvals, 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 = list(onc.variables[sndrefvar].dimensions)
            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
              dimvariables[getdims[2]], dimvariables[getdims[3]]]

        # 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 sndv == sndrefvar:
            ojvar = onewnc.variables['jpoint']
            ojvar[:] = minji[0]
            ojvar = onewnc.variables['ipoint']
            ojvar[:] = minji[1]
            if dimvariables['H'] != 'None':
                ohvar = onewnc.variables['fheight']
                ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
            onewnc.sync()

        # Computing a single time the diagnosted variables
        if gen.searchInlist(NONcheckingvars,fvn):
            if ist == 0:
                if fvn[0:3] == 'WRF':
                    fvn_diag = fvn[3:len(fvn)]
                    if not Wcomputed[fvn_diag]:
                        # Using a given series of prepared diagnostics
                        Wcomputed[fvn_diag] = True
                        vardiag = diag.W_diagnostic(fvn_diag, onc,                   \
                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
                        diagvars[fvn] = vardiag.values
                elif fvn[0:2] == 'C_':
                    fvn_diag = fvn[2:len(fvn)]
                    if not Ccomputed[fvn_diag]:
                        # Using a given series of prepared diagnostics
                        Ccomputed[fvn_diag] = True
                        vardiag = diag.C_diagnostic(fvn_diag, onc, sfcvariables,     \
                          sndvariables, CFdims)
                        diagvars[fvn] = vardiag.values

                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)